Response of the rare biosphere to environmental stressors in a highly diverse ecosystem (Zodletone spring, OK, USA)

Sampling

Sediment samples were collected from Zodletone spring, a sulfide and sulfur-rich spring in Southwestern Oklahoma (35°0′9″N, 98°41′17″W). The spring geochemistry has been previously described in detail (Bühring et al., 2011; Senko et al., 2004). The source of the spring is a contained area (1 m²) with anaerobic, biomass-laden, and sulfide-rich black viscous sediments, and anoxic, sulfide-rich (8.4 mM) 40-cm water column. Spring source water was also sampled for use in enrichments as detailed below. Samples were stored on ice until returning to the lab, where they were used for the experimental procedures described below within 24 h of sampling.

Enrichments preparation

The overall experimental setup is shown in Fig. 1.

Enrichments were prepared in Balch tubes under anaerobic conditions in an anaerobic chamber (Coy Laboratories, Grass Lake, Michigan, USA) by adding 1 g sediment to 4 ml autoclaved anoxic spring water. For this purpose, spring source water was autoclaved, cooled under a stream of N₂ to maintain anaerobic conditions, and used in all enrichments. To mimic environmental stress, either salinity or temperature was varied in the enrichments. Salinity was adjusted by adding NaCl to increase the concentration to 1%, 2%, 3%, 4%, or 10% above ambient concentration (measured at 0.9%) and these enrichments were incubated at room temperature. Temperature was varied by incubating enrichments at 28 °C, 30 °C, 32 °C, or 70 °C. For each condition, triplicate enrichments were incubated for 60 days in the dark. To account for any changes that could occur due to the incubation process itself, triplicate control enrichments were prepared by incubating 1g sediment in 4 ml autoclaved anoxic spring water with no further changes in salinity or temperature. This no-stressor enrichment control was subsequently used as the baseline for studying the effect of temperature and salinity change on the microbial community.

DNA extraction, PCR amplification and sequencing

Triplicates were pooled, centrifuged at 10,000 × g for 20 min and DNA was extracted from the obtained pellet using MoBio FastDNA Spin kit for Soil (MoBio, Carlsbad, California, USA) following the manufacturer’s instructions, and subsequently quantified using Qubit^® fluorometer (Life Technologies, Grand Island, New York, USA). Variable regions V1 and V2 of the 16S rRNA gene were then amplified using barcoded primers for multiplex sequencing using the FLX technology. The forward primer was constructed by adding 454 Roche FLX adaptor A (GCCTCCCTCGCGCCATCAG) to the 27F primer sequence (AGAGTTTGATCCTGGCTCAG) as previously described (Youssef, Couger & Elshahed, 2010). The forward primer also contained a unique barcode (octamer) sequence for multiplexing. The reverse primer was constructed by adding 454 Roche FLX adaptor B (GCCTTGCCAGCCCGCTCAGT) to the 338R primer sequence (GCTGCCTCCCGTAGGAGT). PCRs were conducted in 100 µl volume. The reaction contained 4 µl of the extracted DNA, 1 × PCR buffer (Promega, Madison, Wisconsin, USA), 2.5 mM MgSO₄, 0.2 mM dNTPs mixture, 0.5U of the GoTaq flexi DNA polymerase (Promega, Madison, Wisconsin, USA), and 10 µM each of the forward and the reverse primers. PCR was carried out according to the following protocol; initial denaturation at 95 °C for 5 min, followed by 35 cycles of denaturation at 95 °C for 45 s, annealing at 52 °C for 45 s, and elongation at 72 °C for 30 s. A final elongation step at 72 °C for 5 min was included. All PCR reactions were run in at least triplicates, and the resulting products of the expected size were combined and purified using QIAquick PCR cleanup kit (Qiagen Corp., Valencia, California, USA). Purified combined PCR products (11–15 µg, equivalent to ∼1 µg of each enrichment condition) were sequenced using FLX technology at the Environmental Genomics Core facility at the University of South Carolina.

Sequence quality filtering, OTU identification, and phylogenetic assignments

Sequence quality control was handled in mothur (Schloss et al., 2009) as described previously (Youssef, Couger & Elshahed, 2010). Briefly, sequences with an average quality score below 25, sequences that did not have the exact primer sequence, sequences that contained an ambiguous base (N), sequences having a homopolymer stretch longer than 8 bases, and sequences shorter than 80 bp were considered of poor quality and removed from the data set. High-quality reads from each treatment condition were aligned against the Greengenes alignment database using a Needleman-Wunsch pairwise alignment algorithm. Filtered alignments were used to generate an uncorrected pairwise distance matrix, followed by binning the sequences into operational taxonomic units (OTUs) at 3, 6, 8, 10, and 15% cutoffs. For phylogenetic placement, representative OTUs defined at the 3% cutoff (OTU_0.03) were classified with Greengenes taxonomy scheme using the PyNAST pipeline (Caporaso et al., 2010). Phylum level affiliation of sequences were determined according to the classifier output, and sequences with less than 85% similarity to their closest relative in Greengenes database were considered unclassified.

Defining the rare biosphere

The cutoff for defining rare biosphere is arbitrary and the methods used include relative abundance cutoffs, as well as frequencies of occurrence in a dataset (Youssef, Couger & Elshahed, 2010). For the sake of this study, we used a more relaxed definition for the rare biosphere, where rare members of the microbial community were identified using an empirical cutoff of n ≤ 10.

Examining microbial community response to environmental stressors

Our analysis had two main goals: To identify differences in the overall microbial community structure and diversity pre and post-enrichment, and to examine and quantify the contribution of various members of the rare biosphere in response to environmental stressors.

For the first goal (i.e., identification of differences in the overall microbial community structure and diversity due to environmental stress), microbial diversity was quantified using various diversity indices e.g., Shannon-Weiner, and Simpson diversity indices, Chao, and ACE estimators of species richness across enrichments. As well, beta diversity across samples was examined using rarefaction curve ranking since this method is not sensitive to sample size as described before (Youssef & Elshahed, 2009). Finally, variations in community structure due to stressors were calculated using Morisita-Horn index. This index was chosen due to its insensitivity to sample size variations (Anderson & Millar, 2004). Morisita-Horn index values obtained were subsequently employed to construct non-metric multidimensional scaling (NMDS) plots for community comparisons using the command nmds in mothur. The proportion of variance (r²) among communities was estimated from the NMDS plots by first calculating the Euclidean distance between all pairs of data points using the equation; $d=\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2}$, where d is the Euclidean distance between 2 points of coordinates (x₁, y₁) and (x₂, y₂) in ordination space. The Euclidean distance was then regressed on the beta-diversity index to estimate r². In order to study the effect of the stressors on the overall phylum-level community composition, likelihood-ratio-Chi-squared test was used to examine the significant difference between the relative abundances of phyla in the no-stressor enrichment control vs. the various enrichments.

Our strategy to achieve the second goal (i.e., understanding the contribution of the rare biosphere to environmental stressors) is based on the identification of the abundant members of the community post-stressor application in different enrichments, and tracing their origins to the various fractions of the community in the no-stressor control incubation. This allows for quantifying the contribution of the rare biosphere in the no-stressor control incubation to the microbial community that developed post application of stressors. The promotion of specific members of the rare biosphere in the control incubation to become abundant members in the post-stressor enrichments communities suggests a role for such rare biosphere members in responding to environmental stressors. Also, the relative abundances of such promoted members in the post-stressor enrichment community could be regarded as an indicator of the magnitude of importance of this promotion process.

Details of the analysis conducted are shown in Fig. 1. First, each sequence affiliated with the no-stressor control incubation was classified into either “abundant,” “unique rare,” or “non-unique rare” classes using the classification criteria detailed before (Youssef, Couger & Elshahed, 2010). Briefly, “abundant” refers to all sequences binned into OTUs with >10 representatives, “unique rare” refers to all sequences binned into OTUs with ≤10 representatives and that are phylogenetically distinct from more abundant members of the community (>85% sequence similarity), while “non-unique rare” refers to all sequences binned into OTUs with ≤10 representatives and that are phylogenetically similar to more abundant members of the community (Youssef, Couger & Elshahed, 2010). We then queried all abundant (n > 10) OTU_0.03 representatives from individual post-stressor treatments against all sequences recovered from the no-stressor control experiment (14,071 sequences) using local Blastn. Identification of the best “hit” for each abundant OTU_0.03 representative in the post-stressor incubations allowed us to trace its origin in the no-stressor control incubation into either “originating from abundant,” “originating from unique rare,” or “originating from non-unique rare” sequence. Abundant post-stressor OTUs originating from an “abundant” control sequence represent the “no change” effect, where an abundant member remained abundant post-stressor application. Abundant post-stressor OTUs originating from a “rare” control sequence represent the “promotion” effect, where a rare member in the no-stressor incubation control was promoted to a more abundant OTU post-stressor. This latter group could be further subdivided into two distinct categories: Abundant post-stressor OTUs “originating from non-unique rare” control sequence, and abundant post-stressor OTUs “originating from unique rare” control sequence. Based on these results, we calculated the percentage of post-stressor abundant OTUs that were “originating from abundant,” “originating from unique rare,” or “originating from non-unique rare” sequences, and correlated these values to the severity of enrichment (salinity or temperature) using Pearson correlations.

Controlling the false rate of discovery using Benjamini–Hochberg adjustment

P-values were calculated for all Pearson correlation coefficients. Due to the large number of correlations conducted in this study, and to avoid any spurious correlations that might arise by chance (rejecting the null hypothesis due to a high p-value, when in fact the null hypothesis is true), we used the Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995) to adjust for the false discovery rate. The command p.adjust in R was used.

Overall sequencing results

A total of 74,017 high quality sequences were obtained in this study. Estimates of diversity, species richness, evenness, and beta diversity at the species level (OTU_0.03) in various enrichments are shown in Table 1. As expected, multiple sequences affiliated with various sulfur-metabolizing lineages were identified within the datasets, including members of the sulfate- and sulfur-reducing δ-Proteobacteria, chemolithotrophic sulfur-oxidizing β- and γ-Proteobacteria, as well as purple (α- and γ-Proteobacteria) and green (Chloroflexi and Chlorobia) sulfide- and sulfur-oxidizing anoxygenic phototrophs (Table S1).

Effect of environmental stressors on overall microbial community

In general, the increase in temperature negatively affected the community diversity (measured as rarefaction curve rank, Fig. S1A and Fig. S1B) (Pearson correlation coefficient = − 0.67, P = 0.17, Fig. S2A), with the magnitude of diversity loss directly correlated to the prevailing temperature. On the other hand, all salinity treatments resulted in an increase in community diversity compared to the no-treatment control, as evident by diversity ranking estimates obtained (Table 1). However, that increase was not directly correlated to increments in salinity (Pearson correlation coefficient = 0.06, Fig. S2B).

At the phylum level, 50 distinct bacterial phyla and candidate phyla were identified in this study. The highest phylum level diversity was observed in sediments obtained directly from the spring prior to enrichment (50 phyla), followed by the no-treatment control enrichment (35 phyla, Table S1 and Fig. 2). Similar to diversity statistics at the putative species level, the number of phyla identified was negatively correlated to the enrichment temperature (Pearson correlation coefficient = − 0.66, P = 0.17, Fig. S1A and Fig. S1B), while salinity had no clear effect on phylum level diversity (Pearson correlation coefficient = − 0.09) (Table 1).

Likelihood-ratio-Chi-squared test was used to examine the significant difference between the relative abundances of phyla in the no-treatment enrichment control vs. the various enrichments. This method failed to identify any significant difference (likelihood ratio χ2 = 99.4, p = 1). Nevertheless, Pearson correlations between specific phyla relative abundance in the various enrichments, and the enrichment condition (salinity or temperature) identified the following patterns (Fig. 2 and Table S1). The increase in the enrichment temperature was positively correlated to the percentage abundance of Firmicutes, Aminicenantes (candidate division OP8), Parcubacteria (candidate division OD1), and Thermotogae (Pearson correlation coefficients = 0.98 (P = 0.0005), 0.9(P = 0.014), 0.98 (P = 0.0003), and 0.75 (P = 0.1), respectively, Fig. S3A, Fig. S3B, Fig. S3C and Fig. S3D), and negatively correlated to the percentage abundances of Bacteroidetes, Marinimicrobia (candidate division SAR406), Verrucomicrobia and Latescibacteria (candidate division WS3) (Pearson correlation coefficients = − 0.9 (P = 0.016), −0.75 (P = 0.1), −0.82 (P = 0.06), and −0.75 (P = 0.1), respectively, Fig. S3A, Fig. S3B, Fig. S3C and Fig. S3D). The increase in enrichment salt concentration was positively correlated to the percentage abundance of Firmicutes, and Thermotogae (Pearson correlation coefficients = 0.53 (P = 0.258), and 0.52 (P = 0.269), respectively, Fig. S3A, Fig. S3B, Fig. S3C and Fig. S3D), and negatively correlated to the percentage abundances of Verrucomicrobia and Latescibacteria (candidate division WS3) (Pearson correlation coefficients = − 0.6 (P = 0.18), and −0.73 (P = 0.08), respectively, Fig. S3A, Fig. S3B, Fig. S3C and Fig. S3D). All P-values shown in boldface were significant following Benjamini–Hochberg adjustment (Benjamini & Hochberg, 1995).

Shifts in dominant microbial populations post-stressor application

To examine the impact of stressors on the microbial community, we examined the occurrence and magnitude of shifts in the structure of abundant community members post enrichments, compared to the no-treatment control. Multiple lines of evidence suggest a shift in the abundant community post-stressor application: The proportion of sequences belonging to abundant OTUs (n > 10) decreased post enrichment (Table 1). In addition, high levels of beta diversity were observed between the no-treatment incubation and all salinity and temperature enrichments. The abundant no-treatment community showed Morisita Horn indices of 0.63 ± 0.07, and 0.57 ± 0.3 for all possible pair-wise comparisons with the salinity and temperature post-stressor abundant communities, respectively. Indeed, non-metric multidimensional scaling plot using Morisita-Horn indices clearly shows a shift in the abundant community structure following stress application (Fig. 3). Finally, phylogenetic affiliations of abundant members of the community following application of stressors showed marked differences from those in the no treatment control (Fig. 4).

Origins of the abundant community members following enrichment support a role for the rare biosphere in responding to environmental stressors

We hypothesized that the differences observed in dominant members of the community following application of stressors are due to recruitment of organisms from the rare biosphere to become part of the abundant members of the communities. To this end, we sought to identify the origin of all members of the abundant community post-stressor application (Fig. 1) and determine their origin (abundant, rare unique, and rare non-unique) in the no-treatment control community. Our analysis identified three distinct patterns (Table 2): First, post-enrichment abundant OTUs that were similar to “abundant” sequences in the no-treatment control. This group constituted 36.7–91.1% of the community in various enrichments. As expected from the community structure analysis, the highest percentages were encountered with the 28 °C and 30 °C incubations, since the abundant community from these temperature incubations was very similar to the no-treatment control (Fig. 3). The percentage of this group decreased as the enrichment salinity, and temperature increased (Pearson correlation coefficient = − 0.6, and −0.86 (P = 0.07) for salinity, and temperature, respectively, Fig. S4A and Fig. S4B). The second group constitutes abundant OTUs in the enrichments (8.9–63.3%) that were recruited from the no-treatment control rare biosphere, providing direct evidence that the rare biosphere could act as a backup system to respond to environmental stressors. Within this group, we differentiate between two distinct fractions: (A) Post-enrichment abundant OTUs that were promoted from “rare non-unique” sequences in the no-treatment control, i.e., rare members of the original community phylogenetically similar to more abundant members of the community. Of the total number of abundant OTUs, this fraction constituted 2.6–10.2% in various conditions. This percentage decreased as the enrichment salinity, and temperature increased (Pearson correlation coefficient = − 0.86 for both salinity (P = 0.035), and temperature (P = 0.076), Fig. S4A and Fig. S4B). And (B) Post-enrichment abundant OTUs that were promoted from “rare unique” sequences in the no-treatment control, i.e., rare members of the original community phylogenetically distinct from more abundant members of the community. Of the total number of abundant OTUs, this fraction constituted 1.8–60.7% of abundant OTUs in various conditions. This percentage increased as the enrichment salinity, and temperature increased (Pearson correlation coefficient = 0.85 (P = 0.039), and 0.86 (P = 0.072) for salinity, and temperature, respectively Fig. S4A and Fig. S4B). Therefore, we argue that, while both factions of the rare biosphere are important in responding to changes in environmental conditions and both act to recruit members to the abundant community, the magnitude of contribution of the “non-unique” rare biosphere to the promotion process was less significant. On the other hand, the “unique” rare biosphere seemed to contribute to the promotion process both when the magnitude of environmental stressors applied was slight, e.g., similar to what would happen during diurnal variation in salinity and temperature, as well as severe, e.g., similar to what would be encountered during seasonal variation in temperature or salinity, or following a drastic change in environmental condition, e.g., drought or fire. This latter effect is so evident in the 10% salt, and the 70 °C incubation, where 55%, and 61%, respectively of the abundant community was recruited from the “unique” rare biosphere (Table 2).