Uropygial gland microbiota of nearctic−neotropical migrants vary with season and migration distance

Abstract

Symbiotic microbiota are important drivers of host behaviour, health, and fitness. While most studies focus on humans, model organisms, and domestic or economically important species, research investigating the role of host microbiota in wild populations is rapidly accumulating. Most studies focus on the gut microbiota; however, skin and other glandular microbiota also play an important role in shaping traits that may impact host fitness. The uropygial gland is an important source of chemical cues and harbours diverse microbes that could mediate chemical communication in birds, so determining the factors most important in shaping host microbiota should improve our understanding of microbially-mediated chemical communication. Hypothesizing that temporal, geographic, and taxonomic effects influence host microbiota, we evaluated the effects of season, migration distance, and taxonomy on the uropygial gland microbiota of 18 passerine species from 11 families. By sampling 473 birds at a single stopover location during spring and fall migration and using 16S rRNA sequencing, we demonstrate that season, followed by migration distance, had the strongest influence on uropygial gland microbial community composition. While statistically significant, taxonomic family and species had only weak effects on gland microbiota. Given that temporal effects on gland microbiota were nearly ubiquitous among the species we tested, determining the consequences of and mechanisms driving this seasonal variation are important next steps.

Background

Symbiotic microbes affect the behaviour, health, and fitness of their hosts [1, 2]. Foundational research investigating host-microbe co-evolution and the effects of symbiotic microbes on their hosts has focused on humans and model species [3,4,5], but our knowledge of wildlife microbial ecology is growing rapidly [6,7,8,9,10]. Given the importance of the gut microbiota for host health and behaviour, these communities are increasingly well studied in humans and animal models [11,12,13,14] as well as in wild animal populations [15,16,17,18,19]. More recently, the skin microbiome has been recognized as important in regulating host infection dynamics [20,21,22] and as a source of variation upon which selection may operate [23, 24]. By comparison, there is a paucity of data on other microbial reservoirs that can impact host behaviour, health, and fitness, such as the reproductive microbiome [25] and microbiota of the salivary glands [26, 27], mammary glands [28], tarsal glands [29], subcaudal gland [30, 31], anal glands [32, 33], and the uropygial gland [34,35,36].

The avian uropygial gland (preen gland) is a large holocrine integumentary gland located near the base of the tail in most bird species [37]. The uropygial gland is rich in oily secretions (uropygial or preen oil, often considered to be the main source of avian body odour [38, 39]) and harbors diverse microbiota [40, 41]. Birds frequently engage in ‘preening’ behaviour in which uropygial oil is regularly groomed throughout the feathers, and experimental evidence implicates microbes associated with the uropygial gland in the production or modification of chemical cues in uropygial oil that may be important for social communication in birds [42, 43].

Host microbiota are affected by multiple factors, including genetic [44,45,46], temporal [47,48,49], and spatial variation [24, 32, 41], and host age, sex [23, 29, 50], diet [51, 52], and evolutionary history [53,54,55,56]. In birds, the uropygial gland is an important source of odour cues involved in mediating both interspecific [57, 58] and intraspecific interactions [59,60,61]. Uropygial gland microbiota vary among populations [41], between captive and free-living birds [62, 63], with major histocompatibility complex (MHC) immune genotype [35, 64], and breeding group membership [65], but the consequences of such variation remains poorly understood. Determining the factors that drive variation in the uropygial gland microbiota is a key first step towards understanding the importance of microbially-mediated chemical communication in birds.

Comparative analyses can be instrumental in uncovering broad relationships that may reflect shared evolutionary history, selection pressures, or both. By looking at interspecies differences, comparative analyses can also inform us about species-specific adaptations, offering insights into the generation and maintenance of biological diversity [66]. Rapid technological advances have made microbial DNA sequencing more accessible and affordable than ever before. Recent comparative analyses have revealed incredible diversity in the gut microbiota of birds, providing evidence that host phylogeny and ecology (e.g., diet) explain some of the variation in gut microbiota [53, 54, 67,68,69] (but see [70]). Gut microbial community composition differs between resident and migratory avian subspecies [71], and the gut microbiome undergoes extensive remodelling during migration, presumably to meet the physiological and energetic demands of migration [72]. To determine whether similar patterns of variation are found for other avian microbial communities, large scale, multi-species sampling is needed.

We hypothesized that temporal, geographic, and taxonomic effects influence uropygial gland microbiota. To test this, we sampled bacteria from the uropygial gland of passerine birds at a single stopover site in Ontario, Canada during spring and fall migration. We used 16S rRNA gene sequencing to characterize the uropygial gland microbiota of 473 birds from 18 passerine species belonging to 11 families, then evaluated the effects of season (spring versus fall), migration distance (intercontinental versus intracontinental migration), and taxonomy (family and species) on microbial community composition and diversity. We predicted that all four factors would affect uropygial gland microbial community composition and diversity. Specifically, we predicted that intercontinental migrants would have greater microbial diversity than intracontinental migrants, particularly during spring migration, as these birds travel greater distances and likely encounter greater habitat diversity on the journey from the wintering grounds to the stopover site at which we sampled. We did not have directional predictions about the effects of season or taxonomic family and species on microbial community composition and diversity, but we expected to see broad differences in microbial community composition between seasons and among avian families and species.

Materials and methods

Study site and sample collection

All birds were captured at the Old Cut Field Station of the Long Point Bird Observatory (hereafter LPBO; 42°34′58.5″ N, 80°23′54″4 W) near Port Rowan, Ontario, Canada. We captured individuals using mist nets during fall (5 Sept–4 Nov 2020) and spring (9 Apr–26 May 2021) migration as part of the station’s standard, long term migration monitoring efforts at this location. We selected 18 species (order Passeriformes) for sampling, based on three main criteria: i) typical abundance at the study site during migration, ii) taxonomic diversity, and iii) selection of species pairs within each taxonomic family that differ in migration distance. We scored species as either intercontinental (long-distance, species that overwinter in South America) or intracontinental (short-distance, species that overwinter in North or Central America) migrants (Table 1). Ten of the 18 species studied breed in the mixedwood plains ecoregion, which includes the study site, with all but two of the species sampled (House Sparrow, Northern Cardinal) breeding as far north as the boreal shield. The intracontinental migrants overwinter as far south as Texas, Florida, Mexico, and Central America, with the shortest distance migrants (Black-capped Chickadee, House Sparrow) overwintering in the mid-eastern United States. All intercontinental migrants overwinter in South America.

Table 1 Species sampled for analysis and assessment of variation in uropygial gland microbiota according to season (fall, spring), migration distance, taxonomic family, and species

We attempted to collect uropygial gland swabs from 20 individuals of each species in each season (fall, spring) and to sample adult birds whenever possible. All birds sampled in spring were adults, since young are not yet born, but we typically sampled more juveniles than adults in fall because they outnumber adults at this time of year. We also attempted to balance the number of samples collected from each sex within species, but sex could not be determined in all cases. Full sampling details are available in Table S1.

All field technicians sanitized their hands using 70% ethanol prior to handling each bird, including prior to mist net extraction and again before handling for sample collection. Each captured bird was placed in a clean cotton bag, and all bags were washed between uses. One researcher (LAG) performed all sample collection. From each bird, we gently probed the uropygial gland with a capillary tube to express preen oil for use in a separate study. We then swabbed the gland by dipping a sterile medical grade swab into sterile phosphate buffered saline (PBS) then rubbed the swab around the gland using a continuous circular motion ten times each clockwise and counterclockwise. Given the small size of the gland in these species, it is impossible to noninvasively sample microbes from directly within the gland. However, by collecting preen oil first, our external swabbing method allowed us to collect microbes living immediately outside the gland as well as those inhabiting the gland that would have been excreted along with the oil [41, 62, 65]. Consistent with typical operations at the banding station, each bird was then banded, aged, sexed, weighed to the nearest 0.1 g, assigned a standardized score of fat reserves, and its wing chord length was measured to the nearest mm before release. Samples were stored on ice in the field for up to 6 h, then stored at −20 °C for 4–7 months (fall samples) or 2–3 months (spring samples) pending DNA extraction.

DNA extraction and 16S amplification

We extracted bacterial DNA from uropygial gland swabs using Norgen soil DNA isolation plus kits, consistent with prior work [41, 62, 65]. Extractions were carried out in batches of 24 or 48, with samples from different species distributed approximately equally among batches and a water-only negative control (a sterile swab dipped in PBS) in each batch. We amplified the V4 region of the bacterial 16S rRNA gene using the universal primers F518 [73] and R806 [74]. Each primer included an Illumina MiSeq adaptor sequence, four randomized nucleotides, and a unique ‘barcode’ of twelve nucleotides. We performed PCR in a total volume of 25 µL, with final concentrations of 1X buffer, 1.5 mM MgCl₂, 0.2 µM dNTPs, 0.2 µM of each primer, 1.2 U Platinum Taq Green Hot Start (Invitrogen), and 2 µL of DNA template. The thermocycling profile was 3 min at 94 °C; 30 cycles of 45 s at 94 °C, 60 s at 50 °C, and 90 s at 72 °C; and a 10 min final extension at 72 °C.

Sequencing and pipeline

We pooled PCR products of the expected size (approx. 300 nt) into a library and sequenced with 250 nt paired-end reads on an Illumina MiSeq at the London Regional Genomics Centre. Using the R package dada2 [75], we overlapped reads, removed ambiguous reads, chimeras, and singleton sequences, and assigned reads to samples. Sequences rarer than 0.1% in the full dataset were removed as they contain little information and removing them has no impact on downstream analyses [76]. The workflow and parameters used are available at github.com/ggloor/miseq_bin/. We obtained an initial dataset of 69 885 unique sequences (i.e., amplicon sequence variants; hereafter ASVs) from 574 samples and then assigned ASVs to taxon using the naïve Bayesian Ribosomal Database Project (RDP) Classifier [77].

Most of the ASVs in the initial dataset were rare. Rare sequences occurring in only a few samples are generally uninformative and samples with very low read counts likely represent undersampling, so we filtered sequences by the minimum proportion, minimum occurrence, and minimum sample count of reads. Sequences found in less than 0.5% of reads (consistent with MiSeq instrument error rates reported in [78]), fewer than 10% of samples, and samples with fewer than 5000 reads were removed (following [41, 62, 65]), resulting in the retention of 39.7% of reads and 473 samples (101 samples were removed due to low read counts). While this approach led to the removal of a large proportion of sequences, less stringent filtering conditions that removed sequences found in fewer than 5% and 1% of samples, respectively, led to the retention of a large number of sequences that were subsequently recommended for removal based on decontam ([79], see below), so we deemed filtering sequences found in fewer than 10% of samples to be an appropriate threshold. After removing 17 ASVs classified as organelles (chloroplast or mitochondria), we obtained a final data set of 134 ASVs (Table S2) from 473 samples from 18 species (x̄ retained reads per sample = 7 805 ± 189 SE, see Table S1 for sample sizes retained for each species after filtering).

We tested for signatures associated with external contamination in the full (unfiltered) dataset using the frequency and prevalence methods in the R package decontam [79]. The frequency method identifies contaminants by comparing the frequency distribution of each ASV as a function of input DNA concentration. In the contaminant model, the expected frequency varies inversely with total DNA concentration. In the non-contaminant model, the expected frequency is independent of the total DNA concentration [79]. Using the frequency method, we identified 9 of the initial 69 885 SVs as candidate contaminants, all of which were removed by our filtering steps (of these ASVs, 2 and 9 were retained by filtering sequences found in fewer than 5% and 1% of samples, respectively). The prevalence method identifies contaminants by comparing the presence/absence of each ASV in samples to the presence/absence of each ASV in negative controls and is appropriate for low biomass samples such as ours [79]. We did not identify any candidate contaminants using this method (at three different thresholds: 0.5, 0.1, 0.05). Thus, we retained all ASVs that passed our filtering and quality control steps.

High throughput sequencing generates relative abundance data that have a constant sum where the number of reads is imposed by the capacity of the sequencing instrument (and is therefore irrelevant to data interpretation), rendering the data compositional [80]. Compositional data provide information about the relationship among components [81], so we used a compositional data analysis approach that examines the read ratios between sequences [80, 82, 83]. In most such datasets, due to missing components, observed and actual totals are not equal. Small values like those below the detection limit of an instrument are typically observed as zero; however, in such cases these zero counts reflect sampling or equipment limitations rather than true zeros. While the value is certainly below some threshold, the true value is unknown (i.e., these are left-censored data). Discarding or replacing these values with zero can result in estimation bias, so we imputed the true values using Bayesian-multiplicative replacement to impute values for zero count sequences (following [76, 84]) using the R package zCompositions [85]. We then applied a centered log-ratio transformation to the zero-replaced data set, rendering the use of Euclidean (Aitchison) distances meaningful and straightforward for downstream analyses [82, 86].

Data analysis

We performed statistical analyses in R version 4.0.3 [87]. We conducted a PCA of the centred log-ratio transformed data using zero-centered rotated variables and the ‘prcomp’ function in base R. Based on visual analysis of the PCA scree plot and the cumulative variance explained by the principal components (PCs), we retained the first three PCs, which together accounted for 12.2% of the variance. Visual assessment of qq-plots and residuals indicated that the data and residuals were distributed approximately normally and the residuals showed no evidence of heteroscedasticity.

Analysis of high throughput sequencing data often combines multiple analytical approaches, and consistency of results across methods should increase confidence in the results [86]. Accordingly, we conducted analysis of variance (ANOVA) tests to evaluate the similarity of group means and permutational analysis of variance (PERMANOVA) tests to evaluate the similarity of group centroids and group dispersion. To test for differences among samples between seasons (fall, spring), with migration distance (intercontinental, intracontinental), taxonomic family, and species, we conducted three two-way ANOVAs using season, migration, family, and species as predictor variables and factor scores from each of the three retained PCs as dependent variables. We visualized the distribution of microbiota using a PCA biplot. We also conducted a PERMANOVA on the Euclidean distance matrix using the ‘adonis’ command in vegan [88] to test for an effect of season, migration distance, taxonomic family, and species on uropygial gland microbiota.

To control for any influence of our sampling design with respect to bird age and sex, we also tested for an effect of bird age and sex on uropygial gland microbiota using the subset of data for which we had sufficient information. We performed ANOVA tests on PC 1–3 to test for an effect of age on birds sampled in fall for five species: Backpoll Warbler, Eastern Wood-Pewee, Magnolia Warbler, Rose-breasted Grosbeak, and Swainson’s Thrush (Table S1). We performed ANOVA tests on PC 1–3 to test for an effect of sex on birds sampled in spring and fall for six species: Common Grackle, Golden-crowned Kinglet, House Sparrow, Magnolia Warbler, Northern Cardinal, and Rose-breasted Grosbeak, and on Baltimore Orioles sampled in spring (Table S1). For both sets of analyses, we pooled species and seasons (where we had data from both seasons) and tested only for a main effect of age or sex.

Next, to test for differentially abundant taxa between seasons and with migration distance (i.e., between groups), we conducted differential abundance tests on the raw read count data (retaining ASVs with an average of two counts across all samples) using t-tests and a generalized linear model with the ALDEx2 (v1.6.0) package in Bioconductor [86, 89, 90]. Next, we used coda4microbiome [91] to identify the maximally discriminant balances; groups of taxa that were maximally predictive of season. ALDEx2 tests each ASV, assuming independence of groups, identifying those that are maximally different between groups, while coda4microbiome tests groups of ASVs to identify subgroups that are maximally predictive of group membership.

To evaluate microbial diversity among samples, we calculated Shannon (alpha) diversity using the ‘diversity’ function in vegan [88]. Using diversity as the response variable, we then used a linear model to test whether season, migration, or species predicted microbial diversity.

To compare our compositional method with proportional methods commonly used in the ecology literature, we also analyzed the 16S rRNA gene sequencing data using another approach. We converted the raw read count data to proportions (rather than performing centered log-ratio transformations on the zero replaced data set), then conducted PERMANOVA on a Bray–Curtis distance matrix constructed from the proportional 16S rRNA gene read count data, again testing for an effect of season, migration distance, taxonomic family, and species on uropygial gland microbiota.

Finally, we conducted additional analyses of seasonal differences in uropygial gland microbiota separately for each species for which we had paired data (N = 16). We used the same PCA and ANOVA methods as described above for the full dataset, with season as the sole predictor variable, but we used two different filtering approaches. First, we used the full dataset that was filtered when all species were pooled (removing sequences found in less than 0.5% of reads, fewer than 10% of samples, and samples with fewer than 5000 reads), then split the data by species before running PCA and two-way ANOVA tests on the first two PCs. Second, we split the data by species and then filtered the dataset (removing sequences found in less than 0.5% of reads, fewer than 10% of samples, and samples with fewer than 5000 reads) before running PCA and two-way ANOVA tests on the first two PCs. We evaluated different filtering methods (i.e., filtering by sequences found in fewer than 5% and 1% of samples) for each species, but these less stringent filtering approaches again resulted in the retention of many sequences that were marked for removal by the package decontam. Across all species, our filtering approach resulted in the retention of 35–52% (x̄ = 40.7%) of reads and the first two PCs accounted for 15.4–25.4 (x̄ = 18.9%) of the variance.

Results

The 134 uropgyial gland amplicon sequence variants (ASVs) retained from across 18 passerine species were distributed among 6 bacterial phyla, 9 classes, and at least 23 orders, 40 families, and 74 genera (Fig. 1; Fig. S1; Table S2).

Fig. 1

Barplot showing bacterial taxa identified from uropygial glands of 16 neotropical passerine species. The relative abundances of ASVs that shared the same level of identification (e.g., multiple ASVs identified as belonging to bacterial species within the same genus) were pooled. Data were filtered to include only taxa with a relative abundance of at least 5% in at least one species within season; taxa with less than 5% abundance in any one species within season are pooled as ‘Other’. Taxa denoted F- are those identified to family level, all other taxa were identified to genus. Fig S1 shows all taxa in 1% or greater abundance in at least one species within season

Two-way ANOVA tests on the full dataset indicated significant differences in the uropygial gland microbial community composition of passerine birds between seasons (spring, fall; PC 1–3), with migration distance (intercontinental, intracontinental migrant; PC 1 and 3), taxonomic family (PC 1–3), and species (PC 3), with the largest effect being for season, followed by migration distance (Table 2). Component loadings of the first three principal components are available in the supplementary material (Table S3). Statistical results were generally supported by visual assessment of PCA biplots (Figs. 2 and 3; see Figs S2 and S3 for ASV loadings most strongly associated with each PC).

Table 2 Results of ANOVA tests using factor scores from the first three principal components of a PCA to test for differences in passerine uropgyial gland microbiota depending on season, migration distance, taxonomic family, and species (see Table S3 for factor loadings)

Fig. 2

PC1 and PC2 scores derived from relative abundances of passerine uropgyial gland bacterial amplicon sequence variants (ASVs) indicating variation between spring (green) and fall (gold) migration and with migration distance (circle: inter = intercontinental migrant, triangle: intra = intracontinental migrant). A version of this figure with loadings based on ASV relative abundances most strongly associated with PC1 and PC2 is available in the supplementary material (Fig. S2)

Fig. 3

PC2 and PC3 scores derived from relative abundances of passerine uropygial gland bacterial amplicon sequence variants (ASVs) indicating variation between spring (green) and fall (gold) migration and with migration distance (circle: inter = intercontinental migrant, triangle: intra = intracontinental migrant). A version of this figure with loadings based on ASV relative abundances most strongly associated with PC2 and PC3 is available in the supplementary material (Fig. S3)

A PERMANOVA test on the Euclidian distance matrix (compositional analysis) also identified significant differences in the uropgyial gland microbiota of passerine birds depending on season, migration distance, taxonomic family, and species (Table 3), and a PERMANOVA test on the Bray–Curtis distance matrix (proportional analysis) yielded similar results (Table S4).

Table 3 Results of PERMANOVA using a Euclidean distance matrix to test for differences in passerine uropygial gland microbiota depending on season, migration distance, taxonomic family, and species

Two-way ANOVA tests on the subset of five species sampled in fall for which we had a mix of adult and juvenile birds indicated no effect of bird age on uropygial gland microbiota (PC 1: F = 0.18, p = 0.67; PC 2: F = 0.03, p = 0.87, PC 3: F = 0.04, p = 0.84). Two-way ANOVA tests on the subset of six species sampled in both spring and fall, and one species sampled only in spring, for which we had sufficient data on the sex of sampled birds indicated no effect of sex on uropygial gland microbiota (PC 1: F = 0.45, p = 0.63; PC 2: F = 1.13, p = 0.32; PC 3: F = 1.17, p = 0.31).

Using the ALDEx2 general linear model, we detected only two differentially abundant microbial taxa between seasons, both in the family Beijerinckiaceae, with a false discovery rate cutoff of 5%; these were also recovered if the cutoff was dropped to 1%. Coda4microbiome identified eight taxa that collectively formed a maximally distinguishing log-contrast with four taxa in each season (spring, fall) (Fig. 4). The ratio between these groups of taxa was highly predictive in separating the seasons and had a mean cross-validated AUC of 0.90 (Fig. S4). A permutation test showed a mean of 0.51 (95% CI 0.44–0.56), indicating that the observed AUC was highly distinguishable from random permutations of the data. In contrast, there were no statistically significant ASVs observed for migration distance, nor were there any predictive log-contrast groups that separated by migration distance. Finally, there was no significant interaction between season and migration distance.

Fig. 4

Bacterial taxa composing the microbial signature that best differentiates among uropygial gland swabs collected during spring and fall migration based on coda4microbiome analysis. The magnitude of the coefficients represents the contribution of each variable to the model (gold = positive, green = negative). Where possible, ASVs were identified to genus, but two could only be identified to order or family (O-Enterobacteriales, F-Rhodocyclaceae)

The mean ± SE microbial Shannon (alpha) diversity of fall and spring samples was 3.3 ± 0.04 and 3.3 ± 0.03 respectively (t = 0.49, p = 0.62; Table S5 Fig. S5. Similarly, there was no difference in Shannon diversity among inter- and intracontinental migrants (t = −0.79, p = 0.43; Table S5 Fig. S5). Among the 18 species in this study, only the Hermit Thrush showed a significantly different (lower) Shannon diversity compared to other species (Table S5, Fig. S6).

Two-way ANOVA tests run for each species individually indicated significant differences in the uropgyial gland microbial community composition between seasons for 87.5% (14/16) of species (Fig. 5; Fig. S7; Table S6).

Fig. 5

Box and whisker plots of PC1 and PC2 scores derived for each species from relative abundances of passerine uropygial gland bacterial amplicon sequence variants (ASVs) indicating variation between spring and fall migration. Midline indicates median values, whiskers indicate minimum and maximum values. Asterisks denote significant differences at α = 0.05. See Table S3 for ANOVA output. Species photos courtesy of Brock and Sherri Fenton and used with permission. A version of this figure presented as PC biplots is available in the supplementary material (Fig. S5)

Discussion

Uropygial gland microbial community composition significantly differed between spring and fall migration in both our full model and in 87.5% (14/16) of within-species comparisons. In our full model, we also detected an effect of migration distance, with inter- and intracontinental migrants exhibiting significant differences in uropygial gland microbial community composition, and a weak, albeit significant, effect of taxonomic family and species on gland microbial community composition. Based on a subset of species for which we could assign age and sex, we found no effect of either variable on uropygial gland microbial community composition, consistent with prior studies [41, 65]. Shannon (alpha) diversity was similar between seasons, migration distances, and species, but Hermit Thrush had significantly lower alpha diversity than the other 17 species we evaluated.

Our ALDEx2 analysis identified two taxa in the family Beijerinckiaceae that were differentially abundant between spring and fall. Beijerinckiaceae includes generalist and specialist aerobic methanotrophs (methane metabolizers) [92, 93]. The source of these bacteria and their potential relationship to avian hosts is currently unknown. The microbial groups identified by coda4microbiome that best explained seasonal differences were the genera Rickettsiella and Actinomycetospora, order Enterobacteriales and family Rhodocyclacea (elevated in fall) and the genera Massilia, Buchnera, Limnobacter, and Bartonella (elevated in spring). Rickettsiella contains obligate intracellular bacteria species that are widespread invertebrate pathogens [94]. Species in the genera Actinomycetospora and Massilia are primarily associated with plants [95, 96], while species in the genus Buchnera are mutualistic intracellular symbionts of aphids [97], and Limnobacter species are thiosulfate-oxidizing bacteria predominantly associated with aquatic environments [98]. The Bartonella genus contains species that are facultative intracellular parasites, considered opportunistic pathogens of both human and non-human animals [99, 100]. While the plant and aquatic-associated bacteria (Actinomycetospora, Massilia, Limnobacter) are most likely environmentally derived, the potentially pathogenic bacteria Rickettsiella, Buchera, and Bartonella warrant further consideration, particularly in light of the potential role of migratory birds as disease vectors [101,102,103]. Given that Buchnera species are specialized aphid endosymbionts, these bacteria may also be most likely derived from the environment of migrating birds, potentially picked up during foraging. On the other hand, Rickettsiella and Bartonella contain species that include important human and non-human animal pathogens.

Bartonella species have been found in a variety of avian hosts, including Neotropical migrants such as the Purple Martin (Progne subis), Tree Swallow (Tachycineta bicolor), Eastern Bluebird (Sialia sialis), Northern Mockingbird (Mimus polyglottos), Red-winged Blackbird (Agelaius phoeniceus), Red-bellied Woodpecker (Melanerpes carolinus), Common Loon (Gavia immer), Ross's Goose (Anser rossii), and Lesser Snow Goose (Chen caerulescens caerulescens) and avian ectoparasites such as ticks, fleas, mites, and flies [104,105,106]. Recent evidence shows that Bartonella can jump between mammalian and avian hosts, mediated by generalist ectoparasite species and raising important concerns about their zoonotic potential [104, 105]. Bartonella have also been found in a broad range of avian families in the Pantanal region of Brazil [107] and Tomsk region of Russia [108], suggesting this genus is likely widespread across bird species.

Rickettsiella species have been found in a wide variety of invertebrates, including arthropods, arachnids, and crustaceans [94, 109, 110]. This genus has recently been identified in wild populations of Atlantic salmon (Salmo salar) and samples collected from the Salton Sea, suggesting the prevalence of Rickettsiella, typically considered an invertebrate specialist, in different environments and hosts has likely been underestimated [94]. Rickettsiella species have also been identified in the poultry red mite (Dermanyssus gallinae), which is endemic in many commercial poultry farms and has significant impacts on hen health, reducing both egg production and quality [111]. Bacteria in the family Diplorickettsiaceae, which contains the genus Rickettsiella, were recently detected on the feathers, but not in the uropygial gland, of a New Guinean bird, the Regent Whistler (Pachycephala schlegelii) [112]. To the best of our knowledge, it is currently unknown whether Rickettsiella can be pathogenic in vertebrate hosts.

Seasonal variation in vertebrate gut microbiota has frequently been demonstrated in the literature, commonly associated with changes in host diet [48, 113,114,115,116], hibernation physiology [117], and migration physiology [72]. For example, Swainson’s Thrush and Gray Catbird (Dumetella carolinensis, not included in our study) exhibit seasonal shifts in gut microbiota during spring and fall migration, but sampling location was not held constant between the two time points in this study, so geographic variation could at least partially explain the observed variation [118]. A recent study of four Neotropical migrant thrush species, including Swainson’s and Hermit Thrush, found that sampling year had the greatest effect on intestinal microbial community composition and diversity, followed by season (diversity) and species (community composition), with no significant effect of age or sex on either measure [119]. In this study, birds sampled in summer exhibited higher microbial diversity than those sampled in fall and spring [119], which could be consistent with our finding of similar diversity between spring and fall across 16 passerine species.

Seasonal variation in skin microbiota has also been documented, particularly in amphibians, and may be related to changes in diet, habitat type, temperature, and susceptibility to infection [47, 120, 121]. While few studies have focused on the avian uropygial gland, seasonal variation in gland microbiota has been shown in the Eurasian Hoopoe (Upupa epops), with bacterial density increasing from the nonbreeding to breeding season in females [122], suggesting a possible role in antimicrobial defense during nesting when the likelihood of bacterial infections is highest [43, 122]. Our study took place during spring and fall migration, and did not reveal any differences in alpha diversity between the pre- and post-breeding periods. More detailed studies measuring shifts in the uropygial gland microbiota throughout the annual cycle may reveal patterns of microbial diversity in other species consistent with findings in Eurasian Hoopoe, providing insight into the potential benefits of host-associated microbiota.

We detected an overall effect of season (spring versus fall) on the uropygial gland microbiota of passerine birds sampled at a single stopover site in Ontario, Canada, with almost all species (14/16) sampled exhibiting seasonal changes in community composition but not diversity. Although our retained PCs (1–3) explained only 12.2% of the variance in our full model, we conclude that the temporal effect on uropygial gland microbiota is likely a general trend across passerine birds. Indeed, when evaluating seasonal effects on individual species, season explained an average of 18.9% (range 15.4–25.4%) of the variance in the data. In one study investigating the uropygial gland microbial communities of Dark-eyed Juncos, individual identity explained nearly 50% of variation in the microbial profiles of sampled birds [42]; thus, individual differences may similarly account for the relatively low variance explained by our variables of interest, particularly when pooling data across many taxonomically diverse species.

Given that seasonal shifts in diet are associated with seasonal variation in gut and skin microbiota, it is tempting to assume that seasonal shifts in uropygial gland microbiota may also be due to dietary shifts. Diet can affect the composition of uropygial secretions (preen oil) [123], which may in turn affect the uropygial gland microbiota [42, 43] (but see [41, 62]). We did not analyze the diet of our study species directly; however, based on data available online [124], we sampled a mix of insectivores, granivores, and omnivores, and most exhibit seasonal shifts in diet, typically either increasing consumption of fruit or other plant material during fall migration, or shifting the invertebrate taxa consumed throughout the annual cycle (Table S7) [124]. Data on dietary variation are lacking for three of our study species that exhibited seasonal shifts in gland microbiota. Of the two species that did not exhibit seasonal shifts, no diet data were available for the Eastern Wood-Pewee, but the Blue-headed Vireo does exhibit seasonal shifts in diet (Table S7). In a study of captive versus free-living Song Sparrows (Melospiza melodia), dietary treatments administered to captive birds did not affect uropygial gland microbial community composition or diversity [62]. We are thus unable to link seasonal changes in uropygial gland microbiota to seasonal changes in diet, but this warrants further study.

Seasonal variation in the chemical composition of preen oil secreted from the uropygial gland is nearly ubiquitous [125], and microbes likely influence preen oil composition [42, 43]. However, the reverse may also be true; preen oil chemistry could influence uropygial gland microbiota. Studies have failed to find an overall relationship between uropygial gland microbiota and preen oil chemical composition, but such relationships could be masked by the complexity, diversity, and multifunctionality of both uropygial gland microbial and preen oil chemical communities [35, 126]. The bacterial genera Pseudomonas and Staphylococcus have been associated with production of preen oil chemicals that are likely involved in intraspecific sexual communication in Dark-eyed Juncos, a Neotropical migratory passerine bird [42]. These genera were detected in all species included in our study, with relative abundances ranging from 1 to 8.7% (Pseudomonas) and 1.4–24.1% (Staphylococcus) (Fig. 1). However, based on our analyses, neither genus is likely to be driving the seasonal differences we observed. The most abundant genera found across species in our study were Rickettsiella (3.8–54.5%) and an unresolved ASV placed in the genera Burkholderia-Caballeronia-Paraburkholderia (4.8–13.8%) (Fig. 1). These genera contain species known to produce volatile compounds [127] and thus have the potential to be involved in avian chemical communication. More work is needed to evaluate the drivers and consequences of seasonal variation in uropygial gland microbiota and their relationship with seasonal variation in preen oil chemical composition.

We detected an effect of migration distance, broadly categorized as intercontinental (long distance) or intracontinental (short distance) migration, on uropygial gland microbial community composition but not diversity. While the mechanisms by which uropygial gland microbiota might be affected by migration patterns have rarely been investigated, we reasoned that seasonal shits in diet, physiology, and geography, known to affect avian fecal and gut microbiota [16, 41, 72], are likely involved given that these factors can also affect the chemistry of uropygial secretions [125], which might in turn affect uropygial gland microbiota. Environmental exposure during migration may be limited depending on species-specific movement patterns. For example, species may use the same stopover, wintering, and breeding sites annually, thereby limiting their microbial exposure to predictable locations or environments [18]. We compared long and short distance migrants, reasoning that species travelling further would necessarily have more opportunities to encounter microbes than species travelling shorter distances. However, we did not find support for our prediction that microbial diversity would be higher in long distance migrants. Importantly, our characterization of migration distance was broad, and does not take into account the stopover frequency or duration of the species we sampled, which could exert a greater influence over host microbiota than the absolute distance travelled. Nevertheless, we found a significant influence of migration distance on uropygial gland microbial community composition, likely driven by environmental factors associated with migration.

A potential explanation for our finding is that species wintering farther south (i.e., in South America) encounter a different suite of microbiota than species wintering farther north (i.e., in North America), despite encountering a similar diversity of microbes overall. To test this, species should be sampled on both the wintering and breeding grounds. For example, in Kirtland’s Warblers (Setophaga kirtlandii), gut microbial community composition and diversity differed between breeding and wintering sites, and alpha diversity and relative abundance changed over time within a single sampling location [128], suggesting that microbial turnover is rapid and dynamic. While this result conflicts with our suggestion, it is important to note that the factors affecting gut and uropygial gland microbiota may differ. A major advantage of our approach is that by sampling a wide variety of species at the same location at two different time points, we detect broad trends that likely apply across passerine species, providing a valuable starting point for future research.

Microbial turnover is poorly studied in birds [18], so it is unknown whether the microbial communities we detected during spring and fall migration are more likely reflective of local differences at the stopover site we sampled at, or whether they reflect microbiota associated predominantly with the breeding (fall samples) or wintering (spring samples) grounds. Conducting repeated sampling of birds throughout their migratory range could resolve some of these uncertainties [128], but would require a coordinated research effort across international borders and careful study design to reduce confounders. An alternative approach could be to hold birds captive under semi-natural outdoor conditions (i.e., in aviaries), offering seasonally appropriate wild-type diets that may support the maintenance of wild-type microbiota [17]. Repeated sampling throughout the annual cycle, with effects of geography removed, would then allow us to evaluate the effects of host physiology and dietary shifts on microbial turnover and seasonal variation in microbiota.

Conclusions

The uropygial gland microbial community composition, but not diversity, of passerine birds was best explained by temporal (seasonal) variation, followed by migration distance. We also detected a weak but significant effect of taxonomic family and species on microbial community composition. Temporal effects on the uropygial gland microbiota were nearly ubiquitous across the species we tested and should be taken into account when designing sampling protocols. Future research should focus on determining the mechanisms driving seasonal variation in host microbiota and evaluate the consequences of this variation on host fitness.