Experimental evidence that symbiotic bacteria produce chemical cues in a songbird

Animals are inhabited by microbial communities that can provide services for their hosts, influencing nutrient uptake, development, immune response and behaviour. Among the most frequently studied links between animal behaviour and microbiota is microbial involvement in host intraspecific communication, particularly chemical signalling (Archie and Theis, 2011; Davis et al., 2013; Ezenwa and Williams, 2014). Many vertebrate integumentary glands, especially specialized scent glands, are warm, moist places that are conducive to microbial growth. The fermentation hypothesis for chemical recognition (Albone et al., 1974; Gorman et al., 1974) posits that bacteria in scent glands produce volatile compounds that can be used as recognition cues by the host animals, and that variation in these resident bacterial communities can translate to variation in chemical recognition cues among hosts, and, through host–microbiota ritualization, potentially to reliable chemical signals. Researchers have demonstrated correlations between scent gland bacterial community structure and host animal odour profiles, including co-variation with host sex, age and social group membership (Leclaire et al., 2014a; Theis et al., 2012, 2013), but no studies have yet systematically demonstrated, through experiment, that variation in bacterial community structure translates to variation in host odour profile.

Chemical odourants produced by symbiotic bacteria may provide information about an animal host's species, sex, group membership, individual identity or condition (Archie and Theis, 2011; Davis et al., 2013). Group- or sex-specific behaviours affecting the acquisition and transmission of microbes can ultimately affect microbial community composition and, consequently, the host animal's odour profile, thereby providing reliable information about group membership or sex (Leclaire et al., 2017a, 2014a; Theis et al., 2012, 2013). Host genotype, particularly immune system genes such as the major histocompatibility complex (MHC), can regulate microbiota composition via the tolerance or rejection of particular microbes (Bolnick et al., 2014; Hooper et al., 2012; Kubinak et al., 2015). Chemical cues reflect MHC genotype in many vertebrates (Penn and Potts, 1999; Setchell et al., 2011), including birds (Leclaire et al., 2017b; Slade et al., 2016). Because MHC molecules can exclude bacterial species from, and thus directly shape, the host's microbiota (Penn, 2002), odour-producing symbiotic microbes are a strong candidate mechanism to explain the relationship between MHC genotype and chemical signal structure.

Although mammals are the best-studied vertebrates in this area of research (Leclaire et al., 2017a; Sin et al., 2012; Theis et al., 2013), recent research has also focused on the avian uropygial gland. Birds spread preen oil (comprised largely of aliphatic monoester waxes and fatty acids; Kolattukudy and Rogers, 1978) secreted by this gland on their feathers while preening, where it provides a variety of functions (Moreno-Rueda, 2017), including maintaining feather condition, coloration and brightness (López-Rull et al., 2010; Piersma et al., 1999), and protecting against ectoparasites (Clayton et al., 2010) and pathogenic bacteria (Shawkey et al., 2003). Like secretions from mammalian scent glands, preen oil gives off volatile compounds, which function as intraspecific chemical signals, indicating species, sex and breeding condition (Mardon et al., 2010; Whittaker et al., 2010, 2011b). Potential mates and rivals are attracted to preen oil odours (Amo et al., 2012; Whittaker et al., 2011a; Zhang et al., 2010), and birds may use them for mate and kin recognition (Bonadonna and Nevitt, 2004; Krause et al., 2012). Relative abundances of particular preen oil volatile compounds have been shown to predict reproductive success in dark-eyed juncos (Junco hyemalis) (Whittaker et al., 2013). Additionally, similarity among individuals' MHC genotype is reflected in the similarity of their preen oil odour profiles (Leclaire et al., 2017b, 2014b; Slade et al., 2016), suggesting that these profiles may be a reliable indicator of genotypes potentially important in mate choice.

Symbiotic bacteria are associated with the uropygial gland opening, and may contribute to preen oil's beneficial qualities. For example, Enterococcus bacteria in hoopoe (Upupa epops) uropygial glands produce bactericidal compounds that help protect against pathogenic infections (Ruiz-Rodríguez et al., 2013). Recently, bacteria have been hypothesized to produce the volatile compounds found in preen oil (Whittaker and Theis, 2016). These compounds include linear alcohols, methyl ketones and carboxylic acids, many of which can be end products of bacterial metabolism (Madigan et al., 2010). A survey of bacterial communities associated with dark-eyed junco uropygial glands identified several odour-producing genera capable of producing the same volatile compounds that are prominent in junco preen oil (Whittaker and Theis, 2016). Microbial communities and volatile profiles reflected group membership in juncos, but significant covariation between microbiota and volatile profiles could not be detected, perhaps owing to high microbial community diversity and/or the potential for uropygial gland bacteria to provide services other than chemical signalling (Whittaker et al., 2016). Thus, we must move beyond correlative studies, and conduct experiments to understand how symbiotic microbes contribute to avian chemical cues and signals.

Here, we first addressed this question by administering a broad-spectrum antibiotic directly to the uropygial glands of captive dark-eyed juncos and characterized the effects on uropygial microbiota and preen oil volatile profiles. Next, we collected preen oil from free-living juncos and cultivated bacteria from the oil. We found that all of the cultivated bacteria produced volatile compounds known to be important in junco reproductive behaviour. Together, these results suggest a significant role for symbiotic bacteria in the production of junco chemical cues and signals.

We captured 20 (8 male, 12 female) adult northern dark-eyed juncos (Junco hyemalis hyemalis Linneaus 1758) in Bloomington and Muncie, IN, USA, in February and March 2013, and housed them at the Kent Farm Bird Observatory in Bloomington, IN, USA. The experiment described for this study took place in July 2013, when birds were in full breeding condition. We determined individual sex using wing length and the presence/absence of a brood patch, and later confirmed via molecular sexing (Griffiths et al., 1998). The birds were housed in individual cages in four indoor aviary rooms with light cycles mimicking the natural light:dark cycle for the date. Birds were fed ad libitum a diet of millet, sunflower seeds, orange slices, and a mixture of puppy chow, hardboiled eggs and carrots.

We randomly assigned birds to control or experimental groups (N=4 males and 6 females in each group). Prior to administering treatment on day 1, we took swab and preen oil samples from each bird. Preen oil was collected by gently rubbing the uropygial gland with a 100 µl glass capillary tube (Drummond Scientific), which stimulates the gland to secrete 1–3 mg of oil (Whittaker et al., 2010). We did not sterilize the outer surface of the gland before sampling; thus the samples also contained biological material from the outside of the gland (including older preen oil and any bacteria present). We sampled bacterial communities by rubbing the tip of the uropygial gland using a sterile cotton swab, ensuring that our samples included preen oil and microbes from inside and outside the gland, which represent the mixture that birds collect on their bills in preparation for preening. We stored all samples at −80°C until analysis.

On days 1–5, we injected a broad-spectrum antibiotic, Enrofloxacin (Baytril), directly into the uropygial gland of birds in the experimental group once per day (five times total) (Martín-Vivaldi et al., 2010). Using a 26-gauge needle and a syringe, we injected 0.01 ml of Baytril 2.5% injectable solution (25 mg ml⁻¹) into each of the two uropygial gland lobes, following the recommended dose of 20 mg kg⁻¹ (Carpenter, 2012). Birds in the control group were injected with sterile saline. Post-treatment preen oil and uropygial gland swab samples were collected on day 7.

We extracted volatile compounds from preen oil samples using stir bar sorptive extraction (Soini et al., 2005): we placed the capillary tube containing the sample into a vial containing 2 ml of high-purity water, 100 mg of ammonium sulfate and a Twister^® (Gerstel GmbH, Mülheim an der Ruhr, Germany) stir bar. An internal standard (8 ng of 7-tridecanone in 5 μl of methanol) was added to each sample. The samples were stirred at ≥800 rpm for a 60-min extraction time. We performed quantitative analysis with an Agilent 6890N gas chromatograph connected to a 5973i MSD mass spectrometer (Agilent Technologies, Inc., Wilmington, DE, USA) with a Thermal Desorption Autosampler and Cooled Injection System (TSDA-CIS 4 from Gerstel) (Whittaker et al., 2010). We identified all major compounds by comparison to standards (Sigma-Aldrich), using mass spectra and retention times. We normalized peak areas of the compounds of interest by dividing each peak area by that of the internal standard (7-tridecanone) in corresponding runs, yielding relative concentrations (i.e. relative amounts per 1.0 mg of preen oil). Relative standard deviation (RSD, a measure of reproducibility) of the internal standard peak area was 11% (N=19).

We focused on junco preen oil volatile compounds that significantly increase in birds in breeding condition (Soini et al., 2007), as these compounds are most likely to play a role in sexual signalling and mating behaviour. Because the abundance measurements were not normally distributed, we used a natural log transformation, resulting in normal distributions for all but four variables (pre-treatment 1-tridecanol and dodecanoic acid, and post-treatment 1-undecanol and dodecanoic acid). Data for each volatile compound before and after treatment are shown in Table S1 and Fig. S1. To counteract potential problems associated with correlated variables and multiple comparisons, we used a principal components analysis (PCA) with varimax rotation to reduce the number of variables, retaining PCs with an eigenvalue greater than 1. Three PCs explained a cumulative 93.5% of the variance (Table 1). PC1 was most heavily influenced by nine linear alcohols, PC2 corresponded to three methyl ketones and PC3 corresponded to three carboxylic acids. Using one-way ANOVA, we observed no significant pre-treatment sex differences in any of the three PCs (P>0.05, all comparisons), so we grouped the sexes together for subsequent analyses. Nevertheless, we included sex as a variable in our models.

To analyse the effects of the antibiotic treatment on the volatile compound PC scores, we fit a linear mixed-effects model with random intercepts using the lme function of the R package nlme (https://CRAN.R-project.org/package=nlme). Fixed effects included in the model were sex, treatment (antibiotic or saline), time (pre- or post-treatment), and the interaction between time and treatment. We included bird ID as a random effect. We implemented these analyses in R 3.5 (https://www.r-project.org/) using RStudio 1.1.447 (https://rstudio.com/). We also conducted additional analyses with random slopes (Schielzeth and Forstmeier, 2009) and additional interaction terms (sex×treatment and sex×time), but the results were not substantively different from what we report here.

We extracted genomic DNA from swabs using the MO BIO Powerlyzer PowerSoil^® DNA Isolation kit (MO BIO Laboratories, Inc.), with minor modifications to the manufacturer's protocol: (1) swabs were immersed in 500 µl of bead solution and 200 µl of phenol:chloroform:isoamyl alcohol for 10 min prior to bead beating; (2) 100 µl of solution C2, 100 µl of solution C3 and 1 µl of RNase A were added and samples were incubated at 4°C for 5 min prior to centrifugation rather than being added to samples over two centrifugation steps; (3) lysates were combined with 650 µl of solution C4 and 650 µl of 100% ethanol prior to loading samples onto the spin column rather than being combined with 1200 µl of solution C4 alone; and (4) DNA was eluted in 60 µl, rather than 100 µl, of solution C6. The two extractions for each bird (i.e. pre- and post-treatment) were completed in the same batch. To assess any potential issues with background DNA contamination, we also conducted six extractions of blank DNA extraction kits. These background technical controls were processed and sequenced alongside the junco samples.

Owing to the low microbial biomass and high concentration of background host DNA in these swab samples, we amplified and sequenced the bacterial DNA using a nested PCR approach (Fan et al., 2009; Yu et al., 2015). The first round of amplifications included 15 cycles, with each reaction containing 1 µl each of the universal 16S rRNA gene primers 27f-CM (5′-AGA GTT TGA TCM TGG CTC AG-3′) and 1492R (5′-ACG GCT ACC TTG TTA CGA CTT-3′), 12.5 μl of 2X GoTaq Green Master Mix (Promega, Madison, WI, USA), and 3.0 μl purified DNA. Following an initial 5 min incubation at 95°C, the cycling parameters were 94°C for 30 s, 50°C for 30 s and 72°C for 120 s. Amplification products were diluted in nuclease-free water by a factor of 1:15, and submitted to the University of Michigan's Center for Microbial Systems (Ann Arbor, MI, USA) for sequencing. We sequenced the V4 region of the 16S rRNA gene on the Illumina MiSeq platform using a V2 500 cycle MiSeq Reagent Kit (Illumina MS102-2003), and the dual indexing sequencing strategy developed by Kozich et al. (2013). The primers targeting the V4 region of the 16S rRNA gene were 515F (5′-GTG CCA GCM GCC GCG GTA A-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′). PCR was performed under the following conditions: 95°C for 2 min, followed by 30 cycles at 95°C for 20 s, 55°C for 30 s and 72°C for 5 min, with an additional elongation at 72°C for 10 min. MiSeq sequence files for all samples have been deposited on the NCBI Sequence Read Archive (BioProject ID PRJNA554313, accession numbers SAMN12262297–SAMN12262336: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA554313).

Sequences were processed using mothur software (Schloss et al., 2009). Reads were assembled, trimmed, filtered (no ambiguous base calls and no homopolymer runs exceeding eight bases), aligned (SILVA reference database, release 102; Pruesse et al., 2007; Quast et al., 2013), screened for chimeras (UCHIME; Edgar et al., 2011), taxonomically classified (SILVA; confidence threshold 80%; sequences classified as Archaea, chloroplasts or mitochondria were removed) and clustered into operational taxonomic units (OTUs) based on a percent nucleotide similarity of 97%. The average length of all bacterial sequences in the final dataset was 253 bases. Two male juncos were excluded from all microbiota analyses because their pre-antibiotic treatment samples did not yield sequence data.

Collectively, the junco uropygial gland and background technical control samples (i.e. blank DNA extraction kits) yielded 839,064 quality sequences. These sequences clustered into 327 OTUs, including 113 singletons and 44 doubletons. The remaining junco samples had at least 7285 sequences and Good's coverage values exceeding 99.9%.

Of the six technical controls addressing potential background DNA contamination, four yielded at least 1827 sequences with Good's coverage values exceeding 99.6%. The two other technical control samples were excluded from analyses. The bacterial profiles of the 36 samples from the 18 juncos (12 female, 6 male) were distinct from those of the background technical controls [non-parametric multivariate ANOVA (NPMANOVA: Bray–Curtis; F_1,38=1.460, P=0.009; Fig. 2A]. Nevertheless, the R package decontam (Davis et al., 2018) was used to statistically identify OTUs that were potential contaminants based on their pattern of occurrence in biological versus technical control samples. Using a threshold of 0.4, 21/327 (6.4%) OTUs in the dataset were identified as likely contaminants. The taxonomic characterizations of these OTUs were Actinobacteria (1 OTU; uncl. Solirubrobacterales), Bacteroidetes (2 OTUs; Bacteroides, Chryseobacterium), Firmicutes (6 OTUs; Abiotrophia, Blautia, Caldalkalibacillus, Gemella, Lactobacillus, Paenibacillus) and Proteobacteria (12 OTUs; uncl. Proteobacteria, uncl. Comamonadaceae, uncl. Sphingomonadaceae, Achromobacter, Bradyrhizobium, Caulobacter, Escherichia, MN_122.2a, Novosphingobium, Phyllobacterium, Sphingomonas, Suttonella). These OTUs were removed from the dataset prior to microbiota analyses.

After removal of these 21 OTUs, all but one of the junco samples had at least 7221 sequences and Good's coverage values exceeding 99.9%. The exception was a pre-saline treatment sample from a female junco (32 sequences). This female (ID 239121749) was thus excluded from all microbiota analyses.

For alpha diversity analyses (i.e. diversity within samples), we subsampled each sample's dataset to 7221 sequences, and calculated OTU richness and Simpson evenness (1–D). Statistical comparisons were made using t-tests, Mann–Whitney tests, paired t-tests and Wilcoxon tests, as appropriate. Sequential Bonferroni corrections were applied to account for multiple comparisons. For beta diversity analyses (i.e. diversity between samples), we removed singleton and doubleton OTUs from the dataset and converted each sample's OTU count data to percent relative abundance data. Beta diversity between sample microbial profiles was characterized using the Bray–Curtis index. Alpha diversity indices and matrices of beta diversity similarities were calculated using PAST software (v 3.14, Hammer et al., 2001). We visualized beta diversity among samples through principal coordinates analyses (PCoA) in PAST (Hammer et al., 2001), and statistically evaluated variation in beta diversity through NPMANOVA (Anderson, 2001) using the adonis function (999, or the complete set of available, permutations, as is applicable) within the vegan package (v 2.5-4) in R (https://CRAN.R-project.org/package=vegan). The strata command within the adonis function was used to control for individual identity when assessing treatment effects within a repeated-measures framework. As previously observed (Whittaker et al., 2016), there was not a consistent difference in the pre-treatment bacterial profiles of female (N=11) and male (N=6) uropygial glands (NPMANOVA: F_1,15=0.996, P=0.428). Furthermore, within the saline-treated group following treatment, there was also no consistent effect of sex on glandular bacterial profile (N=5 females and 4 males; F_1,7=0.767, P=0.903). Therefore, sex was not included as a factor in the microbiota analyses. We visualized the relative abundance of the most prominent OTUs (i.e. those OTUS with >1% average relative abundance) in post-treatment uropygial gland samples by creating a heatmap using the heatmap.2 function in the gplots package in R (Warnes et al., 2016). Similarity percentage (SIMPER) analysis was used to elucidate the contribution of specific OTUs to the overall observed variation in uropygial gland bacterial profiles between saline and antibiotic treatment groups post-treatment (Clarke, 1993).

To assess which bacterial OTUs best explained variation in the volatile profiles of preen oil, we used information theoretic comparisons of models containing all possible combinations of predictive variables. For each volatile principal component, we developed a global model including sex, treatment, time, the interaction between time and treatment, and abundance of the most prominent (>1% average relative abundance) OTUs. This list of OTUs included those that were most prominent in both the pre- and post-treatment samples, as well as two that were only prominent in post-treatment samples (OTU 12, unclassified Sphingomonadaceae, and OTU 19, Corynebacterium) and one that was only prominent in pre-treatment samples (OTU 11, Stenotrophomonas). To facilitate the interpretation of model estimates, we standardized the input variables using the standardize function of the R package arm (Gelman, 2008; Grueber et al., 2011). For each global model, we plotted the residual and fitted values, and Q-Q norm plots to check model assumptions. We then used the dredge function of R package MuMIn (https://CRAN.R-project.org/package=MuMIn) to generate comparison models, including the intercept-only null model, and defined the set of top models as those within two AICc units of the best-fit model. Finally, we model-averaged these top models using the natural average method (Burnham and Anderson, 2002). We report the conditional average results for each parameter, which only averages over the models where the parameter appears, and is appropriate when variables may have a weak effect obscured by the inclusion of zeroes in the model average (Grueber et al., 2011).

We captured 11 adult Carolina dark-eyed juncos (J. hyemalis carolinensis Brewster 1886) (7 male, 4 female) in baited mist nets at Mountain Lake Biological Station near Pembroke, VA, USA, in August 2017. This subspecies is closely related and ecologically very similar to the northern dark-eyed junco population used in experiment 1, and the two subspecies have an overlapping distribution.

We collected two preen oil samples from each subject using the same methods as experiment 1. The ends of the tubes containing the preen oil were placed into transport media, wherein they remained until plating. We stored one sample in an Anaerobic Transport Medium tube (Anaerobe Systems) and the other in an Amies Transport Medium tube (Thermo Scientific). Samples were stored in a cool dark drawer at room temperature until they were delivered to Wayne State University, Detroit, MI, USA, for cultivation (range: 29–78 h between collection and plating).

Preen oil samples were plated on tryptic soy agar (TSA) medium and incubated under both aerobic and anaerobic (95% N, 5% H) conditions at 37°C for 2 weeks. The preen oil of nine of 11 juncos yielded growth on TSA plates under aerobic conditions. Two of these samples also yielded growth under anaerobic conditions. The plates with growth bore fewer than 10 total colonies, typically with a single morphologically distinct cultivar. Colonies of visually distinct morphotypes were streaked for isolation in pure culture. If plates did not yield bacterial cultures within 2 weeks, they were discarded. Control TSA plates in both conditions yielded no cultivars.

A single isolated colony from each pure-culture plate was placed into nuclease-free water, mixed and used directly in colony polymerase chain reaction (PCR). Each sample was amplified in 25 µl reactions containing 1 µl of template, 1 µl each of 27F and 1492R bacterial 16S rRNA gene primers at 10 µmol l⁻¹, 9.5 µl of nuclease-free water, and 12.5 µl of SybrGreen PCR Supermix. PCR was performed per the following conditions: 95°C for 5 min, followed by 32 cycles of 95°C for 30 s, 46°C for 30 s and 72°C for 90 s, with an additional 5 min elongation at 72°C after the last cycle. PCR products were submitted to GeneWhiz (South Plainfield, NJ, USA) for Sanger sequencing of the 5′ end of the PCR product using the 515F (V4 region) 16S rRNA gene primer. Sequencing data from GeneWhiz were quality trimmed via Chromatogram Explorer Lite (V 5.0.2), and resulting fasta sequence files were queried using the Ribosomal Database Project's Naive Bayesian rRNA Classifier (V 2.11) with an 80% confidence threshold for taxonomic classification.

Several morphologically identical colonies from each plate streaked for pure culture were inoculated in 3.75 or 7.5 ml (based on colony size) of lysogeny broth (LB) in 15 ml centrifuge tubes and incubated at 37°C under applicable atmospheric conditions. Aerobic cultivars (N=12) were shaken at 60 rpm throughout their incubation while anaerobic cultivars (N=2) were manually inverted two to three times during incubation. Tubes were visually inspected daily for turbidity or pelleting of cells, at which point a glycerol/PBS solution was added to the broth to bring the total volume up to 5 or 10 ml, respectively, at 17.5% glycerol. The glycerol cryoprotectant was added to avoid cell lysis during freezing and transport to Indiana University (Bloomington, IN, USA), where volatile profile analysis was performed. Sterile aerobic and anaerobic samples from two batches of LB with 17.5% glycerol were submitted as baseline controls for volatile profile analyses. Semiquantitative analysis was conducted using the same methods described in experiment 1 for preen oil volatile profiling. We extracted volatile compounds from 2.0 ml aliquots of broth samples using headspace stir bar sampling at 40°C for 60 min in 10 ml ultraclean Gerstel headspace vials following our previously described procedure (Zhang et al., 2005). The glycerol added to the LB resulted in a large band eluting between 12 and 26 min, obscuring smaller volatiles with retention times below 26 min. Therefore, we report only volatiles with higher retention times, which includes all the behaviourally relevant volatile compounds we focused on in the preen oil analysis in experiment 1. To correct for volatiles that may be emitted from the growth medium itself, we averaged non-zero values from the controls for each volatile compound measured and subtracted this average from each cultivar's measured peak areas. We visualized these data through a heatmap created in R using the heatmap.2 function in the gplots package (https://CRAN.R-project.org/package=gplots).

All work was conducted in compliance with the Bloomington Institutional Animal Care and Use Committee guidelines (BIACUC protocols 12-050 and 15-026) and with permission from the US Department of Fish and Wildlife, the Virginia Department of Game and Inland Fisheries, and the US Forest Service.

Using a linear mixed-effects model, we found significant differences in PC1 between pre- and post-treatment samples (F_1,15=6.913, P=0.019), but detected no effect of or interaction with treatment, suggesting that PC1 decreased in both groups over the course of the experiment (Fig. 1). However, PC2 was significantly different for the antibiotic and control groups (F_1,15=6.373, P=0.023) and between the sexes (F_1,15=7.268, P=0.017). Observation of the data suggests that PC2 decreased only in the antibiotic group (Fig. 1), but we did not detect a significant effect of or interaction with time. PC3 was significantly different before and after treatment (F_1,15=−4.90, P=0.043), and appears to have increased in both groups (Fig. 1), with no effect of treatment group or sex.

There were no significant differences in richness (i.e. number of OTUs) or heterogeneity (i.e. Simpson index) between the bacterial profiles of juncos pre- and post-treatment in either the saline control (N=9; richness: paired t-test, t=1.446, P=0.186; heterogeneity: Wilcoxon test, W=42, P=0.021) or antibiotic treatment (N=8; richness: W=21, P=0.743; heterogeneity: W=19, P=0.945) groups after sequential Bonferroni corrections were applied. Furthermore, there were no significant differences in the alpha diversities of the bacterial profiles of juncos receiving saline and antibiotic treatments either before or after the treatments were administered (pre-treatment richness: Mann–Whitney U-test, U=32.5, P=0.771; pre-treatment heterogeneity: U=34, P=0.885; post-treatment richness: t-test, t=−0.375, P=0.713; post-treatment heterogeneity: U=21, P=0.163).

Prior to any treatments being administered, >99% of the bacterial sequences obtained from junco uropygial glands were from three phyla: Proteobacteria (60.5%), Firmicutes (24.3%) and Actinobacteria (14.3%). Ninety percent of sequences from these pre-treatment samples were from 20 OTUs. Seventeen of these OTUs were confidently classified at the genus level. In order of greatest to lowest relative abundance, these OTUs were classified as Staphylococcus, Pseudomonas, Duganella, Rothia (2 OTUs), Comamonas, Janthinobacterium, Stenotrophomonas, Herbaspirillum, Neisseria, Micrococcus, Ralstonia, Sphingomonas, Acinetobacter (two OTUs), Methylobacterium and Paenibacillus.

At the OTU level, the individual identity of juncos accounted for one half of the overall variation in the bacterial profiles of uropygial glands (N=17; R²=0.491). Nevertheless, controlled for individual identity (i.e. strata=ID), there were significant effects of time (two-way NPMANOVA, F_1,30=1.321, P=0.032) and treatment group (F_1,30=1.038, P=0.001) on uropygial gland bacterial profiles. Furthermore, there was a significant interaction between treatment type and time (F_1,30=1.632, P=0.001). The reason for the significant interaction was that both saline [one-way NPMANOVA (strata=ID), F_1,16=1.666, P=0.008] and antibiotic (F_1,14=1.290, P=0.031) treatments affected the bacterial profiles of junco uropygial glands; however, while the bacterial profiles of subjects receiving saline and antibiotics did not significantly differ prior to treatment (one-way NPMANOVA: F_1,15=1.232, P=0.087), they did differ after treatment (F_1,15=1.442, P=0.022; Figs 2B, 3).

SIMPER analyses indicated that seven OTUs were responsible for 75% of the overall variation between the bacterial profiles of saline- and antibiotic-treated birds after treatment (Table 2). Of those seven OTUs, OTU 1 (Pseudomonas) and OTU 2 (Staphylococcus) were less relatively abundant among antibiotic- than saline-treated birds. These OTUs were the two most relatively abundant OTUs pre- and post-treatment, and together they accounted for 36% of the overall variation among sample bacterial profiles post-treatment (Table 2). Most notably, seven out of eight (88%) birds receiving antibiotics exhibited a decreased relative abundance (≥27% of pre-treatment levels) of OTU 2 (Staphylococcus; Wilcoxon test: W=35, P=0.019) within their bacterial profiles. In contrast, the relative abundance of OTU 2 did not differ pre- and post-treatment among saline-treated birds (W=35, P=0.166). The relative abundance of OTU 1 (Pseudomonas) did not consistently change in response to either saline control (W=39, P=0.057) or antibiotic treatment (paired t-test: t=0.486, P=0.642). Nevertheless, OTU 1 was more relatively abundant post-treatment among birds treated with saline than antibiotics (t-test, t=2.145, P=0.049). Only two out of eight birds receiving antibiotics had bacterial profiles in which OTU 1 had a relative abundance >25%, whereas seven out of nine birds receiving the saline control had profiles in which OTU 1 exceeded 25% (Fig. 3).

Given the high degree of variation in uropygial gland bacterial profiles among individual juncos, there was significant, yet not complete, separation of glandular bacterial profiles before and after antibiotic treatment (Fig. 2C). Nevertheless, there was a highly consistent shift in the structure of uropygial gland bacterial profiles in response to antibiotic treatment (Fig. 2C). Specifically, PCoA showed that the glandular bacterial profiles of each antibiotic-treated bird (N=8) shifted to the right along the axis reflecting Coordinate 1, and that this axis explained 33% of the variation among bacterial profiles (Fig. 2C). Among the OTUs constituting >5% of the sequences in the bacterial profiles of these birds, the relative abundance of OTU2 (Staphylococcus; r_S=−0.91, P<0.0001; Fig. 2D) was negatively correlated with values of Coordinate 1, and the relative abundances of OTU 12 (uncl. Sphingomonodaceae; r_S=0.73, P=0.001) and OTU 3 (Herbaspirillum; r_S=0.73, P=0.001) were positively correlated with values of Coordinate 1 (Spearman's correlations with sequential Bonferroni corrections applied).

We detected associations between specific OTUs and each of the three PCs, controlling for treatment (antibiotic and saline control) and time (pre- and post-treatment). For PC1, our information theoretic analysis retained five models from the full set of submodels (Table S2), which we averaged to obtain final model parameters. The averaged model revealed potential influences of OTU 4 (Duganella) relative abundance and time on PC1. The model also retained OTU 5 (unclassified Comamonadaceae) and OTU 3 (Herbaspirillum) relative abundances as variables, as well as sex (Table 3A).

We next tested which bacterial types best predicted variation in PC2. Our analysis identified nine models within 2 AICc units of the top model (Table S3). The final averaged model revealed a significant effect of treatment, and the time–treatment interaction, plus five bacterial OTUs: OTU 1 (Pseudomonas), OTU 2 (Staphylococcus), OTU 10 (Comamonas), OTU 4 (Duganella) and OTU 7 (Janthinobacterium) on PC2 (Table 3B). The model also retained sex, time and the relative abundances of OTU 11 (Stenotrophomonas) and OTU 5 (unclassified Comamonadaceae).

Variation in PC3 was best predicted by sex and the relative abundances of OTU 6 (Rothia) and OTU 11 (Stenotrophomonas) (Table 3). The final model was based on five averaged models (Table S4) and also retained sex and OTU 5 (unclassified Comamonadaceae) and OTU 7 (Janthinobacterium) relative abundance as predictive variables.

A total of 14 distinct cultivars were recovered, 12 under aerobic conditions and two under anaerobic conditions. Aerobic cultivars included Bacillus, Cronobacter, Curtobacterium, Dermacoccus, Kocuria, Lactococcus, Paenibacillus, Pseudomonas, Staphylococcus (two samples), Stenotrophomonas and an unclassified Enterobacteriaceae. The two anaerobic cultivars were identified as Enterobacter and an unclassified Enterobacteriaceae. Notably, Pseudomonas and Staphylococcus, both here cultured from preen oil, were detected, and were frequently abundant, in all 16S rRNA gene sequencing surveys of pre-treatment uropygial gland samples in experiment 1 (N=17). We detected Stenotrophomonas and Enterobacteriaceae in about half of the experiment 1 samples, and detected Kocuria, Lactococcus and Paenibacillus in at least some of those samples.

We successfully characterized the volatile compound profiles of 11 of the 14 cultivars (Fig. 4). Each targeted volatile compound was found in the profiles of an average of 5.25 of the 11 cultivars, and all compounds were associated with at least one cultivar (Fig. 4). Linear alcohols 1-dodecanol, 1-tridecanol, 1-pentadecanol and 1-hexadecanol, and methyl ketones 2-tridecanone, 2-tetradecanone and 2-pentadecanone were the most common volatiles identified, with each being found in the profiles of at least half of the cultivars. Paenibacillus, Kocuria, Curtobacterium and Pseudomonas had the most compound-rich volatile profiles (Fig. 4).

We present two exploratory studies evaluating whether symbiotic bacteria associated with uropygial glands produce volatile compounds used as songbird chemical cues. Our results provide evidence that the composition of resident uropygial gland bacterial communities affects the volatile profiles of preen oil, and that these bacteria can produce behaviourally relevant volatile compounds independent of their host. Specifically, by administering antibiotics to dark-eyed junco uropygial glands and culturing and characterizing bacteria from these glands, we found the following: (1) Pseudomonas and Staphylococcus were less abundant in the uropygial glands of birds treated with a broad-spectrum antibiotic compared with control birds that received saline; (2) when comparing preen oil volatile profiles, PC2 – the principal component score most strongly influenced by the abundance of 2-tridecanone, 2-tetradecanone and 2-pentadecanone – was the only variable significantly affected by antibiotic treatment; (3) PC2 was also significantly related to the relative abundance of Pseudomonas and Staphylococcus in the uropygial gland; specifically, reduced relative abundances of these two bacteria were strongly associated with a reduction in PC2; and (4) bacteria cultivated from junco preen oil, including Pseudomonas and Staphylococcus, independently produced junco preen oil volatile compounds, including 2-tridecanone, 2-tetradecanone and 2-pentadecanone.

Staphylococcus and Pseudomonas are commonly detected residents in junco uropygial gland and feather bacterial communities. In previous studies, we detected these two bacterial genera in the glands of over 75% of the juncos surveyed (Whittaker et al., 2016; Whittaker and Theis, 2016). They have also been cultivated from junco feathers (Dille et al., 2016) and in the present study, we cultivated them from junco preen oil. Avian isolates of these two genera have been shown to be highly sensitive to enrofloxacin (Eid et al., 2019), further corroborating our conclusion that the change in volatile profiles was likely due to antibiotic treatment reducing the relative abundance of these genera in the uropygial glands. Although enrofloxacin is broadly effective against both Gram-negative and Gram-positive bacteria, in the present study it did appear to have a greater effect on Staphylococcus (a Gram-positive bacterium) than on Pseudomonas (Gram-negative), which may be related to increased antibiotic resistance in Gram-negative bacteria including Pseudomonas (Hancock, 1998).

Many of the bacteria identified through molecular surveys in experiment 1 and through culture in experiment 2 have previously been identified as members of plumage or uropygial gland microbiota in other bird species. Staphylococcus was one of the most common bacteria found on blue petrel (Halobaena caerulea) plumage (Leclaire et al., 2019), and it has also been found on chicken (Gallus gallus), least auklet (Aethia pusilla) (Shawkey et al., 2006) and house finch (Carpodacus mexicanus) (Shawkey et al., 2003) feathers. Bacillus, Pseudomonas and Stenotrophomonas have been found on crested auklet (Aethia cristatella) (Shawkey et al., 2006) and eastern bluebird (Sialia sialis) plumage (Shawkey et al., 2005). Kocuria, identified here in experiments 1 and 2, is also associated with birds, and recently described species were isolated from the uropygial glands of the great spotted woodpecker (Dendrocopos major) (Braun et al., 2018) and American barn owl (Tyto furcata) (Braun et al., 2019a; Braun et al., 2019b).

Each of the bacterial cultivars from preen oil in experiment 2 produced one or more of the 16 junco volatile compounds measured independent of their junco hosts. Notably, these included compounds that have been shown to be behaviourally relevant for juncos: specifically, the methyl ketones 2-tridecanone, 2-tetradecanone and 2-pentadecanone (Whittaker et al., 2013, 2011a), which were also the three compounds that most strongly influenced PC2 in experiment 1. In a previous study, these three compounds predicted individual reproductive success among juncos, in opposite directions for the sexes: males with higher proportions of these compounds sired more offspring, while females with lower proportions produced more offspring (Whittaker et al., 2013). Here, in experiment 2, 2-tridecanone, 2-tetradecanone and 2-pentadecanone were produced by 10, seven and six, respectively, of the 11 bacterial cultivars whose volatile profiles were characterized. The Pseudomonas isolate produced all three of these compounds while Staphylococcus produced 2-tridecanone.

The results for volatile profile PC1 (influenced by linear alcohols 1-decanol through 1-octadecanol) and PC3 (dodecanoic acid, tetradecanoic acid and hexadecanoic acid) are less clear. Most of the linear alcohols and carboxylic acids that contribute to PC1 and PC3 were produced by several (mean 4.9, range 1–10) bacteria cultured in experiment 2 (Fig. 4), and most of these cultured bacteria were also detected in the molecular surveys of experiment 1 (Pseudomonas and Staphylococcus were abundant across all individuals; Stenotrophomonas, Enterobacteriaceae, Kocuria, Lactococcus and Paenibacillus were present in some individuals). However, the bacteria most closely associated with PC1 in experiment 1 were Duganella and Herbaspirillum, both of which are recognized as common contaminants of DNA extraction kits and/or PCR and sequencing reagents (Glassing et al., 2016; Salter et al., 2014). Similarly, PC3 was most closely associated with Rothia, another common contaminant. Although these three OTUs were not statistically identified as contaminants in our samples (i.e. using the program decontam; Davis et al., 2018), we urge caution in interpreting these results. PC3 was also significantly related to the relative abundance of Stenotrophomonas, which we were able to culture from preen oil and which was detected among experiment 1 sequences. However, we were unable to analyse the volatile profile of the Stenotrophomonas cultivar and thus could not directly demonstrate the production of these compounds by this bacterium in preen oil.

The small sample size in experiment 1, coupled with the high degree of inter-individual variation in uropygial gland microbial profiles in general, resulted in high levels of variation in phenotypic responses to the experimental treatment with antibiotics. Pre-treatment differences between the volatile profiles of treatment and control groups may have impeded our ability to detect more effects of antibiotic treatment on odour. Furthermore, we found that bacteria strongly influence the abundance of certain volatile compounds that have been previously shown to be associated with male and female reproductive success in opposite directions between the sexes; unfortunately, our sample size was not sufficient to detect sex differences in either volatile abundance or response to treatment among juncos in this study. We also cannot rule out the possibility that antibiotic treatment had effects on the birds' physiology, which may have indirectly affected their volatile profiles. However, we injected the antibiotic solution directly into the preen gland rather than administering it through the diet, thereby reducing the likelihood of systemic effects. Finally, we did not quantify absolute bacterial loads in uropygial glands before and after the experiment, and thus we cannot confirm that antibiotic treatment reduced bacterial presence in the glands.

The finding that bacteria associated with avian preen glands are capable of producing primary components of avian chemical cues and signals raises questions about the relationships between birds and their bacteria, and how these bacterially generated cues could honestly reflect host traits (Archie and Theis, 2011; Carthey et al., 2018; Maraci et al., 2018). Group members and social partners in general share microbes through social interactions and physical contact (Archie and Tung, 2015; Ross et al., 2017). In a study of junco pairmates and nestlings, the uropygial gland microbiota of nestlings more strongly resembled those of their mothers, which brooded and fed them, than those of their attendant fathers, which only fed them; furthermore, the paternity status of attendant males did not affect their microbial similarity to the nestlings (Whittaker et al., 2016). Similarly, in a cross-fostering experiment, nest environment, not genetic relatedness, shaped microbial profiles in hoopoe nestling uropygial glands (Ruiz-Rodríguez et al., 2014). Because social behaviour can influence the structure of symbiotic microbial communities, which can, in turn, influence the structure of host chemical cues and signals, more work is needed to understand how symbionts may be able to reliably signal static information (such as underlying genotype) about their host. Given that host animals and their microbial symbionts are often co-diverged, host genotype itself may play an important role (Brooks et al., 2016; Goodrich et al., 2014; Hooper et al., 2012). Specifically, the host's MHC genotype may mediate the effects of social behaviour and microbial transfer between individuals and thereby promote honest information of the underlying genotype. MHC genotype is known to influence mate choice in many vertebrates, including songbirds, and it can be communicated via chemical signals (Penn and Potts, 1999; Setchell et al., 2011). These relationships among host association patterns, genotypes, microbiomes, chemical signals and mating behaviours are fruitful areas for future research.

We thank the University of Virginia, Mountain Lake Biological Station, and the Mountain Lake Lodge. We thank Rachel Hanauer and Nicole Gerlach for assistance with fieldwork.