Sampling locality is more detectable than taxonomy or ecology in the gut microbiota of the brood-parasitic Brown-headed Cowbird (Molothrus ater)

Hypotheses

The four hypotheses we evaluated concerning the structuring of cowbird microbiota are not mutually exclusive. The Nature Hypothesis posits that cowbird gut microbiota are determined predominantly by their biological cowbird parents (Fig. 1A), with the brood host contributing little to the microbiota of the cowbird young. The Nature Hypothesis would be corroborated if we found that the gut microbiota of cowbird samples are most similar to those of other cowbirds, rather than to the microbiota of brood host species. Under the alternative Nurture Hypothesis, the brood host exerts the most influence on its brood’s gut microbiota, regardless of its genetic relationship to them (Fig. 1B). This hypothesis predicts that the microbiota of each cowbird should be more similar to that of its brood host species than to other cowbirds. The Nurture Hypothesis would be supported if clusters of cowbird and brood host individuals having similar gut microbiota are detected.

The Environment Hypothesis (Fig. 1C) posits that the local environment (e.g., climate, flora, fauna, etc.) accounts for similarity of gut microbiota. It predicts that the microbiota of birds in closer geographic proximity will be most similar, despite their genetic background, ecology or evolutionary history. Finally, cowbirds may have different gut microbiota assemblages at different life stages (the Convergence Hypothesis, Fig. 1D). One prediction is that younger cowbirds have a more diverse microbiota assemblage that is able to utilize a variety of diets but which converges to a stable cowbird-like microbiota as the cowbirds mature. This will be supported if the microbiota of adult cowbirds form their own clade (to the exclusion of the juveniles) or if the gut microbiota of juveniles are taxonomically more diverse than adults. (“Juvenile” is defined here as an individual having the combination of juvenile plumage, a completely unossified cranium, presence of a Bursa of Fabricius and undeveloped gonads.)

Sampling

We sampled from birds available in the LSU Museum of Natural Science’s Collection of Genetic Resources for this investigation, including 32 cowbirds (Icteridae: Molothrus ater) from Louisiana (LA, near Baton Rouge) and California (CA, near Weldon) and 16 individuals from nine known parasitized passerine species, all but one from LA (Table 1, more detailed specimen information in Table S1): Northern Cardinal (Cardinalidae: Cardinalis cardinalis), House Finch (Fringillidae: Haemorhous mexicanus), Orchard Oriole (Icteridae: Icterus spurius), Indigo Bunting (Cardinalidae: Passerina cyanea), Blue-gray Gnatcatcher (Polioptilidae: Polioptila caerulea), Prothonotary Warbler (Parulidae: Protonotaria citrea), Carolina Wren (Troglodytidae: Thryothorus ludovicianus), White-eyed Vireo (Vireonidae: Vireo griseus), Hooded Warbler (Parulidae: Setophaga citrina). Birds were collected under a general collection permit for another project (USFWS scientific collecting permit MB679782; LDFW permit LNHP-13-032) and frozen within two hours of their collection. Throughout the manuscript, individuals are identified by their common names and an individual identifying number. One individual (NorthernCardinal4) served as a replicate to assess PCR/sequencing bias; these samples are identified as NorthernCardinal4.1 and NorthernCardinal4.2 and bring the total number of samples to 49.

DNA extraction, amplification, sequencing and quality control

The entire digestive tract was removed when the bird was thawed for museum specimen preparation. Total DNA was immediately extracted from the contents of the large intestine, halfway between the ceca and cloaca, using a QIAamp DNA Stool Mini Kit (QIAGEN). Following Gloor et al. (2010), we used combinatoric primers and massive multiplexing of PCR amplicons for sequencing on one lane of an Illumina Hi-Seq. This method uses paired-end sequencing technology to generate pairs of sequences with 100% overlap across variable region 6 (V6) of the 16S component of rRNA; primer sequences align to positions 967–985 and 1078–1061 on Escherichia coli 16S rRNA (Gloor et al., 2010).

We used several measures of sequence quality control. First, both reads of a given pair had to match across 100% of the bases. The pairs also must have no errors in tag sequence or priming sequence. We used the Bellerophon (Huber, Faulkner & Hugenholtz, 2004) function within the mothur program (Schloss et al., 2009) to identify and discard potentially chimeric sequences. Finally, we used mothur to discard sequences that did not blast to the domain Bacteria. The reads passing these filters were included in the final dataset.

Clustering analyses

Individuals were partitioned into four datasets in order to test the hypotheses outlined above: All Birds (N = 49), All Louisiana Birds (N = 35), Cowbirds Only (N = 32), Hosts Only (N = 17).

The microbial ecology package QIIME (Caporaso et al., 2010) was used for the following analyses. First, the reads were assigned to phylotypes at 97% sequence similarity because 3% is frequently cited as the “species” level of microbial taxonomy (Schloss & Handelsman, 2005). Next, we assigned taxonomies to OTUs using the RDP Classifier Program (Wang et al., 2007), with the default confidence threshold of 80%. A pairwise matrix of distances between each gut microbial community (i.e., each bird specimen) was constructed using UniFrac (Lozupone & Knight, 2005). UniFrac distances are calculated based on the amount of branch length in a phylogenetic tree that is unique to either of two environments (versus how much of the tree is shared by the environments). These distances can be weighted by abundance or be based on presence–absence of OTUs (unweighted); we use both distance matrices for analyses because it is unclear which method is best for representing gut microbiota samples and is more transparent than choosing one arbitrarily. Our microbial phylogenetic tree was constructed with FastTree (Price, Dehal & Arkin, 2009). To reduce the effects of sampling (sequencing) bias, all individuals were randomly reduced to 5 018 reads, equal to the lowest number of reads for any bird in the dataset.

We constructed UPGMA dendrograms based on both the unweighted and weighted UniFrac distances to visually represent the relatedness of the gut microbiota for all five datasets and test the hypotheses. As a confidence metric, we jackknifed the trees using the QIIME recommendation of 75% of the reads used in the rarefaction (3 760) with 10 replicates. Principal coordinates analysis (PCoA) was also performed on both the weighted and unweighted UniFrac distance matrices.

As a complement to the phylogenetic-based methods, we visualized the data with nonmetric multidimensional scaling (NMDS). We square root-transformed the percentage of each sample that belonged to each bacterial phylum, then created a pairwise distance matrix using Bray–Curtis dissimilarity, applied through the vegdist function of the vegan package (Oksanen et al., 2011) in R (R Development Core Team, 2010). The nmds function of the ecodist package (Goslee & Urban, 2007) was then used to calculate the two-dimensional positions of the samples (such that closer samples are more similar), the stress and R² value of the plot. Stress values >0.3 should not be considered valid whereas values <0.2 can be considered a good representation of the data (Quinn & Keough, 2002).

To specifically test the Convergence Hypothesis, we compared UniFrac distances between and within adult cowbirds, juvenile cowbirds and brood hosts. If adult cowbirds converge on a cowbird-specific microbiota, adults will have lower pairwise distances than within juveniles or either category to brood hosts. Both weighted and unweighted UniFrac distances were assessed.

Categorical variable significance

To further explore the factors correlated with gut microbiota, we tested for statistical associations between categorical metadata associated with each bird and the UniFrac distance matrices (Table S1); we used the statistical tools Adonis (McArdle & Anderson, 2001) and Anosim (Clarke, 1993) implemented in QIIME. The categorical variables included family, genus, species, age (either “juvenile” or “adult” based on percent of skull ossification), locality (LA or CA), diet (the known dietary specialization of the species: mostly plant material, mostly animal material, both animal and plant material) and stomach contents (what was actually found in the stomach, e.g., “insects” or “white millet”). We also tested the total number of bacterial phyla detected per bird (richness) to see if the phylum-level diversity of the established microbial community was associated with the communities. We calculated significance of all variables for both the weighted and unweighted UniFrac distance matrices with 999 iterations; we also repeated analyses with a higher number of sequences (and fewer individuals) for each dataset to see if the signal changed when more data were included. Datasets were rarefied to 17 000 and 42 000 sequences, reducing the number of birds to 46 and 35, respectively.

Based on results, we ran additional Adonis tests to partition the variation in the samples as due to the age and locality variables. We used the All Birds dataset and analyzed the weighted UniFrac distance matrix, unweighted UniFrac distance matrix and a sites (birds) by species (bacterial phyla) matrix, where cells were assigned the value of the number of sequences belonging to each phylum for each bird. We used the adonis function of the vegan package in R and performed 999 iterations, constraining resampling to be within species. Since the order of variables being tested matters, we tried age then locality as well as locality then age for each of the datasets.

Comparison to mammals and insects

To put avian gut microbiota in a broader context, we compared our results to a mammal dataset (Ley et al., 2008) containing 56 individuals from 56 species across 13 orders (Table S1) and an insect meta-analysis dataset containing 85 individuals from 62 species across seven orders (Colman, Toolson & Takacs-Vesbach, 2012). Although the mammal and insect datasets were collected with different methods than those outlined above, most sequence fragments contained the V6 region. We pruned all reads to the same homologous region and length for analysis. We only analyzed samples with greater than 200 sequences. To increase coverage for some species we combined mammals belonging to the same species into single samples. This treatment should not skew the results of our analysis, because Ley et al. (2008) found that individuals from the same species clustered together. We taxonomically assigned reads using RDP Classifier Program. For PCoA, we rarefied all samples to 200 reads and used the unweighted UniFrac distances as input. We also performed NMDS on samples, as described above. To test for significant associations between class, order and diet categories (herbivore, carnivore, omnivore), we tested each variable against both the weighted and unweighted UniFrac distance matrices in the same manner as above.

Clustering analyses

Cowbird samples did not cluster together in the UPGMA dendrogram of All Birds (Fig. 3). Brood host species having more than one individual also did not cluster together, even when cowbird samples were excluded from the analyses (Fig. S1). NMDS representation of All Birds showed little segregation by age or locality (Fig. 3), although the stress of the plot was low (0.1395). In general, UPGMA dendrograms and NMDS plots of all datasets showed little clustering by bird species and high levels of variation (Fig. S1). The two replicate samples, NorthernCardinal4.1 and NorthernCardinal4.2 were most closely related to each other in every analysis.

Pairwise distances were assessed between and within adult cowbirds, juvenile cowbirds and brood hosts to test whether adult cowbirds converged on a cowbird specific microbiota (i.e., if adult cowbird microbiota were more similar to each other than they were to other group comparisons or other groups were to themselves). All pairwise comparisons between and within groups had largely overlapping distributions (Fig. S2) for both weighted and unweighted UniFrac distances, and thus, adults were not more similar to each other than other comparisons.

Categorical variable significance

There were a total of 12 statistical tests computed for each categorical variable for each dataset (Fig. 4E). Figure 4 shows how frequently significant each of the variables was in each dataset; generally speaking, locality was most frequently significant, followed by age and then diet. Taxonomic categories, stomach contents and bacterial richness were generally not significant.

The multifactorial Adonis tests to assess locality and age were run twice, varying the order of the variables, since this can have an affect on the results on three distance matrices and we had no a priori reason to prefer one variable as being more important than the other. Age was significant in all six tests and locality was significant in one of six tests (Table 2).

Comparison to mammals and insects

For all analyses, we only included individuals with more than 200 sequences; for the PCoA, we randomly pruned all individuals to 200 sequences. This reduced the mammal dataset to 38 samples belonging to 11 orders and the insect dataset to 22 individuals belonging to four orders (Table S1). Consistent with Ley et al. (2008), mammals were dominated by Firmicutes and Bacteroidetes (Fig. 2C), whereas bird samples were predominately Firmicutes, with some samples having mostly Proteobacteria or Bacteroidetes (Fig. 2A). Insects also generally contained a majority of Firmicutes although individual samples varied between 100% Proteobacteria and 98% Firmicutes (Fig. 2B). The PCoA showed birds as distinct from mammals and insects, which largely clustered into their respective groups but contained some overlap (Fig. 5A); birds spanned a greater portion of PC2 than all Mammals. NMDS was broadly overlapping but birds clustered together in the middle of the plot (Fig. 5B). Independent Adonis tests for significance of class, order and diet categories revealed highly significant associations (p < 0.01) between the gut microbiota and all three variables (Table 3), except for diet and the weighted UniFrac distance matrix, which was not significant (p = 0.206).

Hypotheses

Across analyses and datasets, cowbirds and hosts were interspersed (Fig. 3)—neither cowbirds nor individuals from each brood host species clustered together. We therefore reject the Nature Hypothesis because there is no evidence for a specific cowbird gut community shared by all cowbirds. We are also unable to appropriately evaluate the Nurture Hypothesis, since an underlying assumption was that brood host species would be most similar to one another. The Environment Hypothesis is most strongly supported, since the locality variable is most frequently correlated with cowbird gut microbiota. Finally, we reject the Convergence Hypothesis, because adult cowbirds were no more similar to each other than they were to any other group.

Several aspects of these hypotheses warrant further discussion. The importance of locality and diet are consistent with an important role of ecology in shaping gut microbiota, so the Nurture Hypothesis remains a possible explanation of cowbird gut microbiota. For example, if the factors that shape gut microbiota are drawn from largely overlapping ecological niche space, a lack of bird taxonomic signal and generally low levels of clustering would be expected, in addition to the importance of parameters like locality and diet. An ideal experiment to test the Nurture Hypothesis would be to sample entire brood families from a single nest (parents, offspring, brood parasites) plus the cowbird mother.

Another caveat with the Nurture Hypothesis comes from the sampling, which included only samples that LSUMNS had previously collected. While some juveniles are most similar to brood host individuals, we have no way to evaluate whether that is because they were raised by that species. The Brown-headed Cowbird parasitizes approximately 250 bird species, only nine of which are represented here and of those sampled, there is an uneven distribution across age and locality (Table 1). This sampling is not ideal, however one priority of this study was to explore the amount of information gained through the collection and analysis of gut microbiota of museum specimens collected for other research purposes; therefore, we did not design our own sampling scheme so much as utilize previously sacrificed animals. In other words, by taking advantage of recent technological (sequencing) advances, can we increase the amount of data gained from a given sampling effort? In the future, field biologists and microbial ecologists may work together to derive increasingly more data from museum collecting efforts and address interdisciplinary questions. Even modest microbiota analyses can provide the foundation for and preliminary answers to novel questions. For example, based on our results, targeted future studies should include as much taxonomic breadth as possible, with replicates from each locality. It was not anticipated that brood hosts would not cluster together, because they are genetically, ecologically and presumably environmentally most similar. We recommend including more than one individual from each species, when possible. Finally, regarding sampling, it is important to note that the single adult we have from California is a House Finch, which is parasitized by cowbirds but does not successfully rear their young (Kozlovic, Knapton & Barlow, 1996; Lowther, 2013). The inclusion of this cowbird “victim” has little affect on our interpretation of the data and could even be expected to be an outgroup to the other samples which are ecologically similar enough to raise cowbird chicks. However, House Finch does not appear significantly different from the other parasitized species. Its inclusion in the study has no effect on the major findings, that (1) cowbirds do not form their own cluster and (2) locality is the strongest correlate to cowbird gut microbiota.

Although we reject the Convergence Hypothesis, we leave room for a transition occurring sometime between hatching and fully-grown adult cowbirds. Furthermore, the scale of our analyses may have missed the critical transition period between young and adult microbiotas—which may vary by both bird and microbe species (Moreno et al., 2003; Scupham, 2007; Scupham, 2009)

High inter-individual variation appears to be a hallmark of microbiota studies (e.g., Turnbaugh et al., 2010; Dethlefsen et al., 2006; Ruiz-Rodriguez et al., 2009b), so much so that the concept of a “core microbiota” is in doubt (Lozupone et al., 2012). The birds in this study belong to a single order, the Passeriformes, and appear to have more variation than any of the non-avian orders we analyzed, with the exception of Hymenoptera (Fig. S3). Cowbirds in particular appear to have a highly variable gut microbiota, which may be a species level trait. The relative contribution of Proteobacteria and Firmicutes to cowbirds spans nearly all brood-host species and the cowbirds span nearly the entire dendrogram and NMDS plot (Fig. 3). However, the four Northern Cardinal samples appear to contain as much variation as the cowbirds (Figs. 2 and 3), so high variation in gut microbiota may be an order level trait.

Despite high levels of variation, taxonomic signals of the vertebrate hosts are frequently detected in microbiota studies (e.g., Ochman et al., 2010; Dewar et al., 2013; Ley et al., 2008), as are dietary classifications (Ley et al., 2008; Anderson et al., 2012). Our comparison of birds, mammals and insects revealed the taxonomic classes to be visibly and significantly different (Table 3, Fig. S3). Even at the order level, clustering was observed in all classes. No clustering by dietary classification was observed across or even within the classes (Fig. S3).

This study was conducted using a single marker and it relies on OTUs delimited from these genetic data. Metagenomic studies randomly sequence many loci across a microbiota sample and show that while taxonomic identity can vary widely across individuals, functional groups are highly conserved (Lozupone et al., 2012). The overlapping ecology, lack of taxonomic signal and significant effect of sampling locality indicate this would be an interesting application of metagenomics. Does the brood parasite’s gut contain more functional categories than a traditional bird? What, if any, functional categories are most represented? Retesting all the above hypotheses with metagenomic data instead of fingerprint data would undoubtedly be valuable.