Comparative Analyses of Vertebrate Gut Microbiomes Reveal Convergence between Birds and Bats

1The mammalian gut microbiome has emerged as a key regulator of host physiology, and coevolution between host and microbial lineages has played a key role in the adaptation of mammals to their diverse lifestyles. Diet, especially herbivory, is an important correlate of microbial diversity in mammals. Most mammalian microbiomes are also strongly correlated with host phylogeny, despite profound shifts in diet. This suggests host factors that themselves change across host phylogeny, such as gut physiology, play an important role in structuring the gut microbiomes across mammals. The vertebrate adaptive immune system is even speculated to have evolved as just such a factor for selective maintenance of symbiotic homeostasis.

2The importance of phylogeny-correlated factors to the diversity of vertebrate microbiomes more generally is still poorly understood. Phylosymbiosis, or the observation that more closely related host species have more similar microbiomes, has been described in a number of nonmammalian taxa. Other analyses have found substantial variation in phylosymbiotic signals among mammalian taxa, sometimes with conflicting results. The presence of a robust phylosymbiotic correlation implies that host factors control microbial assembly. Even if the specific mechanisms are unknown, variation in the strength or presence of a measurable phylosymbiotic signal across host phylogeny could prove useful for identifying such mechanisms through comparative studies. However, most studies to date have focused on just a few taxa at a time, and variable methods for both surveying the microbiome and measuring phylosymbiosis and host specificity (or the restriction of microbes to specific host lineages) have made generalizations difficult.

3Without broader evolutionary context, it is unclear how universally conserved patterns of host-microbe phylosymbiosis actually are. Growing evidence indicates that the strong patterns identified in mammals are the exception rather than the rule in vertebrates. Meta-analyses of fish and birds have failed to detect the strength of correlations to diet and phylogeny reported in mammals. A recent analysis of samples from more than 100 vertebrate species also found the strength of phylogenetic correlation to be much higher in mammals than in birds, reptiles, amphibians, or fish. It is increasingly appreciated in nonvertebrate animals that fundamental aspects of the host’s relationship to its symbiotic community can change drastically between taxa: many insects depend entirely on microbes for key metabolites, while others seem to be devoid of resident gut microbes. The complexity of different factors likely influencing diverse vertebrate gut communities remains a challenge to uncovering the most important causal relationships. The diversity of animal physiologies, habitats, and lifestyles offers opportunities to use convergent “natural experiments” in evolutionary history to sort through this complexity. By filling out our knowledge of microbial diversity across the vertebrate tree of life, we may be able to use major transitions in the host-microbiome relationship to identify elements of animal biology that are most likely to play such causal roles.

4Reasoning that convergences in host phenotypes offer some of the clearest such natural experiments, we have assembled a large sample set using the Earth Microbiome Project standard methods to identify patterns of phylosymbiosis and host specificity across vertebrates. Subsets of these data have previously explored convergences in diet, including convergent evolution of myrmecophagy in mammals, folivory in primates, and blood feeding in birds and mammals, as well as the impacts of captivity on the gut microbiome. Here, we present our first analysis of the complete data set, spanning samples from 315 mammalian species, 491 avian species, and 86 species representing other vertebrate classes (Fig. 1; see Data Set S1 in the supplemental material). Through large-scale vertebrate host sampling using consistent methods, we are able to compare broad patterns in the diversification of microbiomes both within and between major vertebrate lineages with substantially greater power.

5Principal-coordinate analysis (PCoA) of the complete microbiome data set reveals a striking pattern: although microbiomes from each vertebrate class generally cluster in PCoA space (permutational multivariate analysis of variance [PERMANOVA] R2 = 0.04897, P = 0.001), bats group more closely with birds than with other mammals or any other vertebrate class (Fig. 2). Furthermore, both of these clades harbor relatively low proportions of Bacteroidetes but high proportions of Proteobacteria (see Fig. S1C and D in the supplemental material), a phylum of bacteria highly associated with birds and flight (Fig. S1A and B) that is diminished in most nonbat mammalian hosts. Yet, despite these compositional similarities in the gut microbiomes of birds and bats, we find that the strongest similarity between these clades is a lack of strong association with a specific microbiome (Fig. 1).

6In birds, we find gut microbiota in general to be only weakly correlated with host phylogeny (multiple regression on distance matrices [MRM] R2 = 0.02, P = 0.001) (Fig. S2B; Table S1) and not correlated with host diet (P = 0.415) (Fig. S3A) despite diet varying widely among bird lineages; this is juxtaposed with mammals, in which we see a strong correlation to both diet and phylogeny (R2 = 0.17, P = 0.001) (Fig. S2A to S3B and Table S1). Rather, correlations between bird gut microbiota and host phylogeny are more similar to levels in the less extensively sampled nonavian reptiles (R2 = 0.03, P = 0.002) and amphibians (R2 = 0.03, P = 0.021). An analysis of beta diversity through time (BDTT), a technique previously shown to deconvolute the contributions of host diet and phylogeny to microbiome diversity in mammals, also showed little correlation between bird diet and phylogeny at any depth of the microbial tree (Fig. S4C). Similarly, bats had lower correlations between gut microbiota and both diet and phylogeny than other well-sampled mammalian lineages in a BDTT analysis. Correlations between gut alpha diversity and host body mass also held much more weakly for all birds (P = 0.0101) than mammals (P < 0.001) (Fig. S5C), but only held for flightless mammals when considering bats and flightless birds separately (Pbirds, flighted = 0.08, Pbirds, flightless = 0.69, Pmammals, flighted = 0.07, Pmammals, flightless < 0.001) (Fig. S5D).

7The lack of correlation between host factors and microbial composition in birds appears to be partly a consequence of an overall lower level of specificity between microbial and host taxa than is found in mammals. While in mammals, most amplicon sequence variants (ASVs) are only found in a single host order, the opposite is true in birds: most ASVs are shared broadly across host taxa (Fig. 3). Using a standardized effect size measurement of Pielou’s evenness statistic to control for imbalances in our sampling, we find that microbial communities associated with mammals are generally much less evenly distributed across host taxa than are the microbial communities of birds (Fig. 4, inset). In particular, bats and other mammals from insectivorous and carnivorous orders harbored the least specific communities; while among the birds, the flightless orders (ostriches, emus, cassowaries, kiwis, and rheas), and the weak fliers from the related Tinamiformes (tinamous) were among those hosting communities with median specificities approaching those of mammals. The microbiota of these bird orders (constituting the Palaeognathae) were also more likely to occur in mammals, with samples from the Struthioniformes (ostriches), Rheiformes (rheas), and Tinamiformes hosting the least bird-specific microbiota, and samples from the Eulipotyphla (insectivorous shrews and moles), Dasyuromorphia (carnivorous marsupials), Pholidota (ant-eating pangolins), and Chiroptera (bats) hosting the least mammal-specific microbiota (Fig. 4).

8The association between flight and the loss of microbiome specificity is also reflected in measurements of phylosymbiosis, suggesting that the evolution of powered flight profoundly disrupts the pattern of phylosymbiosis between host and gut microbiome. We compared measurements of microbiome dissimilarity to host phylogenetic distance (measured as branch length on the TimeTree) and approximate time-calibrated host phylogeny using the Mantel test (Fig. 1). We see that in mammals as a whole, the correlation between microbial community similarity and host phylogenetic distance is high, especially at recent timescales (Fig. S4A and B), whereas this correlation is much lower in birds (mammals: P = 0.001, r = 0.40, coefficient = 2.61E−4; birds: P = 0.001, r = 0.14, coefficient = 8.92E−5) (Fig. 5, left). In contrast, among mammalian orders, the strength of this correlation varies dramatically, with the steepest relationships between phylogenetic distance and microbiome turnover in the fermentative herbivores; among the well-sampled mammalian orders (≥7 species represented), bats were unique in having little correlation at all (Fig. 5, right).

9Our results suggest that the evolution of powered flight affected interactions between hosts and microbes in birds and bats convergently: birds and mammals with powered flight both display markedly less specificity than most nonflying mammals and somewhat less specificity than nonflying birds. We considered several potential proximate explanations for the convergent loss of microbial specificity in birds and bats, which can be broken down into ecological, evolutionary, and technical factors.

10First, flight could expose individuals to more diverse-source microbial communities, resulting in increased microbial gut diversity if assembly of these communities was a neutral process. This seems unlikely to be an important factor in itself, as vertebrate gut communities have mostly been shown to follow nonneutral dynamics. If this were the case, however, we might expect both birds and bats to host especially high alpha diversity, reflecting colonization from a broader range of environments. Yet, the opposite appears to be the case (see Fig. S5A in the supplemental material): both birds and bats have lower mean alpha diversities than nonbat mammals (pairwise t tests, P < 0.001 for observed operational taxonomic units [OTUs], Shannon diversity, and Faith’s phylogenetic diversity [PD]). Furthermore, we did not see greater alpha diversity in migratory birds, which presumably encounter a greater range of environments than nonmigrators (P = 0.225 for observed OTUs) (Fig. S5B).

11Similarities in microbiome composition and patterns of association suggest the presence of flight-correlated host factors exerting similar selective pressures on the assembly of the gut microbiome. Birds and bats both tend to have reduced intestinal lengths and shorter intestinal content retention times, perhaps as a by-product of selective pressure to decrease mass for more efficient powered flight. Shorter guts, with correspondingly reduced anaerobic volumes, could also present less of a barrier to microbial exchange through an aerobic environment. Accordingly, relative to that in nonflying mammals, both birds and bats have fewer obligate anaerobes (P < 0.001) and more facultative anaerobes (P < 0.001), and flightless birds have more obligate anaerobes than flighted birds (P = 0.001) as a proportion of their gut microbial communities, as measured by predicted high-level phenotypes using BugBase (Fig. S6).

12It is possible that some other aspect of adaptation to flight has had a similar net effect on microbiome specificity. The metabolic demands of powered flight are profound and appear to have led to extensive and surprising degrees of convergence. Birds have the smallest genomes of all amniotes, and birds and bats both have more compact genomes than their flightless relatives, most likely due to extensive DNA loss—hypothesized to be related to powerful ongoing selection to reduce mass and enhance flight-related correlates of metabolism. It is therefore plausible that pressure to decrease mass may extend to microbial biomass. Both birds and bats also have greater rates of intestinal paracellular absorption than nonflying vertebrates, meaning that a higher proportion of simple nutrients are absorbed directly by the host, potentially decreasing the role for symbiotic microbial metabolism. Indeed, some studies of microbial biomass across animal hosts have reported that both birds and bats carry much lower numbers of microbial cells in feces than nonflying mammals. A selection toward microbial reduction may also explain in part why birds and bats both seem to have lost an association with Bacteroidetes but retained an association with Proteobacteria, a group of bacteria suggested to have high functional variability: to maximize microbial function while also reducing diversity and mass. Even if selective pressure to reduce overall microbial biomass did not, in itself, lead to the loss of tight evolutionary ties with specific symbionts, consistently lower resident microbial biomass might increase the proportion of transient environmental microbes relative to the total pool of microbes in the sample, decreasing the observed degree of host specificity on average. Indeed, this could also help explain the dominance of Proteobacteria, which make up a large proportion of the airborne microbiome, a source to which flighted animals are constantly exposed.

13If the convergent loss of microbial specificity in birds and bats is ultimately related to selection toward a decreased reliance on microbes, it raises the possibility that the proximate drivers include parallel changes in mechanisms mediating host-microbe interactions. Consistent with this possibility, both birds and bats have lost large numbers of genes involved in immune function, concomitant with a reduction in genome size, including reduced interferon (IFN) locus copy number. In bats, genes thought to be involved in minimization of DNA damage as a metabolic consequence of flight have been linked to immune function and are also under positive selection, suggesting that adaptation to flight may have affected key aspects of immunity. Continuous and dampened expression of certain immune genes has been proposed to allow bats to simultaneously manage flight-induced immune damage and to tolerate a large number of viral pathogens without experiencing disease. Interestingly, our results suggest that bats also host significantly more potentially pathogenic bacteria than other mammals (Tukey’s post hoc test, P < 0.001), implying that links between flight and antiviral immunity may affect antibacterial immunity as well.

14We cannot completely discount the possibility that the similarities we observed between bird and bat gut microbiomes are due at least in part to factors also correlated with flight. For example, birds and bats both excrete uric acid, in contrast to most other mammals. If uric acid or associated compounds interfere with DNA extraction or somehow compromised the microbial DNA prior to extraction, our observed sequence profiles might not accurately represent the relevant communities in vivo. However, nonavian reptiles, which also excrete uric acid, resemble nonflying mammals in their microbiome taxonomic compositions more than they resemble either birds or bats (Fig. 2). In addition, in birds, unlike bats, uric acid and feces mix directly in the cloaca before excretion. Finally, a subset of our bird samples was collected via intestinal aspiration, prior to the addition of uric acid in the cloaca, and analyses performed with this subset of samples did not differ substantially from those performed on feces (MRM fecal: R2 = 0.02, P = 0.001 for host phylogeny, P = 0.778 for diet; intestinal: R2 = 0.02, P = 0.001 for host phylogeny, P = 0.123 for diet). Consequently, this technical artifact explanation appears unlikely.

15The observed bird-bat similarities may also be due to the captivity state of the source animals, in that much higher proportions of bird and bat samples were from wild individuals than those of nonbat mammals. However, we found that captivity only explains a small amount of the variation in our data (see Table S2), consistent with a previous study showing that the effects of captivity on the animal gut microbiome is not strong across mammals but rather varies markedly across host orders. Furthermore, a recent study of strictly wild animals also found higher levels of phylosymbiosis in mammals than in other vertebrate groups, including birds, consistent with captivity being an unlikely explanation for our observations.

16Our data also suggest insectivory could be playing a role in microbiome specificity, at least in mammals. Although diet explained only a very small proportion of variance in bird microbiomes, and insectivorous mammals were not tightly clustered in our analysis (Fig. S3B), insectivorous mammals harbored some of the least mammalian-specific microbes (Fig. 4) and, like bats and birds, harbored higher proportions of Proteobacteria (Fig. S1). Most individual lineages of insectivores were not sampled deeply enough to accurately determine the strength of phylosymbiosis, but insect-eaters like armadillos (Cingulata) did have some of the weakest signals in our data set (Fig. 5). Previous work has reported some degree of convergence in the gut microbiomes of insectivorous mammals. Genomic analysis of insectivorous mammals also indicates convergent retention of functional chitinase genes expected to play a role in digesting the exoskeletons of their prey, potentially indicating a shift in reliance away from microbial and toward host function. If this is the case, the major shifts in host specificity we observe along independent branches on the vertebrate phylogeny in our data set could potentially be used to identify parallel changes in host genes mediating gut microbiome diversity.

17Collectively, our results suggest a striking convergent loss of host-microbe associations with the evolution of powered flight, breaking a conserved and consistent pattern of phylosymbiosis observed in nonflighted mammals. We propose that this loss of microbiome specificity represents not just a passive shift in neutral exposure to microbial diversity but rather a fundamental and convergent shift in the physiological mechanisms responsible for maintaining host-specific gut microbiomes. Though these mechanisms are not yet well understood, the hypothesis of convergence between bats and birds offers novel opportunities for their discovery and validation via comparative approaches. Flight need not be the only mediator of such decreased reliance; indeed, the other mammalian taxa shown here to have relatively less-specific microbiomes offer additional opportunities for such comparisons. Our findings also raise a question not often voiced: what are the evolutionary and metabolic costs of maintaining a specific microbiome?

18Through a large collaborative network of researchers and zoo directors, samples were collated from museum collections, zoos, reserves, and populations in their natural habitats (see Data Set S1 in the supplemental material). Collections were conducted under approved IACUC protocols, and appropriate permits were obtained for sample collection and export where necessary. Sampling of fecal material from animals in zoos and reserves involved collecting up to 2 g of fresh fecal material using sterile swabs (BD CultureSwab or equivalent) within minutes to hours of deposition, which was immediately frozen, and stored at −20°C until DNA extraction. Samples from museum collections included fecal or intestinal contents subsampled from frozen specimens using sterile implements or swabs. The 214 Smithsonian specimens were aspirations of large intestine (colon) from freshly killed wild birds. Samples taken directly from wild populations included fecal contents stored in RNAlater, on Whatman FTA cards, or in 95% ethanol (EtOH) (Data Set S1).

19Samples collected for this project were processed using the Earth Microbiome Project standard processing protocols. A subset of samples originally collected for previously published studies were reprocessed using these protocols. Briefly, sterile cotton-tipped swabs were used to transfer approximately 50 mg of sample to 96-well PowerSoil PowerMag DNA extraction plates (Qiagen), which were homogenized using a TissueLyser beadbeater (Qiagen); extractions were completed on a KingFisher magnetic bead transfer robot (Thermo Fisher). From eluted DNA, triplicate PCRs using the 515f/806r EMP primers amplified the V4 region of the 16S rRNA gene, and pooled amplicons were sequenced on Illumina MiSeq and HiSeq instruments.

20To ensure that data from new samples were processed consistently with samples from prior studies, we uploaded all new sequence data to the Qiita web-based microbiome analysis platform for initial processing. All sequences were demultiplexed and quality filtered using Qiita defaults, and forward reads were trimmed to 100 bp prior to processing with Deblur to remove sequencing errors. To avoid potential artifactual sequences, we used the “positive-filtered” output table from Deblur, which retains only those sequences that approximately match the Greengenes 13_8 16S reference database.

21Deblurred ASV tables and sample metadata were downloaded from each study and further processed using Qiime2. Per-study ASV tables were combined, and a phylogeny was estimated with SEPP as implemented in the q2-fragment-insertion Qiime2 plugin. Taxonomy was assigned using the Qiime2 naive Bayes feature classifier trained against the Greengenes 13_8 reference. ASVs classified as from mitochondria or chloroplasts were excluded from further analysis. Samples were rarefied to 10,000 reads. To prevent large imbalances in sample number among species from influencing ordinations, we randomly subsampled the rarefied data set to a maximum of 5 individuals per species, and these tables were used throughout the analysis.

22Metadata tables for all samples were combined, and host-level metadata were added from several sources. To accomplish this, we first curated host species names and taxonomies against the NCBI taxonomy database using the Taxon Names Resolver python package. Samples for which the provided species name could not be automatically resolved were manually curated, as were taxa for which individual taxonomic levels were missing from NCBI taxonomy (e.g., Cetartiodactyla, which is an unranked level in NCBI). Curated species names were then matched against the EltonTraits ecological trait database, and any missing taxa checked and manually curated if necessary. Metadata on bird migration were compiled by K.W. from sources cited in Data Set S1.

23To obtain an approximate host phylogeny, we matched curated species names to the TimeTree database. Unmatched taxa were manually curated and matched to the corresponding species binomial in the TimeTree database where possible. For taxa not present in the TimeTree database, we substituted a close congeneric species if (i) one was present in the TimeTree database, and (ii) there were no other congeneric species present in our sample set.

24Beta diversity measures were calculated for the complete data set in Qiime2, using the SEPP insertion phylogeny for UniFrac phylogenetic metrics. Principal-coordinate analyses were performed in Qiime2 and visualized using Emperor and ggplot2 in R. Tests for categorical differences in beta diversity were performed using PERMANOVA as implemented in R’s vegan package. To compare differences in beta diversity to differences in host diet and phylogeny, we used multiple regression of matrices as implemented in the ecodist package in R. To represent host evolutionary distance in these regressions, we used patristic distances derived from the TimeTree host phylogeny. To represent dietary dissimilarities, we used Bray-Curtis distances derived from the EltonTraits quantitative dietary compositions. To assess how correlations between microbiome dissimilarity and host phylogenetic and dietary dissimilarities were sensitive to microbial phylogenetic resolution, we used the SEPP insertion tree as an estimate of microbial phylogeny and compared correlations at different bacterial phylogenetic resolutions to the above host dietary and phylogenetic distances using Mantel tests in the beta diversity through time algorithm.

25Alpha diversity measures (observed OTUs, Faith’s phylogenetic diversity, and Shannon diversity) were calculated from the ASV tables rarefied to 10,000 sequences per sample in Qiime2. The SEPP insertion tree was used as an estimate of the bacterial phylogeny for the phylogenetic diversity measure. Differences in alpha diversity by categorical metadata variables were tested using ANOVA in R, with pairwise differences between categories assessed with the multcomp package. Differences in alpha diversity by continuous metadata variables were assessed using linear regressions.

26We used several methods to assess host specificity in microbiomes. One measure of specificity is phylosymbiosis, or the correlation between host phylogenetic distance and microbiome dissimilarity. To quantify phylosymbiosis systematically across taxa, we implemented a version of the Mantel test as used by Nishida and Ochman, comparing measures of microbial beta diversity to host patristic distances derived from the TimeTree time-calibrated phylogeny. We used the Jaccard dissimilarity metric for most analyses, as this measure reflects the proportion of shared exact ASVs, and so is more sensitive to recent microbial codiversification than UniFrac. To prevent the possibility of zero-length within-species branches from biasing regressions, we randomly subsampled the data set to a single representative per species; for visualization purposes (e.g., Fig. 5), we plotted all between-species points while excluding within-species points. We calculated the Mantel Pearson correlation between microbiome dissimilarity and host phylogenetic distance and the partial Mantel correlation (conditioned on dietary distances) at every node of the host phylogeny that contained at least seven tips using the EcoPy package. We visualized the strength of the correlation across the host phylogeny using the iToL web-based tree rendering tool (Fig. 1).

27We also tested the specificity of ASVs to particular host taxa using an implementation of the environmental entropy calculation in reference. Briefly, we calculated the host taxonomic distribution of each ASV, represented by a vector of integers corresponding to the number of samples from each host taxon in which the ASV was observed. We only considered ASVs observed in at least ten samples. We transformed these integer vectors to per-ASV proportional values and visually illustrated host taxonomic distributions for ASVs with stacked bar charts (Fig. 3). We quantified specificity of these host taxonomic distributions using Pielous’s evenness statistic: maximally specific ASVs are found in just one host taxon, while minimally specific ASVs are found evenly across all host taxa. Because our sample set was not perfectly balanced across host taxa, biases in sampling could make ASV distribution appear uneven simply due to uneven sampling. Thus, we calculated a standardized effect size (SES) for this statistic. We permuted sample assignment within each ASV in the observation table, such that each ASV was found in the same number of samples in the permuted table as in the raw table. We then calculated the SES as the raw evenness statistic for an ASV, minus the mean of 100 permuted statistics, divided by the standard deviation of the permuted statistics. We calculated a separate Pielou’s evenness SES value for each ASV based on its distribution across mammalian order-level taxa and across avian order-level taxa.

28To assess the degree to which individual ASVs were specific to mammalian or avian taxa, we calculated their likelihood of being observed in mammalian or avian samples, respectively. To do this, we used log risk ratios, or the log of the ratio of an ASV’s prevalence in mammalian samples to its prevalence in avian samples. Because this is undefined for ASVs never found in either mammals or birds, we set values for these taxa to ±3, which just exceeded the range of defined values in our data set. We then calculated an average per sample as a measure of how specific ASVs in each sample were to mammals or birds.

29To predict high level microbial phenotypic characteristics using 16S sequence data, we implemented BugBase (https://bugbase.cs.umn.edu/). We first clustered the deblurred ASVs against the Greengenes reference database (v13_5) at 99% identity and then rarefied the data to 5,000 sequences per sample.

30To rank microbial phyla based on their association with birds, we used a regression approach for detecting differential abundances in microbiome data. ASVs were first collapsed at the phylum level based on taxonomy assigned using the Qiime2 naive Bayes feature classifier trained against the Greengenes 13_8 reference as described above. A model was built testing for differences among host classes, with Mammalia serving as the reference, using a batch size of 10 and an epoch of 1,000,000.

31Sequence data and metadata tables are available without restriction in Qiita (https://qiita.ucsd.edu/study/description/11166; full list of study identifiers [IDs] is in Data Set S1) and EBI (accession no. PRJEB35449). Analysis notebooks are available on Github (https://github.com/tanaes/tetrapod_microbiome_analysis).