Nutrients cause consolidation of soil carbon flux to small proportion of bacterial community
Sample collection and incubationThree 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 Arizona25 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 CO2 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−1 dry soil), or with glucose and nitrogen (C + N treatment; [NH4]2SO4 at 100 µg N g−1 dry soil). All samples for qSIP were incubated with 18O-enriched water (97 atom%) and matching controls necessary to calculate the change in 18O enrichment across the microbial community. We applied water at natural abundance (i.e., no 18O-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, CO2, and microbial biomass measurementsWe analyzed headspace gas of soils for CO2 concentration and δ13CO2 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 K2SO4 solution, while a second sub-sample was first fumigated with chloroform (for 5 days), after which it was similarly extracted. Following K2SO4 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.26).Quantitative stable isotope probingWe performed DNA extraction and 16S amplicon sequencing on 18O-incubated qSIP soils11,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 protocols27,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 18O isotope incorporation into bacterial DNA, we performed density fractionation and sequenced 15–18 fractions separately following methods modified from the canonical publication7. 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−3. 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−1 Phusion HotStart II Polymerase (Thermo Fisher Scientific, Waltham, MA), 1× Phusion HF buffer (Thermo Fisher Scientific), 3.0 mM MgCl2, 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 analysisThe 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 pipeline31. We clustered the remaining sequences into amplicon sequence variants or ASVs (at 100% sequence identity) against the SILVA 132 database32 using an open-reference Naïve Bayes feature classifier33. 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)34. 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 accessible35.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 research7,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 (18O, 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 model10. 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:$$frac{{rm{d}}{N}_{{rm{i}}}}{{{rm{d}}t}}={N}_{{rm{i,t}}}-{N}_{{rm{i,t}}}{e}^{-{g}_{{rm{i}}}t},$$
(1)
Where Ni,t is the number of 16S rRNA gene copies of taxon i at time t (here after 7 days) and gi 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 massIn 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 size13. 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:36$${{{log }}}_{10}({M}_{{rm{i}}})=frac{{{{log }}}_{10}left({G}_{{rm{i}}}right)-9.4}{0.24},$$
(2)
where Mi indicates cellular mass (g) and Gi indicates genome length (bp) for taxon i. We obtained this relationship by digitizing Fig. 436 using DataThief III and re-fitting the trend line in log–log space. We estimated that 20% of the cellular mass was carbon37. 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):$${P}_{{rm{i}}}=frac{{rm{d}}{N}_{{rm{i}}}/{{rm{d}}t}}{C_{{rm{i}}}}cdot {M}_{{rm{i}}}cdot 0.2,$$
(3)
where Pi indicates production of biomass carbon (µg C g dry soil−1 week−1) and Ci 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 modelingWe estimated rates of respiration using qSIP-informed growth rates and community-level carbon use efficiency (CUE). CUE estimates were based on the incorporation of 18O-water into DNA as a measure of gross biomass production38,39 and measured CO2 in headspace gas from soil incubations. We calculated the production of 18O-labeled biomass carbon (18P) at the community-level for each sample by summing the products of per-taxon 18O enrichment (excess atom fraction, EAF) and relative abundance:$${, }^{18}{P}=mathop{sum }limits_{i=1}^{n}({,}^{18}{{{rm{EAF}}}}_{{rm{i}}}cdot {y}_{{rm{i}}})cdot {rm{DN}}{rm{A}}_{0}cdot fleft({{rm{MB}}}{rm{C}}_{0} sim {rm{DN}}{rm{A}}_{0}right),$$
(4)
where 18P indicates the gross production of 18O-labeled microbial biomass carbon per gram of dry soil per week, 18EAFi indicates the enrichment of DNA of taxon i and yi indicates its relative abundance, DNA0 indicates the concentration of DNA per gram of dry soil prior to incubation, and MBC0 indicates the microbial biomass carbon per gram of dry soil prior to incubation. Here, the MBC0 ~ DNA0 function indicates the linear relationship between MBC and DNA concentration. We used the output from Eq. 4 to calculate community CUE for each sample:$${{rm{CUE}}}=frac{{,}^{18}{{P}}}{(!{,}^{18}P+R)},$$
(5)
where R indicates the total CO2 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 CUE14, 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 AICTotal = 2(AICResp) + AICCUE. 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:$${{rm{CUE}}}_{{rm{i}}}=-4({{rm{CUE}}}_{{rm{E}}{rm{:}}{rm{T}}{rm{:}}{{rm{range}}}})cdot {left({g}_{{rm{i}}}-0.5right)}^{2}+({{rm{CUE}}}_{{rm{E}}{rm{:}}{rm{T}}{rm{:}}{max }}),$$
(6)
where CUEi indicates per-taxon CUE, CUEE: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 tradeoff14, 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:$${r}_{{rm{i}}}={r}_{{rm{g,i}}}+{r}_{{rm{m,i}}}=left(frac{{g}_{{rm{i}}}}{{{rm{CUE}}}_{{rm{i}}}}-{g}_{{rm{i}}}right)+left(frac{{g}_{{rm{i}}}}{{{rm{CUE}}}_{{rm{i}}}}-{g}_{{rm{i}}}right)cdot beta,$$
(7)
where ri indicates per-capita respiration for taxon i, rg,i indicates growth-related respiration, rm,i indicates maintenance-related respiration, and β is a constant of 0.01 that represents the maintenance requirements as a proportion of total energy use40. We used these values of per-taxon, per-capita respiration rates to estimate per-taxon respiration per gram of dry soil per week:$${R}_{{rm{i}}}={P}_{{rm{i}}}cdot {r}_{{{rm{g,i}}}}+{P}_{{rm{i}}}cdot {r}_{{{rm{m,i}}}},$$
(8)
where Ri indicates respiration of CO2–C (µg C g dry soil−1 week−1) for taxon i.In addition to per-taxon respiration estimates based on 18O enrichment, we used another model for comparison. Here, respiration was calculated based on 16S abundance alone:$${R}_{{rm{i}}}={N}_{{rm{i}}}cdot f(R sim N+0),$$
(9)
where Ni indicates final 16S abundance for taxon i, R indicates microbial respiration of CO2-C (µg C g dry soil−1 week−1) 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 CO2 respiration and 16S gene concentration across all samples.Diversity, compositional, and statistical analysisFor 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)41 and obtained p-values using the Satterthwaite method in the lmerTest R package (version 3.1-0)42. 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:$${{{log }}}_{10}({C}_{{rm{i}}}) sim {{{log }}}_{10}left({y}_{{rm{i}}}right)ast T+left(1|Eright)+(1|{{rm{Genus}}}),$$
(10)
where Ci indicates the relativized carbon use for taxon i (averaged across all three replicates in a given ecosystem and treatment), yi 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)43. 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 13C-labeled glucose or 12C soil carbon) in response to nitrogen addition, we used Levene’s test with the car package (version 3.0-10)44. Specifically, we analyzed the relationship between 13C use and 12C 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 13C–12C use.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
