Life history strategies among soil bacteria—dichotomy for few, continuum for many

Data were analyzed from samples collected, processed, and published previously [21, 25, 29] and have been summarized here. The present analysis, which consisted of sequence data processing, the calculation of taxon-specific isotopic signatures, and subsequent analyses, reflects original work.

Sample collection and isotope incubation

To generate experimental data, three replicate soil samples were collected from the top 10 cm of plant-free patches in four ecosystems along the C. Hart Merriam elevation gradient in Northern Arizona. From low to high elevation, these sites are located in the following environments: desert grassland (GL; 1760 m), piñon-pine juniper woodland (PJ; 2020 m), ponderosa pine forest (PP; 2344 m), and mixed conifer forest (MC; 2620 m). Soil samples were air-dried for 24 h at room temperature, homogenized, and passed through a 2 mm sieve before being stored at 4 °C for another 24 h. This produced three distinct but homogenous soil samples from each of the four ecosystems that were subject to experimental treatments. Three treatments were applied to bring soils to 70% water-holding capacity: water alone (control), water with glucose (C treatment; 1000 µg C g⁻¹ dry soil), or water with glucose and a nitrogen source (CN treatment; [NH₄]₂SO₄ at 100 µg N g⁻¹ dry soil). To track growth through isotope assimilation, both ¹⁸O-enriched water (97 atom %) and ¹³C-enriched glucose (99 atom %) were used. In all treatments isotopically heavy samples were paired with matching “light” samples that received water with a natural abundance isotope signatures. For ¹⁸O incubations, this design resulted in three soil samples per ecosystem per treatment (across four ecosystems and three treatments, n = 36) while ¹³C incubations were limited to only C and CN treatments (n = 24). Previous analyses suggest that three replicates is sufficient to detect growth of 10 atom % ¹⁸O in microbial DNA with a power of 0.6 and a growth of 5 atom % ¹⁸O with a power of 0.3 (12 and 6 atom % respectively for ¹³C) [30]. All soils were incubated in the dark for one week. Following incubation, soils were frozen at −80 °C for one week prior to DNA extraction.

Quantitative stable isotope probing

The procedure of qSIP (quantitative stable isotope probing) is described here but has been applied to these samples as previously published [17, 21, 25]. DNA extraction was performed on soils using a DNeasy PowerSoil HTP 96 Kit (MoBio Laboratories, Carlsbad, CA, USA) and following manufacturer’s protocol. Briefly, 0.25 g of soils from each sample were carefully added to deep, 96-well plates containing zirconium dioxide beads and a cell lysis solution with sodium dodecyl sulfate (SDS) and shaken for 20 min. Following cell lysis, supernatant was collected and centrifuged three times in fresh 96-well plates with reagents separating DNA from non-DNA organic and inorganic materials. Lastly, DNA samples were collected on silica filter plates, rinsed with ethanol and eluted into 100 µL of a 10 mM Tris buffer in clean 96-well plates. To quantify the degree of ¹⁸O or ¹³C isotope incorporation into bacterial DNA (excess atom fraction or EAF), the qSIP protocol [31] was used, though modified slightly as reported previously [21, 24, 32]. Briefly, microbial growth was quantified as the change in DNA buoyant density due to incorporation of the ¹⁸O or ¹³C isotopes through the method of density fractionation by adding 1 µg of DNA to 2.6 mL of saturated CsCl solution in combination with a gradient buffer (200 mM Tris, 200 mM KCL, 2 mM EDTA) in a 3.3 mL OptiSeal ultracentrifuge tube (Beckman Coulter, Fullerton, CA, USA). The solution was centrifuged to produce a gradient of increasingly labeled (heavier) DNA in an Optima Max bench top ultracentrifuge (Beckman Coulter, Brea, CA, USA) with a Beckman TLN-100 rotor (127,000 × g for 72 h) at 18 °C. Each post-incubation sample was thus converted from a continuous gradient into approximately 20 fractions (150 µL) using a modified fraction recovery system (Beckman Coulter). The density of each fraction was measured with a Reichart AR200 digital refractometer (Reichert Analytical Instruments, Depew, NY, USA). Fractions with densities between 1.640 and 1.735 g cm⁻³ were retained as densities outside this range generally did not contain DNA. In all retained fractions, DNA was cleaned and purified using isopropanol precipitation and the abundance of bacterial 16S rRNA gene copies was quantified with qPCR using primers specific to bacterial 16S rRNA genes (Eub 515F: AAT GAT ACG GCG ACC ACC GAG TGC CAG CMG CCG CGG TAA, 806R: CAA GCA GAA GAC GGC ATA CGA GGA CTA CVS GGG TAT CTA AT). Triplicate reactions were 8 µL consisting of 0.2 mM of each primer, 0.01 U µL⁻¹ Phusion HotStart II Polymerase (Thermo Fisher Scientific, Waltham, MA), 1× Phusion HF buffer (Thermo Fisher Scientific), 3.0 mM MgCl₂, 6% glycerol, and 200 µL of dNTPs. Reactions were performed on a CFX384 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) under the following cycling conditions: 95 °C at 1 min and 44 cycles at 95 °C (30 s), 64.5 °C (30 s), and 72 °C (1 min). Separate from qPCR, retained sample-fractions were subject to a similar amplification step of the 16S rRNA gene V4 region (515F: GTG YCA GCM GCC GCG GTA A, 806R: GGA CTA CNV GGG TWT CTA AT) in preparation for sequencing with the same reaction mix but differing cycle conditions – 95 °C for 2 min followed by 15 cycles at 95 °C (30 s), 55 °C (30 s), and 60 °C (4 min). The resulting 16S rRNA gene V4 amplicons were sequenced on a MiSeq sequencing platform (Illumina, Inc., San Diego, CA, USA). DNA sequence data and sample metadata have been deposited in the NCBI Sequence Read Archive under the project ID PRJNA521534.

Sequence processing and qSIP analysis

Independently from previous publications, we processed raw sequence data of forward and reverse reads (FASTQ) within the QIIME2 environment [33] (release 2018.6) and denoised sequences within QIIME2 using the DADA2 pipeline [34]. We clustered the remaining sequences into amplicon sequence variants (ASVs, at 100% sequence identity) against the SILVA 138 database [35] using a pre-trained open-reference Naïve Bayes feature classifier [36]. We removed samples with less than 3000 sequence reads, non-bacterial lineages, and global singletons and doubletons. We converted ASV sequencing abundances in each fraction to the number of 16S rRNA gene copies per gram dry soil based on qPCR abundances and the known amount of dry soil equivalent added to the initial extraction. This allowed us to express absolute population densities, rather than relative abundances. Across all replicates, we identified 114 543 unique bacterial ASVs.

We calculated the ¹⁸O and ¹³C excess atom fraction (EAF) for each bacterial ASV using R version 4.0.3 [37] and data.table [38] with custom scripts available at https://www.github.com/bramstone/. Negative enrichment values were corrected using previously published methods [17]. ASVs that appeared in less than two of the three replicates of an ecosystem-treatment combination (n = 3) and less than three density fractions within those two replicates were removed to avoid assigning spurious estimates of isotope enrichment to infrequent taxa. Any ASVs filtered out of one ecosystem-treatment group were allowed to be present in another if they met the frequency threshold. Applying these filtering criteria, we limited our analysis towards 3759 unique bacterial ASVs which accounted for a small proportion of the total diversity but represented 68.0% of all sequence reads, and encompassed most major bacterial groups (Supplementary Fig. 1).

Analysis of life history strategies and nutrient response

All statistical tests were conducted in R version 4.0.3 [37]. We assessed the ability of phylum-level assignment of life history strategy to predict growth in response to C and N addition, as proxied by the incorporation of heavy isotope during DNA replication [39, 40]. Phylum-level assignments (Table 1) were based on the most frequently observed behavior of lineages with a representative phylum (or subphylum) as compiled previously [23]. We averaged ¹⁸O EAF values of bacterial taxa for each treatment and ecosystem and then subtracted the values in control soils from values in C-amended soils to determine C response (∆¹⁸O EAF_C) and from the ¹⁸O EAF of bacteria in CN-amended soils to determine C and N response (Δ¹⁸O EAF_CN). Because an ASV must have a measurable EAF in both the control and treatment for a valid Δ¹⁸O EAF to be calculated, we were only able to resolve the nutrient response for 2044 bacterial ASVs – 1906 in response to C addition and 1427 in response to CN addition.

We used Gaussian finite mixture modeling, as implemented by the mclust R package [41], to demarcate plausible multi-isotopic signatures for oligotrophs and copiotrophs. For each treatment, we calculated average per-taxon ¹³C and ¹⁸O EAF values. To compare both isotopes directly, we divided ¹⁸O EAF values by 0.6 based on the estimate that this value (designated as µ) represents the fraction of oxygen atoms in DNA derived from the ¹⁸O-water, rather than from ¹⁶O within available C sources [42]. Two mixture components, corresponding to oligotrophic and copiotrophic growth modes, were defined using the Mclust function using ellipsoids of equal volume and shape. We observed several microorganisms with high ¹⁸O enrichment but comparatively low ¹³C enrichment, potentially indicating growth following the depletion of the added glucose, and that were reasonably clustered as oligotrophs in our mixture model.

We tested how frequently mixture model clustering of each microorganism’s growth (based on average ¹⁸O–¹³C EAF in a treatment) could predict its growth across replicates (n = 12 in each treatment—although individual). We applied the treatment-level mixture models defined above to the per-taxon isotope values in each replicate, recording when a microorganism’s life history strategy in a replicate agreed with the treatment-level cluster, and when it didn’t. We used exact binomial tests to test whether the number of “successes” (defined as a microorganism being grouped in the same life history category as its treatment-level cluster) was statistically significant. To account for type I error across all individual tests (one per ASV per treatment), we adjusted P values in each treatment using the false-discovery rate (FDR) method [43].

To determine the extent that life history categorizations may be appropriately applied at finer levels of taxonomic resolution, we constructed several hierarchical linear models using the lmer function in the nlme package version 3.1-149 [44]. To condense growth information from both isotopes into a single analysis, ¹⁸O and ¹³C EAF values were combined into a single variable using principal components analysis separately for each treatment. Across the C and CN treatments, the first principal component (PC1) was able to explain – respectively – 86% and 91% of joint variation of ¹⁸O and ¹³C EAF values. In all cases, we applied PC1 as the response variable and treated taxonomy and ecosystem as random model terms to limit the potential of pseudo-replication to bias significance values. We used likelihood ratio analysis and Akaike information criterion (AIC) values to compare models where life history strategy was determined based on observed nutrient responses at different taxonomic levels (Eq. 1) against a model with the same random terms but without any life history strategy data (Eq. 2). Separate models were applied to each treatment. To reduce model overfitting, we removed families represented by fewer than three bacterial ASVs as well as phyla represented by only one order. In addition, we removed bacterial ASVs with unknown taxonomic assignments (following Morrissey et al. [21]). This limited our analysis to 1 049 ASVs in the C amendment and 984 in the CN amendment.

Here, life history strategy was defined at each taxonomic level using the mixture models above and based on the mean ¹⁸O and ¹³C EAF values of each bacterial lineage (Supplemental Fig. 2). We compared these models with the no-strategy model (Eq. 2) directly using likelihood ratio testing.

Source: Ecology - nature.com