Decreases in average bacterial community rRNA operon copy number during succession.

Trait-based studies can help clarify the mechanisms driving patterns of microbial community assembly and coexistence. Here, we use a trait-based approach to explore the importance of rRNA operon copy number in microbial succession, building on prior evidence that organisms with higher copy numbers respond more rapidly to nutrient inputs. We set flasks of heterotrophic media into the environment and examined bacterial community assembly at seven time points. Communities were arrayed along a geographic gradient to introduce stochasticity via dispersal processes and were analyzed using 16 S rRNA gene pyrosequencing, and rRNA operon copy number was modeled using ancestral trait reconstruction. We found that taxonomic composition was similar between communities at the beginning of the experiment and then diverged through time; as well, phylogenetic clustering within communities decreased over time. The average rRNA operon copy number decreased over the experiment, and variance in rRNA operon copy number was lowest both early and late in succession. We then analyzed bacterial community data from other soil and sediment primary and secondary successional sequences from three markedly different ecosystem types. Our results demonstrate that decreases in average copy number are a consistent feature of communities across various drivers of ecological succession. Importantly, our work supports the scaling of the copy number trait over multiple levels of biological organization, ranging from cells to populations and communities, with implications for both microbial ecology and evolution.

Introduction

A broader understanding of the forces guiding microbial community assembly is important for predictions about the effects of environmental change on microbially mediated ecosystem processes (Hanson ; Nemergut ). Several studies have documented the roles of species sorting (Langenheder and Székely, 2011; Evans and Wallenstein, 2014), competition (Goberna ; Koeppel and Wu, 2014), priority effects (Fukami ; Andersson ), dispersal (Langenheder ; Lindström and Langenheder, 2012) and ecological drift (Ofiteru ; Stegen , 2013) in microbial community assembly. Further, several studies have shown that the relative importance of these different assembly processes changes with nutrient status (Langenheder ; Knelman ), disturbance (Ferrenberg ) and successional processes (Zhou ). Also, the relative importance of assembly processes may vary between different fractions of the community based on ecological strategies and/or relative abundance (Langenheder and Székely, 2011). However, we lack a comprehensive, theoretical understanding of why different microbial assembly processes vary in time and space. One key step toward developing a more mechanistic appreciation of microbial community assembly is to understand the role of functional traits in community assembly (Violle ; Green ).

Trait-based approaches may greatly inform microbial ecology given how they have advanced our understanding of plant community assembly by providing hypotheses for mechanisms of assembly that apply across ecosystems (Garnier ; Ackerly and Cornwell, 2007; Kraft ; Webb ). For example, Shipley used a statistical approach to identify eight plant traits important for plant community assembly processes, which explained 94% of variation in plant community composition across a successional sere. Likewise, microorganisms exhibit variation in a variety of key traits that may scale from individuals to communities (Fierer ; Goberna ), some of which may correlate with distribution patterns (Fierer ; Lennon ). Indeed, trait-based approaches may be especially important to apply to microbial systems to make sense of assembly mechanisms in light of high levels of biodiversity, horizontal gene transfer and redundancy, which may otherwise complicate taxonomic and/or phylogenetic approaches (Krause ).

One trait that may be important in heterotrophic bacterial community assembly is rRNA operon copy number (Klappenbach ; Shrestha ; Yano ). rRNA operon number varies over an order of magnitude between different bacterial species (Klappenbach ) and has been correlated with rates and efficiencies of a variety of cellular processes including growth rate, lag time, sporulation efficiency, motility and competence (Condon ; Klappenbach ; Stevenson and Schmidt, 2004; Yano ; Lauro ). In general, rRNA operon copy number correlates positively with attributes that could be classified as ‘r-selected' (sensu Reznick ), and genetic studies have revealed the importance of rRNA copy number on microbial physiology under conditions of increasing nutrients (Stevenson and Schmidt, 2004; Lauro ).

As variation in nutrient supply is often a characteristic of successional sequences, and the importance of life history characteristics have long been described in plant community successional theory (Grime, 1977; Matthews, 1992; Chapin ; Walker and del Moral, 2003), here we ask if rRNA operon copy number varies in a predictable pattern in microbial community assembly across succession. We placed microcosms of sterile laboratory media out into the environment and examined bacterial community composition and rRNA operon copy number as a community aggregated trait (Garnier ; Fierer ) at seven time points during the first 99 h of community development. We hypothesized that the average rRNA operon copy number would be high in earlier communities and would decrease throughout succession. We further tested this hypothesis by analyzing bacterial community data from other soil and sediment successional sequences from markedly different environments. Our results suggest that decreases in average rRNA operon copy number are a consistent feature of communities across multiple ecosystem types.

Materials and methods

We arrayed microcosms of heterotrophic media along a 15-km transect beneath the understory of a tropical coastal forest on the Yucatan Peninsula, Mexico (climate data available here: http://www.worldclim.org/). We added 100 ml of autoclaved, filter-sterilized nutrient broth (Difco Laboratories, Detroit, MI, USA) containing 10 mg ml−1 cycloheximide (Sigma-Aldrich, St Louis, MO, USA) as a eukaryotic inhibitor to sterile Mason jars covered with a fine (1 mm mesh) screen. To investigate the role of stochasticity via dispersal, duplicate microcosms (Mason jars) were arrayed along the transect in locations that were at least 1 km apart (Supplementary Figure 1A). At 3, 19, 23.5, 27.5, 43.5, 51.5, 67 and 99 h, microcosms were sampled (4 ml) after stirring the media with a sterile pipet.

Half of each sample was used to measure optical density at 595 nm while the other half was preserved at −20 °C for DNA extraction using the UltraClean Microbial DNA (MO BIO Laboratories, Carlsbad, CA, USA) kit. We amplified 16 S rRNA genes for pyrosequencing using the primers and conditions described in Fierer and Knelman ; none of the samples taken at 3 h yielded PCR products that were visible on an agarose gel and these samples were thus discarded from the study. As well, the sterility of microcosm 9B (all time points) and of sample 8B at time point 5 were compromised and these samples were discarded from future analyses. Pyrosequencing was performed by Engencore (University of South Carolina).

Pyrosequencing data were screened using the QIIME toolkit (Caporaso ) with the following parameters: quality score >25, sequence length >200 and <400, maximum homopolymer of 6, 0 maximum ambiguous bases, and 0 mismatched bases in the primer. First, operational taxonomic units (OTUs) were picked at the 97% identity level using UCLUST (Edgar, 2010) in QIIME with the default settings and data were denoised using Denoiser (Reeder and Knight, 2010). OTUs were randomly subsampled in QIIME such that each library contained 608 sequences, the fewest sequences in a single sample. In addition, OTUs were selected using the closed reference command in QIIME against the gg_13_5otus representative set at the 97% identity level (DeSantis ) for phylogenetic identification and copy number estimation. We then used the web-based implementation of PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States, Langille ) to estimate rRNA operon copy number. OTU data were then normalized (standardized) for copy number by dividing by copy number. For each sample, we calculated the community aggregated trait value (weighted mean) by taking the product of the estimated operon copy number and the relative abundance for each OTU, and summing this value across all OTUs in a sample. Operon copy number variance was similarly calculated as a weighted variance.

We found a strong correlation between standardized and nonstandardized data as well as between UniFrac, a phylogeny-based beta diversity metric (Lozupone ), and Bray–Curtis dissimilarity, an OTU-based beta diversity metric (Mantel tests on pairwise dissimilarity matrices all had correlations >0.8 and were all significant at P<0.001). We chose to perform all subsequent analyses on the Bray–Curtis dissimilarities of standardized data because of our focus on individual taxa. We also supplemented these analyses with an analysis of phylogenetic clustering (Webb ; see below). We used the vegan package (Oksanen ) in R (R Development Core Team, 2013) to generate an ordination plot for visualization of community shifts. The best solution, as determined by the lowest stress and visual examination of stress plots, was found by performing NMDS on square-root transformed Bray–Curtis dissimilarity.

To examine changes in taxonomic community structure based on the Bray–Curtis dissimilarity matrix, we performed an analysis of multivariate homogeneity of group dispersions (Anderson, 2006) to test for changes in community composition over time and used permutational multivariate analysis of variance (PERMANOVA, Anderson, 2001) to test for changes in beta diversity over time. These analyses were conducted in QIIME using the PERMANOVA function and PERMDISP function (which in turn uses the betadisper function of the vegan package in R). We calculated the nearest taxon and net relatedness indices (NTI and NRI, respectively (Webb, 2000)) for each time point using the ses.mpd and the ses.mntd functions in the picante package in R and multiplying values by −1 (Kembel ). For this analysis, we used the filter_tree command in QIIME to trim the gg_13_5 Greengenes tree to contain only OTUs that were observed in this study. We used ANOVA and Kruskal–Wallis post hoc tests (pgirmess package) in R to determine differences in community aggregated trait values for rRNA copy number, NTI and NRI across time. Finally, we generated frequency-abundance plots for successional time points 23.5 h and onwards (19 h showed nondistinguishable patterns from the 23.5 h time point) and implemented the Sloan neutral community model to determine if frequency-abundance patterns were indicative of neutral assembly.

We also examined three additional data sets relating to primary and secondary successional sequences: a salt marsh sediment (Dini-Andreote ), a recently deglaciated soil from the Toklat glacier in Alaska, and a burned forest soil from the Four Mile Canyon Fire outside of Boulder, CO (Ferrenberg ). Soils from the Toklat glacier were sampled in July of 2014. Samples were collected from frost-heave mounds along the glacial chronosequence by collecting the top 2 cm of material. Samples were frozen upon removal from the field and shipped to the University of Colorado at Boulder and kept at −80 °C until processing. DNA was extracted from each sample using the MO BIO PowerSoil DNA extraction kit (MO BIO Laboratories). The 16 S rRNA gene was amplified using the primer set 515 F/806 R (Caporaso ) and sequenced on the Illumina MiSeq platform (Illumina Inc., San Diego, CA, USA) using paired-end 2x150 bp chemistry. A 30% phiX spike was added to the run to compensate for otherwise limited amplicon variability (Caporaso ). Paired-end reads were demultiplexed and assembled using QIIME (Caporaso ).

Soils from the Four Mile Canyon burn site were sampled and sequenced as described in Ferrenberg : 10 soil samples from 29 months after the burn were compared with samples from 4 months after the burn. These data sets were analyzed for rRNA operon copy number according to the workflow previously described. Statistical differences between successional categories were examined using nonparametric Kruskal–Wallis ANOVA (pgirmess package). The Toklat glacier forefield data set was divided into the unvegetated and vegetated chronosequence for this categorical analysis and a linear regression model on the natural log transformed response variable was also performed.

All sequence data have been archived and are publically available at th; doi:10.6084/m9.figshare.1556152.

Results and Discussion

Most bacterial communities grew to stationary phase (Supplementary Figure 2) over the course of the experiment. Early communities were dominated by organisms related to Bacillus (Supplementary Figure 3). Later successional communities showed increases in the relative abundance of organisms related to the Enterobacteriales and Pseudomonadales and significantly changed in taxonomic composition over time (PERMANOVA, P=0.001). Beta diversity in early successional communities was low: NMDS plots showed that communities clustered together and were characterized by low levels of taxonomic dispersion (Figures 1a and b). After 43 h, beta diversity was significantly higher than in the early successional communities (PERMDISP contrasts, P=0.001) and remained high until the end of the experiment (Figure 1b). These results contrast with those of Bell (2010) who found that air-seeded heterotrophic microbial communities converged through time. This may be due to the richer media used in the present study, which could change the relative roles of stochasticity and determinism in community assembly (Keddy ; Chase, 2003, 2010; Kardol ), and/or due to the differences in the time scales of the studies (99 h vs 28 days).

Figure 1

(a) NMDS plot showing an increase in beta diversity over successional time. Different colors represent different microcosms from distinct geographic distances (see Supplementary Figure 1); hatched lines of the same color represent replicate microcosms. Points are joined to show successional trajectories over the course of the experiment; terminal arrows show the community composition in the microcosms at the end of the experiment. The scale represents relative Bray–Curtis dissimilarity values. (b) Boxplot showing an increase in beta diversity through successional time. The plot displays Bray–Curtis dissimilarity of each replicate to the centroid (a measure of beta diversity) for each time point. PERMDISP analysis demonstrates significant increases in beta diversity (measured as multivariate dispersion), with samples after 43 h being significantly more dispersed than earlier time points (PERMDISP, P=0.001).

Mean community aggregated rRNA operon copy number started out high across samples, at more than nine copies (Figure 2). By the end of the experiment, the mean community aggregated copy number had declined to an average of 3.6 copies, significantly lower than the first three successional time points (19–27.5 h; Kruskal–Wallis contrasts P<0.05; Supplementary Table 1). Variance in copy number was lowest at the beginning and the end of the experiment and was highest between 23.5 and 43.5 h. Analysis of copy number trends within individual flasks showed that this variance appeared to be largely driven by differences in the timing of a decrease in mean aggregated copy number, and while most communities showed consistent decreases in copy number through time, this variance was the result of asynchronous decreases rather than stochastic fluctuations. Average copy number therefore appears to be a community-level trait (Fierer ) that serves as an indicator for ecological history. For example, recently DeAngelis demonstrated that soil warming significantly decreased average rRNA copy number, possibly reflecting shifts in microbial substrate availability and thus growth dynamics.

Figure 2

Scatterplot showing a decrease in the mean weighted rRNA copy number over successional time as well as lower variance in copy number in both early and late communities. Copy number was estimated using PICRUSt (Langille ) and weighted values were obtained by multiplying copy numbers by relative abundance for each OTU and taking the summation of these values for each community. Bars represent s.d.

Frequency-abundance plots (Sloan ) revealed that the most locally abundant and widely distributed organisms across samples had high average operon copy number at the beginning of the experiment but that there was a marked transition toward organisms with lower copy numbers that are less widely distributed through successional time (Figure 3). Thus, together, our results support a simultaneous divergence in community composition and a convergence in mean aggregated rRNA operon copy during microbial succession.

Figure 3

Frequency-abundance plots over successional time that show an increase in the relative importance of taxa with lower rRNA copy numbers. The plots at both 19 and 23 h were similar so only the 23-h plot is shown. Each point represents a unique OTU and we have plotted its average relative abundance within a community against its frequency (distribution across samples) among communities. Points are colored by the rRNA copy number as estimated using PICRUSt (Langille ).

Similar patterns of trait convergence have been observed for other microbial communities and for plant communities (Fukami ; Burke ), and the patterns that we observed here allow us to generate more mechanistic hypotheses about the importance of different assembly mechanisms in this experimental system. In other studies, the importance of stochastic assembly processes has been proposed in early microbial succession based on deviation from null models (Ferrenberg ; Zhou ) and based on the analysis of more controlled microcosm manipulations (Langenheder and Székely, 2011). Likewise, here our data show that early community composition conforms to frequency-abundance relationships that would be expected under neutral processes (Figure 3), with data from the first two time points fitting the Sloan neutral community model (NTm= 5.3 and 5.9 for 19 and 27.5 h, respectively), while data from later time points do not.

Yet, this study is important in demonstrating that while the observed taxonomic distribution relationships may reflect stochasticity in the assembly of organisms present within the early successional communities, there is low variance in the mean aggregated rRNA operon copy number in early successional communities (Figure 2). Thus, organisms with high operon copy number may have a selective advantage at the first time points, likely because of their ability to respond rapidly to nutrient inputs (Klappenbach ; Yano ) and/or because of their ability to form spores (Yano ). At the same time, these organisms may not experience significant competition pressure because of combined high nutrient/low biomass levels in early succession. Indeed, calculations of the NRI and NTI, respectively (Webb, 2000; Webb ) reveal phylogenetic clustering of early successional communities and a significant shift to less-clustered communities over time, (Figure 4; ANOVA and Kruskal–Wallis contrasts, P<0.05 for NTI and NRI; Supplementary Table 1). In addition, we know of no outdoor air-sample communities that have been analyzed using molecular gene surveys, which contain Bacilli at such high relative abundance (~90%) as we observed here in early successional communities (Brodie ; Bowers , 2011). Thus, in addition to the support for stochastic processes in early assembly, these data also simultaneously support a role for environmental selection at a potential trait level—rRNA operon copy number—in early succession. The forces acting on rRNA copy number maintenance are more complicated than a simple efficiency-growth rate trade-off (Stevenson and Schmidt, 2004); however, more research will be needed to uncover the relevant evolutionary and ecological pressures on these operons.

Figure 4

(a) Net Relatedness Index and (b) Nearest Taxon Index over successional time. Both indices show a significant decrease in phylogenetic clustering as communities assemble. For these analyses we used the Greengenes gg_13_5 tree that was trimmed in QIIME to contain only OTUs observed in this study. Bars represent s.d.

Our data also reveal a significant increase in OTU-related beta diversity—pairwise dissimilarity—over time (Figures 1a and b), which, combined with the relatively controlled environmental conditions, may lead to the hypothesis that later successional communities are more stochastically assembled. Although this hypothesis may arise at the taxonomic level, the trait-based data suggest that later successional communities are dominated by organisms with relatively low rRNA operon copy numbers. Thus, this low variance in mean aggregated copy number suggests that this trait may engender a selective advantage in the context of later successional communities, pointing to rRNA copy number as a fundamental trait in microbial life history strategies.

Later successional communities were also significantly overdispersed (less phylogenetically clustered) than earlier communities (Figure 4). Overdispersion patterns in phylogeny have been interpreted as support for the role of competition in community assembly (Webb, 2000; Kraft ; Violle ). Stabilizing niche differences require density-dependent growth and thus may become increasingly important as community size (for example, number of individuals) increases (Chesson, 2000; HilleRisLambers ). An alternative (but nonmutually exclusive) hypothesis is that successional time allowed for greater immigration of different species from the regional species pool, which resulted in the development of more phylogenetically diverse communities. However, while both dispersal and competition likely have roles in community assembly in this system, the strong trait convergence that we observed in the later successional communities (Figure 2) clearly supports the supposition that environmental selection is key in community assembly processes in these microcosms.

An important alternative hypothesis that could explain the patterns observed in this particular microcosm study is that the link between average copy number and successional stage may be due to oxygen availability, and because of the low diversity of microcosm communities, any trends in average copy number may be coincidental. Indeed, four dominant taxa: Bacillus (primarily aerobic), Enterobacteriales (facultative), Pseudomonadales (facultative) and Clostridia (anaerobic), all generally shift in dominance with advancing successional stage (Figure 3). To address this hypothesis, we examined trends in copy number in a suite of diverse successional systems that feature much higher levels of microbial diversity and vastly different environmental conditions: a salt marsh sediment primary-successional community (Dini-Andreote ), a recently deglaciated Alaskan soil, and a secondary-successional burned soil from Colorado (Ferrenberg ). In each case, we observed a significant decrease in average copy number with successional time (Figure 5, Kruskal–Wallis contrasts P<0.05). The glacier chronosequence, sampled in a continuous manner, also revealed a significant relationship between rRNA copy number and successional time (linear regression, multiple R2=0.4949, P<0.0001). Although a full metaanalysis is beyond the scope of this paper, it is clear that decreases in average rRNA operon copy number are a consistent feature of a diverse set of microbial successional communities.

Figure 5

Boxplots and scatterplot showing a decrease in the mean weighted rRNA copy number over successional time for the three additional communities examined: a burned forest soil (a) (Ferrenberg ); a salt marsh sediment chronosequence (b) (Dini-Andreote ); and a glacier forefield chronosequence (c), which included linear regression analysis (multiple R2=0.4949, P<0.0001). Copy number was estimated using PICRUSt (Langille ) and weighted values were obtained by multiplying rRNA copy numbers by relative abundance for a given OTU, and summing weighted averages across a given sample.

In the microcosm study, we observed an increase in beta diversity as well as a decrease in phylogenetic clustering and weighted mean rRNA operon copy number during microbial community succession. The low variance in copy number both early and late in assembly supports the role of environmental selection in this microbial trait during succession, leading us to hypothesize that a larger metacommunity of low copy number organisms allowed for the divergence of community composition and simultaneous trait conversion. Thus, while a taxonomy or phylogeny-only approach would have led to hypotheses primarily focused on the role of stochasticity in both early and late community development, patterns of trait dispersion support the importance of environmental selection over the course of succession. Interestingly, the three other systems examined did not reveal consistent successional trends in taxonomic dispersion, with the salt marsh study showing a decrease in dispersion and the soil sequences showing no change in beta diversity over the time period examined. The fact that all four chronosequences displayed consistent decreases in copy number despite stark differences in both environmental conditions and successional trajectories suggests copy number may be a core trait for bacteria underlying basic life history strategies that is generalizable across environments. We expect that other traits may similarly shift over the course of succession in a directional and predictive manner, but because such traits likely relate to more specific conditions of the environment, rRNA copy number is powerful in its apparent importance in assembly across systems.

Trait-based studies are also essential for understanding the implications of assembly processes for ecosystem function. The relationship between the traits that structure communities (that is, response traits) and those traits that affect the environment (that is, effect traits) is key to understanding the implications of community assembly for ecosystem processes (Allison, 2012; Ferrenberg ; Knelman and Nemergut, 2014; Nemergut ). For example, Nemergut hypothesized that stochastic assembly processes could decouple community composition from environmental conditions, which could lead to unpredictable relationships between the environment and ecosystem processes as well as between biodiversity and ecosystem function (Knelman and Nemergut, 2014). Thus, a more trait-centric view on this and other assembly-function hypotheses is essential to predict and understand the impacts of environmental change on ecosystem processes (Lavorel and Grigulis, 2012; Evans and Wallenstein, 2014). Overall, our work supports the scaling of copy number from cells to populations and communities, suggesting that it may be a core trait that can lend insight into the ecological and evolutionary mechanisms driving patterns in microbial biodiversity across multiple levels of biological organization.