Distinct, stable communities assemble in microcosms
The aquatic communities from 10 individual Sarracenia purpurea pitchers were filtered to focus on bacteria and inoculated into a realistic, complex nutrient source: sterilized, ground crickets in acidified water. The in vitro communities were serially transferred every 3 days for 21 transfers, using a low dilution rate of one-part culture to one-part fresh media. We repeated the process with the unfiltered communities as well. Community composition was measured for each transfer using 16S rRNA sequencing (see Methods). From ~8 million sequences across all microcosms and timepoints, DADA2 analysis inferred 889 distinct ASVs (Amplicon Sequence Variants, which we treat as units of diversity). The most abundant phylum was Proteobacteria, followed by Firmicutes and Bacteroidetes. The top twenty ASVs are displayed in Fig. 1a, accounting for 69.4% of the reads; for ASVs without assigned genera, genus names were recovered from the full 16S rRNA genes of our cultured strains that matched 100% with ASV sequences. The most abundant genera (Aquitalea, Pseudomonas, Achromobacter, Comamonas, and Delftia), all contain bacterial species known to live in freshwater, soil, or plant-associated habitats29,30,31,32,33. Using DNA concentrations as a proxy to measure biomass increase, we found that there is a ~10-fold increase in biomass during in vitro assembly (Fig. 1a). Microcosms M03 and M09 stood out as those with the highest biomass yield as well as having the lowest diversity.
a Relative abundances of the top 20 Amplicon Sequence Variants (ASVs) across the ten microcosms during the course of the serial transfer experiment. ASVs are listed one time each on the bar plot, with taxonomic classification in the legend. The DNA concentrations for each timepoint are graphed above as points, with a Loess fit as a solid line. b Communities change quickly and then stabilize in the two-dimensional Non-metric Multidimensional Scaling (NMDS) plot of Bray–Curtis dissimilarities of the microcosm community compositions. The microcosm name is listed in black next to the Day 0 point for each microcosm, and the lines connect the timepoints. Colored numbers indicate the mean effective number of species for each community post Day 21. c The Bray–Curtis dissimilarity of ASV relative abundances between adjacent days decreases over the course of the experiment. The thick black line shows a Loess fit to all data points, and the thin line marks Day 21. d At the end of the experiment, most ASVs were present in only one microcosm (top), but ASVs present in more microcosms tended to have higher mean relative abundances (bottom), shown as black circles with ±1 standard deviation error bars.
Community assembly dynamics in our microcosms were dominated by both ASV loss and dramatic changes in ASV abundances. We assessed coarse-grained dynamics using the Bray–Curtis dissimilarity between subsequent time points, which exhibited two seemingly distinct phases (Fig. 1b, c): a first phase of rapid change, where many species went extinct and the survivors increased in frequency; and a second phase of slow changes that started after seven transfers (21 days, indicated by the thin line in Fig. 1c), in which extinctions were rare and likely caused by competitive interactions. In this last phase, communities slowly approached what seemed to be an equilibrium state (Fig. 1b, c). However, the near-equilibrated communities maintained significant differences both in terms of composition and diversity. For example, no microcosm had the same top three ASVs (Fig. 1a), the NMDS ordination showed generally non-overlapping points for the different microcosms (Fig. 1b), and, even after Day 21 the effective number of species (see Methods for calculation) still differed among the microcosm communities, ranging from ~6 to 16 (small colored numbers in Fig. 1b). Most ASVs were rare (~65% found in ≤2 microcosms) and had low mean relative abundance (<1%), but the few core ASVs found across nine microcosms had higher mean relative abundances (~10%, Fig. 1d). Community composition at the family level remained distinct across our 10 microcosms, and phylogenetic metrics of alpha and beta diversity showed the same dynamics as non-phylogenetic metrics, indicating that observed differences were not simply due to changes in closely related species (Supplementary Fig. 1). In summary, communities first approached convergence due to the fast loss of diversity (primarily along the first NMDS axis), but remain distinct due to historical contingencies (primarily along the second NMDS axis).
Interestingly, the differences in richness near-equilibrium were seemingly pre-determined at early stages of assembly. Communities lost many ASVs between days 0–3 during initial adjustment to the laboratory environment (Fig. 2a), but the richness measured from the first timepoint of the experiment (Day 3) was an excellent predictor of the richness at the end of the experiment (Day 63) in a linear model with R2 = 0.9008 and p < 0.0001 (Fig. 2b). Notably, richness in samples taken directly from the pitcher plant (Day 0) had no significant correlation with that of Day 63 (linear model, R2 = 0.1978, p = 0.1105), consistent with the notion that some ASVs in the pitcher plant fluid were either metabolically inactive or unable to grow in our experimental conditions. Thus, changes in community composition after only three days of adjusting to the lab environment propagated throughout the assembly process, suggesting that historical contingencies played a significant role in structuring communities.
a Richness over time for each microcosm community. b Richness on Day 3 (1st timepoint) is strongly correlated with richness on Day 63 (21st timepoint). Linear model: R2 = 0.9008, p = 1.714 × 10−5. c The probability of going extinct at transfer t. Colored lines are probability densities for individual microcosms. Black points are averages across microcosms, the black line is the maximum likelihood distribution with a common parameter across microcosms (see main text and Methods) given by the inverse mean extinction time. d Richness over time for each microcosm community, normalized by the richness on Day 3. The black line shows the exponential decay curve parametrized by the mean proportion of surviving ASVs (from b) and the common ASV extinction rate (from c).
We found that the effects of historical contingency on community assembly were remarkably consistent and reproducible, implying that the lack of convergence was not due to stochastic dynamics but to deterministic effects. We performed the same experiment on a second set of microcosms from the same inocula that were not subjected to initial filtering, but otherwise underwent the same transfer protocol and amplicon sequencing (Supplementary Fig. 2). These microcosms were very similar to their filtered counterparts in composition, DNA concentrations, community dynamics over time and the equilibration process (Supplementary Fig. 2a–c). Indeed, unfiltered and filtered microcosms were indistinguishable in an NMDS visualization of their assembly dynamics (Supplementary Fig. 2b). Furthermore, in the unfiltered samples we saw the same strong correlation between the richness on Day 3 and 63 (Supplementary Fig. 2d). Taken together, these results reinforce the notion that initial pitcher source environment and composition predetermine the future states of community assembly.
Given the strong correlation between initial and final richness, we wondered whether temporal dynamics of community equilibration were also correlated across microcosms. For example, the temporal dynamics of ASV loss could be driven by external factors such as transfer intervals and the dilution factor, in which case all microcosms should exhibit comparable dynamics. Alternatively, the dynamics might be driven by community context and thus specific to each microcosm. To answer this question, we measured the distribution of ASV extinction times across the microcosm communities (Fig. 2c). As a null model, we tested if the loss of ASVs is a random process where each ASV that is bound to go extinct in a given community context does so with a fixed probability p in each transfer. This assumption implies that the extinction time distribution is described by a geometric distribution, which indeed gave a good fit for all microcosms (Pearson’s chi-squared test for differences was not significant: p > 0.05 in all cases). To investigate if microcosms can be described by a common extinction time distribution, we used the Bayesian Information Criterion to compare two models: using either one parameter per microcosm or a single parameter describing all microcosms. Surprisingly, the single-parameter model was strongly favored (relative likelihood = 7.7 × 1010) indicating that, despite biotic differences in richness and composition across microcosms, there was a common rate of ASV extinction. The process of ASV loss in unfiltered communities was also well-described by a geometric distribution using a single parameter to describe all microcosms (relative likelihood = 1.5 × 109, Supplementary Fig. 2e). This “universal” dynamic of species loss was also revealed when studying relative richness, i.e. normalized by richness at Day 3 (Fig. 2d). After this normalization, all relative richness curves mapped well to our model’s maximum likelihood distribution and approached a common equilibrium relative richness level (~50% of the initial richness, regardless of its absolute value, Fig. 2d). Thus, about half of all initial ASVs within a given microcosm were doomed to go extinct at a random point during the assembly process, while the rest persisted indefinitely. Interestingly, this result has been predicted by theory. Models based on random interaction matrices indeed predict that that the richness of the community after assembly is a fixed proportion of the initial richness of the species pool34. Our work is the first empirical validation of this prediction.
Although microcosms display similar extinction dynamics, the same ASV may still behave differently depending on the community in which it is present. To study the context-dependency of each species, we focus on those shared ASVs, which although representing a small subset of the species list, accounted on average for ~90% of the communities’ relative abundances (Fig. 1d). To ask to what extent species dynamics were dependent on community context, we first looked at the extinction times and found that the majority of the shared ASVs dropped out of different microcosms at different times (Fig. 3a), often with large differences in their persistence time. For example, ASV681 (in Fig. 3b) dropped out by Day 12 in microcosm M08, but persisted at high relative abundance through the end of the experiment in microcosms M04, M06, and M10. Because microcosm compositional dynamics were very similar when starting from the same initial bacterial inoculum (correlation of unweighted UniFrac distances shown in Fig. 3c), we leveraged the unfiltered samples to investigate context dependent behavior of individual ASVs.
a The y-axis lists all ASVs shared by at least two microcosms. Colored points mark which day the ASV was lost from that particular microcosm. Points at the far right of the graph are ASVs that persisted to the end of the serial transfer experiment. Red stars mark the ASVs shown in b. b Three examples of ASV dynamics in the different microcosms over time, where ASVs were lost early in some microcosms but persisted until the end in others. c Unweighted UniFrac distances in filtered and unfiltered communities (repetitions) are strongly correlated, highlighting that community dynamics are reproducible given the same inoculum. Mantel test, r = 0.85, p = 1 × 10−4; Linear model, R2 = 0.78, p = 2.2 × 10−16. Samples from Days 0–21 are included here, to capture dynamics up until equilibration. Hexplot color indicates the density of values in each hexagon. d Probability density of correlation coefficients between the same ASVs in repetitions started from the same inocula, between the same ASVs in distinct microcosms, and between randomly chosen ASVs. Correlation coefficients are large when comparing replicate communities started from the same pitcher plant inocula. Between different microcosms, there is a moderate increase in the number of positive correlation coefficients relative to random pairs. e Mean Bray–Curtis dissimilarities between final community compositions is correlated with mean cosine distances between trajectories of the same ASV.
We focused on ASVs that persisted at least once past equilibration (Day 21) and found that the temporal trajectories of the same ASVs were strongly dependent on community composition. ASVs present in multiple microcosms that survived the early extinction phase were correlated ~40% of the time across microcosms and were more likely to share the same fate (survive vs. go extinct) than a null model (Supplementary Fig. 3). However, while the cross-microcosm correlations of ASVs were moderately higher than those of any two random ASVs (mean = 0.23 and mode = 0.32 vs. random pairs: mean = 0.003 and mode = 0.01, Fig. 3d), species identity alone was not predictive of their dynamics; the community context was far more important. Indeed, when we compared ASV dynamics between filtered and unfiltered microcosms started from the same inoculum, we found much higher correlations (mean = 0.55 and mode = 0.89, Fig. 3d), consistent with the notion that the dynamics in a given community context are deterministic. The same was true when comparing at the family level (Supplementary Fig. 4). More importantly, across all microcosms, similarity in community composition was positively correlated with the cosine similarity (see Methods for details) in ASV dynamics (intercept-free linear model, R2 = 0.63, Fig. 3e). This means that biotic context, and not species identity, was the main factor that explained the ASV dynamics across microcosms. Because community composition is the only factor that changes across microcosms, this result implies that biological interactions (e.g. competition, niche construction via secreted metabolites, etc.) are the main drivers of population dynamics.
ASV dynamics are highly correlated with functional dynamics
Functional redundancy is thought to be widespread in regional pools of bacteria35, and thus communities can have diverse compositions but converge to similar functional activity. Previous studies have generally found stronger convergence in function than in composition, and this was a likely outcome from our experiment as all microcosms experienced the same environment in terms of nutrients, temperature and light. In agreement with this, carbon dioxide production, as measured with the MicroResp system, was highly variable across the different microcosms for the first measurement (Day 0–3), but then quickly converged to a similar level (Fig. 4a). After Day 3, communities could not be reliably distinguished based on their CO2 output. The low variation among microcosms in percent CO2 suggests that the bacterial communities were respiring at about the same rate.
a Variance in percent CO2 among the ten microcosms over the course of the serial transfer experiment. b Two-dimensional NMDS plot of the Bray–Curtis dissimilarities of functional activity (substrate use) as measured by EcoPlates. The microcosm name is listed in black next to the Day 0 point, and the lines connect the timepoints. c Procrustes rotation of the composition NMDS plot with the function NMDS plot for Day 63 (Procrustes correlation = 0.8991, p = 0.001). d Plot of Bray–Curtis dissimilarities in composition between all measured timepoints across all microcosms compared with the dissimilarities in functional profile for the corresponding samples (Mantel test r = 0.6403, p = 1 × 10−4), squared to better illustrate the spread of points. e Endochitinase activity over time for the five microcosms that strains were cultured from (top row) and frequency of the ASVs over time that map to strains with measurable endochitinase activity (bottom row). Endochitinase activity is shown in units/mL for the strains/ASVs in the bottom row by a gradient from gray (low activity) to red (high activity).
In contrast with CO2 production and the expectation of functional convergence, the stabilized microcosm communities showed clear differences in substrate use, as measured across 31 substrates with Biolog EcoPlates (Supplementary Fig. 5, Fig. 4b, c). For example, M09 was the only community able to metabolize 2-Hydroxybenzoic acid (salicylic acid) but unable to metabolize itaconic acid. These particular substrates may not be relevant in the pitcher plant context, but differences in EcoPlate functional fingerprints suggest a likelihood of other functional differences between microcosms. When plotting an NMDS ordination of Bray–Curtis dissimilarities based on EcoPlate functional measures (Fig. 4b) we found a very similar pattern as with species composition: functional activity undergoes an initial large shift and becomes more similar, yet remains distinct across microcosms. A Procrustes test comparing the NMDS plots of composition and function at the end of the experiment recovers a strong and highly significant correlation: 0.8991, p = 0.001 (Fig. 4c). Furthermore, when comparing composition to function across all days and samples using Bray–Curtis dissimilarities, samples with similar composition generally also had similar functional activity (Fig. 4d), and are strongly correlated in a Mantel test (r = 0.6403, p = 1 × 10−4). The correlation is in fact higher when only comparing the final day’s measurements (r = 0.6907, p = 0.002).
To profile the hydrolytic activity of the community, we focused on the activity of chitinases—the enzymes that degrade chitin—since chitin is the main component of insect exoskeletons and is a key carbon and nitrogen source in the pitcher plant system. In two different cricket species, ~7% of the dry weight is chitin and another 5% is the chitin derivative, chitosan36,37, so the ability to degrade these compounds would likely increase bacterial fitness in our microcosms. The chitin hydrolysis rate of the community supernatant was also highly variable across microcosms, with M03 and M09 having the highest measures of endochitinase activity (Fig. 4e). In summary, microcosm communities begin with different compositions as a result of historical contingencies affecting individual pitchers. These compositions shift as bacterial communities are brought into a new laboratory environment, but remain influenced by their starting compositions (Fig. 5, part i). We expected to see functional convergence (ii), since all communities were grown for many generations in the same environment. This convergence was observed in core functions, such as respiration, but not on chitinolytic activity or substrate utilization profiles. These “auxiliary” functions (encoded by transporters and enzymes at the periphery of the metabolic network) were instead strongly correlated with community composition, and thereby also affected by historical contingencies (iii).
Microcosms with the same environmental conditions assemble communities with distinct equilibria due to historical contingencies (i). Core functions such as respiration converge in a common environment, in agreement with the common notion that function and taxonomy are decoupled (ii). However, our results show key functional differences in ‘auxiliary’ metabolic functions, such as chitinase activity or carbon source preference, are strongly coupled to community composition (iii).
We developed an isolate collection to learn more about the mechanisms underpinning the strong differences in substrate utilization and hydrolytic activity across microcosms. By testing the enzymatic activity and substrate metabolism of individual isolates, we hypothesized we could identify bacterial strains responsible for the corresponding community function. We successfully isolated 350 strains and sequenced their full 16S rRNA genes. Of these, 176 mapped with 100% identity to 33 different ASVs from the amplicon sequencing. For the five microcosms we cultured from, 14 out of the combined top 15 ASVs (in terms of relative abundance) on Day 63 matched perfectly with cultured strains, as did 5–7 of the top 10 ASVs per microcosm. Isolated strains accounted for 67-88% of each microcosm’s relative abundance on the final day of the experiment (M03: 88%; M05: 67%; M07: 69%; M09: 87%; and M10: 82%). Our cultured strains had broad representation across different taxonomic groups (Supplementary Data 1).
Consistent with our hypothesis, the chitinase activity of our cultured strains mirrored microcosm activity: strains from M03 and M09 had the highest endochitinase activity out of all measured strains. The activity of individual strains/ASVs mapped well to the activity of the entire microcosms, where M03 and M09 were also highest (Fig. 4e). In M03 only one of the measured strains showed high activity (ASV589, a Chromobacterium species), while in M09 at least three strains had high endochitinase activity (ASV018 Chromobacterium, ASV842 Dyella, and ASV863 Burkholderia). The cumulative endochitinase activity across strains (strain activity multiplied by corresponding ASV frequency over time) correlated with microcosm activity across all microcosms (linear model R2 = 0.2979, p < 0.0001), and within M03, M07 and M09 (R2 = 0.9303, p < 0.0001, R2 = 0.7044, p = 0.0029, and R2 = 0.5142, p = 0.0078, respectively. Supplementary Fig. 6). Interestingly, the strains with the highest chitinase activity also had high protease activity, and all except for ASV842 had high lipase activity (Supplementary Fig. 7), suggesting that these bacteria are generally specialized to degrade insect tissue.
To ask to what extent the pattern of substrate utilization of the community could be reduced to the substrate utilization of its members, we applied the same Biolog EcoPlates test to the individual isolates. We found that substrate utilization differences between communities could be attributed to the presence or absence of particular strains. Out of the 20 strains measured using Biolog EcoPlates, only one showed high metabolic activity when grown with salicylic acid (strain M09D5GC17 corresponding to ASV863 Burkholderia), and it was a strain present only in M09. The M09 community was the only one where growth on salicylic acid was observed. Conversely, the strain with the most growth on itaconic acid (strain M10D5GC19 corresponding to ASV140 Achromobacter) was present at the final timepoint in all the microcosms we cultured from except for M09, and the M09 community was the only one not able to grow on itaconic acid (Supplementary Fig. 8). Individual strains can thus drive functional differences among microcosms, and we were able to culture and analyze a set of these strains—connecting genotypes with their functional phenotypes.
Source: Ecology - nature.com
