Nutrients cause consolidation of soil carbon flux to small proportion of bacterial community

Sample collection and incubation

Three replicates of soil samples were collected from the top 10 cm in of plant-free patches in four ecosystems along the C. Hart Merriam elevation gradient in Northern Arizona²⁵ beginning at high desert grassland (1760 m), and followed at higher elevations by piñon-pine juniper woodland (2020 m), ponderosa pine forest (2344 m), and mixed conifer forest (2620 m). Soils 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. Soil incubations were performed on soils with mass of 20 g of dry soil for measurements of CO₂ and microbial biomass carbon (MBC), while 2 g of dry soil aliquots were incubated separately (but under equivalent conditions) for quantitative stable isotope probing (qSIP). We applied three treatments to these soils through the addition of water (up to 70% water-holding capacity): water alone (control), with glucose (C treatment; 1000 µg C g⁻¹ dry soil), or with glucose and nitrogen (C + N treatment; [NH₄]₂SO₄ at 100 µg N g⁻¹ dry soil). All samples for qSIP were incubated with ¹⁸O-enriched water (97 atom%) and matching controls necessary to calculate the change in ¹⁸O enrichment across the microbial community. We applied water at natural abundance (i.e., no ¹⁸O-enriched water) to the larger soil samples prepared for measurement of carbon flux. All soils were incubated in the dark for one week. Following incubation, soils were frozen at −80 °C for 1 week prior to DNA extraction.

Soil, CO₂, and microbial biomass measurements

We analyzed headspace gas of soils for CO₂ concentration and δ¹³CO₂ three times during the week-long incubation using a LI-Cor 6262 (LI-Cor Biosciences Inc. Lincoln, NE, USA) and a Picarro G2201 (Picarro Inc., Sunnyvale, CA, USA), respectively. Prior to incubation we analyzed soil MBC using the chloroform-fumigation extraction method on 10 g of soil. One sub-sample was immediately extracted with 25 ml of a 0.05 M K₂SO₄ solution, while a second sub-sample was first fumigated with chloroform (for 5 days), after which it was similarly extracted. Following K₂SO₄ addition, we agitated soils for 1 h, filtered the extract through a Whatman #3 filter paper, and dried the filtered solution (60 °C, 4 days). Salts with extracted C were ground and analyzed for total C using an elemental analyzer coupled to a mass spectrometer. MBC was calculated as the difference between the fumigated and immediately extracted samples’ soil C using an extraction efficiency of 0.45 (as per Liu et al.²⁶).

Quantitative stable isotope probing

We performed DNA extraction and 16S amplicon sequencing on ¹⁸O-incubated qSIP soils^11,12,13. The procedures targeted the V4 region of the 16S gene as specified by the Earth Microbiome Project (EMP, http://www.earthmicrobiome.org) standard protocols^27,28_. We used PowerSoil DNA extraction kits following manufacture instructions to isolate DNA from soil (MoBio laboratories, Carlsbad, CA, USA). We quantified extracted DNA using the Qubit dsDNA High-Sensitivity assay kit and a Qubit 2.0 Fluorometer (Invitrogen, Eugene, OR, USA). To quantify the degree of ¹⁸O isotope incorporation into bacterial DNA, we performed density fractionation and sequenced 15–18 fractions separately following methods modified from the canonical publication⁷. We added 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. We separated each sample from the continuous gradient into approximately 20 fractions (150 µL) using a modified fraction recovery system (Beckman Coulter). We then measured the density of each separate fraction with a Reichart AR200 digital refractometer (Reichert Analytical Instruments, Depew, NY, USA) and retained fractions with densities between 1.640 and 1.735 g cm⁻³. We cleaned and purified DNA in these fractions using isopropanol precipitation, quantified DNA using the Quant-IT PicoGreen dsDNA assay (Invitrogen) and a BioTek Synergy HT plate reader (BioTek Instruments Inc., Winooski, VT, USA), and quantified bacterial 16S gene copies using qPCR (primers: Supplementary Table 1) in triplicate. We used 8 µL reactions 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. We amplified DNA using a Bio-Rad CFX384 Touch real-time PCR detection system (Bio-Rad, Hercules, CA, USA) with the following cycling conditions: 95 °C at 1 min and 44 cycles of 95 °C (30 s), 64.5 °C (30 s), and 72 °C (1 min).

We sequenced the 16S V4 region (primers: EMP standard 515F—806R; Supplementary Table 1) on an Illumina MiSeq (Illumina, Inc., San Diego, CA, USA). Sequences were amplified using the same reaction mix as qPCR amplification but cycling at 95 °C for 2 min followed by 15 cycles of 95 °C (30 s), 55 °C (30 s), and 60 °C (4 min). In addition to post-incubation soils, we extracted, amplified, and sequenced DNA of the bacterial community at the start of the incubation.

Sequence processing and qSIP analysis

The raw sequence data of forward and reverse reads (FASTQ) were processed within the QIIME 2 environment (release 2018.6)^29,30, denoising sequences with the available DADA2 pipeline³¹. We clustered the remaining sequences into amplicon sequence variants or ASVs (at 100% sequence identity) against the SILVA 132 database³² using an open-reference Naïve Bayes feature classifier³³. We removed global singletons and doubleton ASVs, non-bacterial lineages, and samples with less than 4000 sequence reads. Removal of global singletons and doubletons resulted in the removal of 2241 unique ASVs from the feature table yielding 115,647 out of 117,888 (a retention of 98% of all ASVs) as well as the loss of 4018 sequences leaving 37,765,678 (a retention >99% of all sequences). We combined taxonomic information and ASV sequence counts with per-fraction qPCR and density measurements using the phyloseq package (version 1.24.2), in R (version 3.5.1)³⁴. Because high-throughput sequencing produces relativized measures of abundance, we converted ASV sequencing abundances in each fraction to the number of 16S rRNA gene copies per g dry soil based on the known amount of dry soil added and the amount of DNA in each soil sample. All data and analytical code have been made publicly accessible³⁵.

To perform qSIP analysis and calculate per-capita growth rates of each ASV, we used our in-house qsip package (https://github.com/bramstone/qsip) based on previously published research^7,10. Because rare and infrequent taxa are more likely to be lost in samples with poor sequencing depth with their absences affecting DNA density changes, we invoked a presence or absence-based filtering criteria on ASVs prior to calculation of per-capita growth rates. Within each ecosystem, we kept only ASVs that appeared in two of the three replicates of a treatment (¹⁸O, C, and C + N) and at that appeared in at least five of the fractions within each of those two replicates. ASVs filtered out of one treatment were allowed to appear in another if they met the frequency threshold.

For all remaining ASVs (1081 representing less than 1% of all ASVs but 58% of all sequence reads), we calculated per-capita gross growth (i.e., cell division) rates observed in each replicate using an exponential growth model¹⁰. We applied these per-capita rates to the number of 16S rRNA gene copies to estimate the production of new 16S rRNA gene copies of each ASV per g dry soil per week using the following equation:

Where N_i,t is the number of 16S rRNA gene copies of taxon i at time t (here after 7 days) and g_i represents the per-capita growth rate (calculated as a daily rate). See Supplementary Fig. 3 for results on the production of 16S gene copies.

Calculation of 16S rRNA gene copy numbers and cell mass

In parallel to taxonomic assignment, we compared quality-filtered 16S sequences against a database of 12,415 complete prokaryote genomes obtained from GenBank. From these genomes, we extracted data on 16S rRNA gene copy number, total genome size, and 16S gene sequence. We used BLAST to find matches against this database to the ASVs generated from QIIME 2 to make per-taxon assignments of 16S rRNA gene copy number and total genome size¹³. For ASVs that did not find an exact match, we assigned 16S rRNA gene copy number values and genome sizes based on the median values observed in the most specific possible taxonomic rank. We estimated the mass of individual cells for each population using published allometric scaling relationships between genome length and cellular mass from West and Brown:³⁶

where M_i indicates cellular mass (g) and G_i indicates genome length (bp) for taxon i. We obtained this relationship by digitizing Fig. 4³⁶ using DataThief III and re-fitting the trend line in log–log space. We estimated that 20% of the cellular mass was carbon³⁷. To validate this approach, cellular mass estimates and initial 16S copy number measurements were used to estimate population-level biomass C values which were summed and compared to initial community-level MBC. We found that these values overestimated initial MBC by an order of magnitude. As such, cellular carbon mass was divided by 10 in our final calculations. We applied cellular mass and 16S copy number estimates to the production of 16S copies to estimate the production of biomass carbon for each taxon during the incubation period (t):

where P_i indicates production of biomass carbon (µg C g dry soil⁻¹ week⁻¹) and C_i indicates 16S copy number per cell for taxon i. The 0.2 coefficient represents an estimate that 20% of cellular mass is composed of carbon.

Efficiency and respiration modeling

We estimated rates of respiration using qSIP-informed growth rates and community-level carbon use efficiency (CUE). CUE estimates were based on the incorporation of ¹⁸O-water into DNA as a measure of gross biomass production^38,39 and measured CO₂ in headspace gas from soil incubations. We calculated the production of ¹⁸O-labeled biomass carbon (¹⁸P) at the community-level for each sample by summing the products of per-taxon ¹⁸O enrichment (excess atom fraction, EAF) and relative abundance:

where ¹⁸P indicates the gross production of ¹⁸O-labeled microbial biomass carbon per gram of dry soil per week, ¹⁸EAF_i indicates the enrichment of DNA of taxon i and y_i indicates its relative abundance, DNA₀ indicates the concentration of DNA per gram of dry soil prior to incubation, and MBC₀ indicates the microbial biomass carbon per gram of dry soil prior to incubation. Here, the MBC₀ ~ DNA₀ function indicates the linear relationship between MBC and DNA concentration. We used the output from Eq. 4 to calculate community CUE for each sample:

where R indicates the total CO₂ respired per gram dry soil per week.

We used the community CUE values from each sample (Eq. 5) to constrain/as upper and lower limits our estimates of per-taxon CUE. For a group of three replicates from a given ecosystem and treatment, we used the minimum and maximum observed community-level CUE values as the acceptable range of per-taxon CUE values. These constraints were used to control the shape of the function of per-taxon CUE and growth rate, though functions were modeled both with and without constraints (i.e., per-taxon CUE values were bounded only by 0 and 0.7). The range of community-level CUE values for each treatment were 0.18–0.53 for control soils, 0.04–0.13 for carbon amended soils and 0.03–0.08 for carbon and nitrogen amended soils and did not vary much between ecosystems. As a result of uncertainty in the literature about the relationship between growth rate and CUE¹⁴, several different relationships were postulated to model per-taxon CUE as a function of per-taxon growth rate: linear increase, linear decrease, exponential decrease, unimodal with peak CUE at growth rate of 0.5, and unimodal with peak CUE at a growth rate of 0.05 (the median of all per-taxon growth rates in the data). Comparisons between functions were made by calculating AIC values from per-taxon respiration, summed, and regressing against measured respiration values. Likewise, for each function, we tested how well per-taxon CUE estimates reconstructed community-level CUE by weighting the CUE value of each taxon by its relative abundance, summing, and regressing against community-level CUE. To select the best per-taxon CUE function, AIC values from both scaling efforts were combined. To make AIC values comparable, all respiration and CUE terms were z-transformed prior to regression scaling. To reflect our priority of estimating per-taxon respiration, AIC values from the respiration scaling regression models were multiplied by two and summed with AIC values from CUE scaling such that AIC_Total = 2(AIC_Resp) + AIC_CUE. Across these comparisons, the best estimate of per-taxon CUE was the unimodal function of growth rate, constrained by community-level CUE and peaking at growth rates of 0.5 (Table 1), such that:

where CUE_i indicates per-taxon CUE, CUE_E:T:max indicates the maximum CUE values observed for a group of replicates within a given ecosystem and treatment (E:T). With this function, higher per-capita growth rate values were parameterized to produce higher CUE values initially and then decrease reflecting a growth-CUE tradeoff¹⁴, here bound by the difference in maximum and minimum CUE values. We applied per-taxon CUE estimates from Eq. 6 to per-taxon growth rates to yield estimates of per-taxon respiration:

where r_i indicates per-capita respiration for taxon i, r_g,i indicates growth-related respiration, r_m,i indicates maintenance-related respiration, and β is a constant of 0.01 that represents the maintenance requirements as a proportion of total energy use⁴⁰. We used these values of per-taxon, per-capita respiration rates to estimate per-taxon respiration per gram of dry soil per week:

where R_i indicates respiration of CO₂–C (µg C g dry soil⁻¹ week⁻¹) for taxon i.

In addition to per-taxon respiration estimates based on ¹⁸O enrichment, we used another model for comparison. Here, respiration was calculated based on 16S abundance alone:

where N_i indicates final 16S abundance for taxon i, R indicates microbial respiration of CO2-C (µg C g dry soil⁻¹ week⁻¹) and N indicates total 16S abundance at the end of the incubation. Here, the R ~ N function indicates the linear relationship, with an intercept of 0, between CO₂ respiration and 16S gene concentration across all samples.

Diversity, compositional, and statistical analysis

For patterns of evenness in bacterial carbon use and relative abundance, we used Pielou’s evenness which is the quotient of Shannon’s diversity and the observed richness. For each sample, we applied Pielou’s evenness to bacterial abundances as well as bacterial carbon use (relativized to sum to one, in both cases).

We created a linear mixed model to test the relationship between the carbon use (the sum of biomass production and respiration) and relative abundance of bacterial genera from the dominant phyla, which accounted for >90% of all C flux. Here, we averaged carbon use and relative abundance for all replicates in a given ecosystem and treatment. We used the lme4 R package (version 1.1-20)⁴¹ and obtained p-values using the Satterthwaite method in the lmerTest R package (version 3.1-0)⁴². To limit pseudo-replication, we accounted for differences in carbon use across ecosystems and due to bacterial Genus by implementing random intercepts. We selected for the optimal random and fixed components by dropping individual terms and comparing models with likelihood ratio tests, disregarding models that failed to converge. Our final model fit was:

where C_i indicates the relativized carbon use for taxon i (averaged across all three replicates in a given ecosystem and treatment), y_i indicates the relative abundance of taxon i (averaged across all three replicates), T indicates soil treatment, and E indicates ecosystem.

For differences in composition, we created species abundance tables using the traditional abundances, as well as measures of carbon use (growth and maintenance respiration) of each ASV in each sample. To account for differences in absolute abundances and flux rates between sites, we relativized all abundance tables. We summarized compositional differences using Bray–Curtis dissimilarities then identified multivariate centroids for all replicates in a site by treatment group. We tested the effect of site and nutrient amendment on the resulting group centroids using PERMANOVA tests implemented with the adonis function in the vegan package (version 2.5-3)⁴³. We related compositional shifts in relative abundance to those in relativized growth and maintenance using Mantel tests with the mantel function in vegan.

To test for changes in the type of soil C preferred by microbial genera (either ¹³C-labeled glucose or ¹²C soil carbon) in response to nitrogen addition, we used Levene’s test with the car package (version 3.0-10)⁴⁴. Specifically, we analyzed the relationship between ¹³C use and ¹²C use (both relativized) on bacterial genera across all replicates and in C and C + N treatments using a linear model. We then extracted model residuals and tested whether variance was significantly different across treatments by focusing on the interaction between individual replicates and treatment. This produced a significance test describing treatment-level differences in ¹³C–¹²C use.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com