Decision letter: Community diversity is associated with intra-species genetic diversity and gene loss in the human gut microbiome

Article Figures and data Abstract Editor's evaluation Introduction Results Discussion Methods Data availability References Decision letter Author response Article and author information Metrics Abstract How the ecological process of community assembly interacts with intra-species diversity and evolutionary change is a longstanding question. Two contrasting hypotheses have been proposed: Diversity Begets Diversity (DBD), in which taxa tend to become more diverse in already diverse communities, and Ecological Controls (EC), in which higher community diversity impedes diversification. Previously, using 16S rRNA gene amplicon data across a range of microbiomes, we showed a generally positive relationship between taxa diversity and community diversity at higher taxonomic levels, consistent with the predictions of DBD (Madi et al., 2020). However, this positive 'diversity slope' plateaus at high levels of community diversity. Here we show that this general pattern holds at much finer genetic resolution, by analyzing intra-species strain and nucleotide variation in static and temporally sampled metagenomes from the human gut microbiome. Consistent with DBD, both intra-species polymorphism and strain number were positively correlated with community Shannon diversity. Shannon diversity is also predictive of increases in polymorphism over time scales up to ~4-6 months, after which the diversity slope flattens and becomes negative – consistent with DBD eventually giving way to EC. Finally, we show that higher community diversity predicts gene loss at a future time point. This observation is broadly consistent with the Black Queen Hypothesis, which posits that genes with functions provided by the community are less likely to be retained in a focal species' genome. Together, our results show that a mixture of DBD, EC, and Black Queen may operate simultaneously in the human gut microbiome, adding to a growing body of evidence that these eco-evolutionary processes are key drivers of biodiversity and ecosystem function. Editor's evaluation This paper analyses meta-genomic human gut microbiome data to understand how biodiversity arises and can be maintained. It makes an important contribution by strengthening the diversity-begets-diversity hypothesis and linking it to signatures of gene loss expected from the Black Queen hypothesis. While only correlative data is used to draw conclusions, the methods are solid and alternative hypotheses are clearly outlined. https://doi.org/10.7554/eLife.78530.sa0 Decision letter Reviews on Sciety eLife's review process Introduction Our understanding of microbial evolution and diversification has been enriched by experimental studies of bacterial isolates in the laboratory, but it remains a challenge to study evolution in the context of more complex communities (Lenski, 2017). Ongoing advances in culture-independent technologies have allowed us to study bacteria in the complex and dense communities in which they naturally occur (Garud and Pollard, 2020). Within a community, individual players engage in many negative and positive ecological interactions. Negative interactions can originate from competition for resources and biomolecular warfare (Hibbing et al., 2010; Mitri and Foster, 2013), while positive interactions can stem from secreted metabolites that are used by other members of the community (cross-feeding) (Venturelli et al., 2018). These ecological interactions can create new niches and selective pressures, leading to eco-evolutionary feedbacks whose nature are yet to be fully understood. Ecological interactions can yield positive or negative effects on the diversification of a focal species. Under the 'diversity begets diversity' (DBD) hypothesis, higher levels of community diversity increase the rate of speciation (or diversification, more generally) due to positive feedback mechanisms such as niche construction (Calcagno et al., 2017; Schluter and Pennell, 2017). Competition for limited niche space could also drive DBD if species diversify into new niches to avoid competition (Meyer and Kassen, 2007; Mitri and Foster, 2013; Schluter, 2000). By contrast, the 'ecological controls' (EC) hypothesis posits that competition for a limited number of niches at high levels of community diversity results in a negative effect on further diversification. Metabolic models predict that DBD may initially spur diversification due to cross-feeding, but the diversification rate eventually slows and reaches a plateau as metabolic niches are filled (San Roman and Wagner, 2021). These theoretical predictions are largely supported by our previous study involving 16S rRNA gene amplicon sequencing data from the Earth Microbiome Project, in which we observed a generally positive relationship (which we call the diversity slope; Figure 1) between community diversity and focal-taxon diversity at most taxonomic levels, reaching a plateau at the highest levels of diversity (Madi et al., 2020). Figure 1 Download asset Open asset Diversity begets diversity (DBD) and ecological controls (EC) hypotheses illustrated. Hypothetical microbial communities are illustrated as gray circles containing assemblages of microbial species, shown in different colors. 'DBD' means that the focal species is more likely to acquire diversity – through de novo mutation, invasion of a different strain of the same species, or a combination of both – in a community with high diversity. This is because new niches are created in a more diverse community. By contrast, 'EC' means that the focal species is more likely to acquire diversity through strain invasion or mutation in a community with low diversity. This is because niches remain unfilled in a low-diversity community, while niche space is saturated in a high-diversity community, impeding further diversification. In this previous study, we found stronger support for DBD in the animal gut relative to more diverse microbiomes such as soils and sediments, which were closer to a plateau of diversity (Madi et al., 2020). While diversity slopes were generally positive at taxonomic levels as fine as amplicon sequence variants (akin to species or strains) within a genus, they were most positive at higher levels such as classes or phyla. A recent experiment on soil bacteria also found evidence of DBD at the family level, likely driven by niche construction and metabolic cross-feeding (Estrela et al., 2022). It therefore remains unclear if the predictions of DBD hold primarily at these higher taxonomic levels, involving the ecological process of community assembly, or if they also apply at the finer intra-species level. Within-host intra-species diversity can arise by co-colonization of a host by genetically distinct strains belonging to the same species or evolutionary diversification of a lineage via de novo mutation and gene gain/loss events within a host. Such fine-scale strain-level variation has important functional and ecological consequences; among other things, strains are known to engage in interactions that cannot be predicted from their species identity alone (Goyal et al., 2022). Although closely-related bacteria are expected to have broadly similar niche preferences, finer-scale niches may differ below the species level (Martiny et al., 2015). For example, the acquisition of a carbohydrate-active enzyme by Bacteroides plebeius allows it to exploit a new dietary niche in the guts of people consuming nori (seaweed) (Hehemann et al., 2010), and single nucleotide adaptations permit Enterococcus gallinarum translocation across the intestinal barrier resulting in inflammation (Yang et al., 2022). Despite their potential phenotypic effects, it is unknown if such fine-scale genetic changes are favored by higher community diversity (due, for example, to niche construction, as predicted by DBD) or suppressed (due to competition for limited niche space, as predicted by EC). Competition could also lead to DBD if focal species evolve new niche preferences to avoid extinction (Mitri and Foster, 2013; Schluter, 2000) – an idea with some support in experimental microcosms (Meyer and Kassen, 2007) but largely unexplored in natural communities. Here, we investigate the relationship between intra-species genetic diversity and community diversity in the human gut microbiome, a well-studied system in which we previously found support for DBD at higher taxonomic levels. We use static and temporal shotgun metagenomic data from a large panel of healthy adult hosts from the Human Microbiome Project (HMP) (Lloyd-Price et al., 2017; Human Microbiome Project Consortium, 2012) as well as from four healthy individuals sampled almost daily over the course of 1 year (Poyet et al., 2019). Using metagenomic data allows us to track change in single nucleotide variation, strain diversity, and gene gain or loss events within relatively abundant species in the microbiome, and to study how these measures of intra-species diversity are associated with community diversity. Although such analyses of natural diversity cannot fully control for unmeasured confounding environmental factors, they are an important complement to controlled experimental and theoretical studies which lack real-world complexity. Results We investigated the relationship between community diversity and within-species genetic diversity in human gut microbiota using two shotgun metagenomic datasets. First, we analyzed data from a panel of 249 healthy hosts (Lloyd-Price et al., 2017; Human Microbiome Project Consortium, 2012), in which stool samples were collected one to three times from each host at approximately 6-month intervals. Second, we analyzed data from four individuals sampled more densely over the course of ~18 months (Poyet et al., 2019). In both cases, we only consider intra-species diversity of relatively abundant species that are well sampled in these metagenomic datasets (Methods). We examined several metrics of community diversity and intra-species diversity and calculated the slope of their relationship, defined as the diversity slope (Figure 1). We note that intra-species diversity can arise within hosts via de novo point mutation, gene gain or loss, or the coexistence of genetically distinct strains that diverged before colonizing the host. To quantify community diversity, we calculated Shannon diversity and richness at the species level. Shannon diversity is relatively insensitive to sampling effort (Madi et al., 2020; Walters and Martiny, 2020) but richness can be underestimated in low sample sizes. We therefore computed richness on data rarefied to an equal number of reads per sample, yielding generally similar results to unrarefied data (described below). In all cases, we included the number of reads per sample (coverage) as a covariate in our models, as this could affect estimates of both community diversity and intra-species diversity. To quantify intra-species diversity, we used a reference genome-based approach to call single nucleotide variants (SNVs) and gene copy number variants (CNVs) within each focal species and computed polymorphism rates, measured as the fraction of synonymous nucleotide sites in a species' core genome with intermediate allele frequencies (between 0.2 and 0.8) within a host (Methods). We also repeated the analysis on nonsynonymous sites, as these are subject to stronger selective constraints. As an additional metric of intra-species diversity, we inferred the number of strains within each species using StrainFinder applied to all polymorphic sites (including those outside the 0.2–0.8 frequency range) (Smillie et al., 2018). Community diversity is positively associated with intra-species polymorphism in the human gut microbiome As an exploratory visualization, we began by plotting the relationship between community diversity and intra-species polymorphism rate calculated at synonymous sites in cross-sectional HMP metagenomes for the nine most prevalent species (Figure 2A and B). The slope of this relationship (the diversity slope; Figure 1) provides an indicator of the evidence for DBD (positive slope) or EC (flat or negative slope). The relationship between polymorphism rate and community diversity was mostly positive in the top nine most prevalent species in HMP hosts (Figure 2A and B). These nine species are used as a simple illustration of the diversity slope, not as a formal hypothesis-testing framework. Figure 2 with 2 supplements see all Download asset Open asset Positive association between community diversity and within-species polymorphism in cross-sectional Human Microbiome Project (HMP) samples. (A) Scatter plots showing the relationship between community Shannon diversity and within-species polymorphism rate (estimated at synonymous sites) in the nine most prevalent species in HMP. (B) Scatter plots showing the relationship between species richness and within-species polymorphism rate in the nine most prevalent species in HMP. These are simple correlations to show the relationships in the raw data. Significant correlations are shown with red trendlines (Spearman correlation, p<0.05); non-significant trendlines are in gray. Results of generalized additive models (GAMs) predicting polymorphism rate in a focal species as a function of (C) Shannon diversity, (D) species richness estimated on all sequence data, and (E) species richness estimated on rarefied sequence data. GAMs are based on data from 69 bacterial species across 249 HMP stool donors. Adjusted R2 and Chi-square p-values corresponding to the predictor effect are displayed in each panel. Shaded areas show the 95% confidence interval of each model prediction. See Supplementary file 1a and Supplementary file 2 section 1 for detailed model outputs. To generalize across species and to formally test the predictions of DBD, we fit generalized additive models (GAMs) to the HMP data. Using GAMs, we are able to model non-linear relationships and account for random variation in the strength of the diversity slope across bacterial species, the uneven number of samples per host, and the non-independence of samples from the same host (Methods; see Supplementary file 1a and Supplementary file 2 section 1 for additional model details). These GAMs included 69 focal species with sufficient coverage to quantify within-species polymorphism (Methods); the results therefore apply to relatively abundant species in the human gut microbiome. GAMs showed an overall positive association between within-species polymorphism and Shannon diversity (Figure 2C, GAM, p=0.031, Chi-square test) as well as between within-species polymorphism and community richness after controlling for coverage as a covariate (Figure 2D, GAM, p=0.017, Chi-square test) or rarefying samples to an equal number of reads (Figure 2E, GAM, p=2.63e-04, Chi-square test). The random effect of species identity is highly significant in all models, indicating that each bacterial species has its own characteristic diversity slope (Supplementary file 1a). It appears that synonymous polymorphism reaches a plateau at high levels of community richness, which is particularly evident when using rarefied data (Figure 2E). Using the same GAMs applied to nonsynonymous polymorphism, we found no significant associations between diversity and within-species polymorphism rate (GAM, p>0.05, Chi-square test) (Supplementary file 1b, Supplementary file 2 section 4). This could be due to lower statistical power, since there are fewer nonsynonymous than synonymous sites, or it could reflect a true difference in the diversity slope between these site categories. These generally positive correlations between focal species polymorphism and species-level measures of community diversity also hold when community diversity is measured at higher taxonomic levels; specifically, synonymous polymorphism rate was significantly positively associated with Shannon diversity calculated at the genus and family levels (GAMs, p<0.05, Chi-square test) (Figure 2—figure supplement 1, Supplementary file 1c). However, synonymous polymorphism rate was not significantly associated with Shannon diversity calculated at the highest taxonomic levels (order, class, and phylum, GAMs, p>0.05, Chi-square test). The positive correlation between polymorphism rate and richness held at all taxonomic levels (GAMs, p<0.05, Chi-square test) (Figure 2—figure supplement 1, Supplementary file 1c, Supplementary file 2 sections 2 and 3). When estimated at nonsynonymous sites, polymorphism rate was not significantly correlated with Shannon diversity at any taxonomic level (GAMs, p>0.05, Chi-square test), but was positively correlated with richness at the highest levels (phyla, class, and order, p=3e-04, p=0.017, and p=6.11e-04, respectively, Chi-square test from GAMs) (Figure 2—figure supplement 2, Supplementary file 1d, Supplementary file 2 sections 5 and 6). Even when not statistically significant, the diversity slopes were generally positive at all taxonomic levels for both synonymous and nonsynonymous polymorphism (Figure 2—figure supplements 1 and 2). Overall, these results are consistent with the predictions of DBD at most taxonomic levels. However, slightly different relationships are observed when considering different measures of community diversity (Shannon or richness) and different components of within-species diversity (nonsynonymous or synonymous). Different measures of community diversity have contrasting associations with intra-species strain diversity Within host polymorphism rates span several orders of magnitude (10–5/bp to 10–2/bp), largely due to the fact that strain content is variable across hosts. As previously argued (Garud et al., 2019), with conservatively high estimates for mutation rate (μ~10−9) (Sung et al., 2012), generation times (~10/day) (Poulsen et al., 1995), and time since colonization (<100 years), polymorphism rates of ~10–2/bp or more are inconsistent with within-host diversification of a single colonizing lineage. Therefore, hosts with relatively high intra-host polymorphism rates are likely colonized by mixtures of multiple strains that diverged long before colonizing a host. Moreover, recent work suggests that the numbers and genetic composition of strains colonizing a host can vary from host to host (Garud et al., 2019; Olm et al., 2017; Russell and Cavanaugh, 2017; Truong et al., 2017; Verster et al., 2017). The associations between polymorphism and community diversity (Figure 2) are likely driven by a combination of de novo mutation and co-colonization by multiple strains. To separate these two sources of diversity and to explicitly account for the strain structure within hosts, we inferred the number of strains per focal species with StrainFinder (Smillie et al., 2018) (Methods) and used strain number as another quantifier of intra-species diversity. The relationship between community diversity and strain number varied depending on the focal species and the measure of community diversity. For example, the inferred number of Bacteroides vulgatus strains increased with community diversity, while Bacteroides uniformis strain count decreased or remained flat (Figure 3A and B). Expanding upon these examples, we used generalized linear mixed models (GLMMs) to investigate the relationship between the number of strains per focal species and community diversity, while taking into account coverage per sample as a covariate and variation between species, hosts and samples as random effects (Methods). GLMMs are a special case of GAMs that can handle overdispersed, zero-truncated count data such as strain counts. The number of strains per focal species was positively correlated with community Shannon diversity (GLMM, p=3.58e-07, likelihood ratio test [LRT]) (Figure 3C, Supplementary file 1e, Supplementary file 2 section 7.1). This suggests that the positive correlation between polymorphism rate and Shannon diversity (Figure 2) is due at least in part to strain diversity. Figure 3 with 1 supplement see all Download asset Open asset Associations between community diversity and strain number in cross-sectional Human Microbiome Project (HMP) samples. (A) Scatter plots showing the relationship between Shannon diversity and the inferred number of strains within each of the nine most prevalent species in HMP. (B) Scatter plots showing the relationship between species richness and the inferred number of strains within each of the nine most prevalent species in HMP. Significant linear correlations are shown with red trendlines (Pearson correlation, p<0.05); non-significant trend lines are in gray. Results of generalized linear mixed models (GLMMs) predicting strain count in a focal species as a function of (C) Shannon diversity, (D) species richness estimated on all data, and (E) species richness estimated on rarefied sequence data. Diversity estimates (X-axis) are standardized to zero mean and unit variance in the models. The Y-axis shows the mean number of strains per focal species predicted by the GLMM. GLMMs are based on data from 184 bacterial species across 249 HMP stool donors. p-Values (likelihood ratio test) are displayed in each panel. Shaded areas show the 95% confidence interval of each model prediction. See Supplementary file 1e and Supplementary file 2 section 7 for detailed model outputs. By contrast, species richness was negatively correlated with strain number (GLMM, p=1.67e-06, LRT) (Figure 3D, Supplementary file 1e, Supplementary file 2 section 7.2). The negative relationship with richness was unlikely to be confounded by sequencing depth, since the same result was obtained using rarefied data (Figure 3E, Supplementary file 1e, Supplementary file 2 section 7.3). The negative strain number-richness relationship also held at all other taxonomic ranks (GLMM, p<0.05, LRT), while the strain number-Shannon diversity relationship was generally positive (Figure 3—figure supplement 1, Supplementary file 1f, Supplementary file 2 sections 8 and 9). These effects also appear to be species-specific: for example, the number of B. vulgatus strains per host is positively correlated with both Shannon diversity and richness (consistent with DBD predictions) whereas Bacteroides ovatus has no relationship with Shannon diversity but a negative correlation with richness (consistent with EC; Figure 2A and B). Together, these results reveal that different components of community diversity can have contrasting effects on the diversity slope. Community Shannon diversity is a predictor of intra-species polymorphism and gene loss in time data Our analyses have only individual time which static of the processes of community assembly and evolution in the microbiome. To these over we analyzed HMP were sampled two to three times months Under a DBD we community diversity at an time point to result in higher within-species polymorphism at a future time point. To test this we defined as the difference between polymorphism rates at the two time (Methods). We also investigated the effects of community diversity on gene loss and gain events within a focal species, as such changes in gene content are known to occur within host gut microbiomes (Garud et al., 2019; et al., and 2020; et al., 2019). Here, a gene was if its copy number was and if As in the cross-sectional analyses we also controlled for sequencing of the sample and genes with coverage or in multiple species (Methods). In HMP polymorphism change showed no significant relationships with community diversity at the time it was estimated with Shannon or species richness (GAM, (Supplementary file 2 section These results that DBD is or over time in the human By contrast, we found that gene loss in a focal species between two time was positively correlated with community diversity at the time point (Figure and for richness and rarefied richness, (Supplementary file Supplementary file 2 section not show any significant relationships with community diversity (GLMM, for gene loss in more diverse communities is a of the Black Queen hypothesis provided that higher community diversity results in more gene functions that for in a focal species et al., species in HMP samples fewer than genes over months – consistent with de novo events of a genes – but of genes were from a host, that strains with were in more diverse communities (Figure and B). Figure Download asset Open asset Positive association between community diversity and gene loss in Human Microbiome Project (HMP) time (A) Scatter plots showing the relationship between Shannon diversity at time point 1 and gene loss between and within each of the nine most prevalent species in HMP. (B) Scatter plots showing the relationship between species richness at and gene loss between and within each of the nine most prevalent species in HMP. Significant linear correlations are shown with red trendlines (Pearson correlation, p<0.05); non-significant trend lines are in gray. The Y-axis is on a for Results of generalized linear mixed models (GLMMs) predicting gene loss in a focal species as a function of (C) Shannon diversity, (D) species richness estimated on all data, and (E) species richness estimated on rarefied sequence data. p-Values (likelihood ratio test) are displayed in each panel. Shaded areas show the 95% confidence interval of each model prediction. The Y-axis is on the which to for negative GLMMs with a count GLMMs are based on data from bacterial species across HMP stool sampled at more than one time point. See Supplementary file and Supplementary file 2 section for detailed model outputs. To study these at higher temporal resolution, we analyzed shotgun metagenomic data from four more sampled healthy individuals from a previous study (Poyet et al., 2019). from was over months with a of 1 between an over months 2 between over 5 months 1 between and over 7 months 2 between In this data, we both polymorphism change and gene and between two time in species with a gene coverage of in at least samples. These species of two two two as well as and Using the we community diversity in the gut microbiome at one time point could predict polymorphism change at a future time point by GAMs with the change in polymorphism rate as a function of the between community diversity at the time point and the number of between the two time Shannon diversity at the time point was correlated with increases in polymorphism (consistent with DBD) up to into the future (Figure supplement but this relationship and (consistent with at time (Figure Supplementary file GAM, Chi-square test). The diversity slope is approximately flat for time between and months, which could no significant relationship was found in samples were collected relationship was observed between community richness and changes in polymorphism (Supplementary file GAM, Figure 5 with 1 supplement see all Download asset Open asset Community diversity is associated with increases in focal species polymorphism over time and gene loss in dense gut microbiome time (A) Results of a generalized additive model predicting polymorphism change in a focal species as a function of the between Shannon diversity at the time point and the time between two time in data from et The response was in the Results of generalized linear mixed models (GLMMs) predicting (B) number of genes and (C) number of genes between two time in a focal species as a function of the between Shannon diversity at the time point and the time between the two time (D) Results of the predicting the number of genes in a focal species as a function of the between rarefied species richness at the time point and the time between the two time The illustrated time to the the and the See Supplementary file and and Supplementary file 2 section for detailed model outputs. These analyses are based on data from bacterial species across four stool from et statistically significant relationships are relationships are not the predicting polymorphism