Elevated temperature and CO2 strongly affect the growth strategies of soil bacteria
Study site, field experiment, and soil samplingThis study relied on an experimental field station that simulates atmospheric CO2 enrichment and warming (Fig. 1a). The experimental site is located at Kangbo village (31°30′N, 120°33′E), Guli Township, Changshu municipality in Jiangsu Province, China. The climate type is a subtropical monsoon climate with a mean annual precipitation between 1100–1200 mm and an annual average temperature of approximately 16 °C. The high rainfall and temperature mainly occur from May through September25. The soils are Gleyic Stagnic Anthrosols derived from clayey lacustrine deposits. The main properties of the topsoil (0–15 cm) before the field experiment were as follows: pH (H2O) 7.0, bulk density 1.2 g cm−3, and contents of organic C and total N of 16.0 g kg−1 and 1.9 g kg−1, respectively.The FACE system has been well described in previous publications26,53. The climate change treatments were set according to the representative concentration pathway (RCP) scenario that modeled CO2 atmospheric concentration and temperature elevation in approximately 30–40 years. Elevated atmospheric CO2 concentration and warming of crop canopy air were maintained 24 h a day during the crop-growing period. Each treatment was replicated in three rings with the same infrastructure, and the rings were arranged in a split row design (Fig. 1b). All rings were buffered by adjacent open fields to avoid any treatment cross-over. Rice and wheat cultivations of all plots were managed with local conventional practices. Soil samples for qSIP incubation were collected in June 2020.Evaluation of nutrient pulses in soil after water additionTo determine if there is a nutrient pulse after soil rewetting, we estimated the dynamics of several biochemical characteristics (i.e., CO2 fluxes, hydrolase activity, DOC, DN, TC, and TN) in soils before and during the incubation. Soil samples for nutrient flux measurements were collected from the free-air CO2 enrichment and warming experimental station in July 2022.To measure the CO2 production rate, soil (2.00 g) was set in a 29 mm diameter glass vial (Volume 23 mL). Water (400 μL) was added evenly to soil in each vial. Gas samples were then collected from each incubation vial at multiple time points (i.e., 3 h, 6 h, 9 h, 12 h, 24 h, 72 h, and 144 h after water addition), and CO2 concentrations were determined in all gas samples with a Trace Gas Chromatograph (Agilent 7890, Santa Clara, CA, USA). Other biochemical characteristics were measured in parallel incubations, with destructive sampling. Briefly, soil (40.00 g) was set in a plastic jar, and 8 mL water was added evenly to the soil. All the incubation vials and jars were placed in the dark at 25 °C and soil moisture content was adjusted twice a day. Destructive sampling was conducted after 1, 3, and 6 days of incubation. Including prewet soil (the soil before water addition), a total of 48 soil samples were obtained. Soil samples were weighed before and after being oven dried at 105 °C, and soil moisture content was measured using the gravimetric method. Fluorescein Diacetate (FDA) hydrolase activity of soil was measured by using soil FDA hydrolase activity assay kit (Solarbio®, Beijing, China) according to the manufacturer’s protocol. Soil total carbon and total nitrogen were measured with a C/N elemental analyzer (multi EA® 5000, analytikjena, Germany). Dissolved organic carbon and dissolved nitrogen were analyzed using a TOC/TN analyzer (multi N/C® 3100, analytikjena, Germany) after extraction with distilled water.Simulating nutrient pulse after rewetting dry soil by adding 18O-waterTo determine whether and how microbial growth varied in response to pulse events after long-term acclimation to warming and CO2 enrichment, we estimated the population-specific growth rates of active microbes by conducting a 18O-water incubation experiment combined with DNA quantitative stable isotope probing (DNA-qSIP) (Fig. 1c). The incubation conditions were similar to those reported in a previous study6. In brief, approximately 60 g fresh soil of each treatment were sieved (2 mm) and air-dried (24 h at room temperature) immediately after transport to the laboratory. Then, triplicate samples of dry soils (2.00 g) were incubated in the dark at room temperature in sterile plastic aerobic culture tubes (17 × 100 mm) with 400 μL of 98 atom% H218O or natural abundance water (H216O) for 6 days, with harvests at four time points (T = 0, 1, 3, 6 d after water addition). DNA was extracted from the soils of 0 d incubation treatment immediately after water addition (~30 s interval), representing the prewet treatment. At each harvest, soils were destructively sampled and immediately stored at −80 °C. A total of 84 soil samples (4 climate change treatments × 3 replicates at 0 h with H216O addition + 4 treatments × 3 subsequent time points × 3 replicates × 2 types of H2O addition) were collected.DNA extraction and isopycnic centrifugationTotal DNA from all the collected soil samples was extracted using the FastDNA™ SPIN Kit for Soil (MP Biomedicals, Cleveland, OH, USA) according to the manufacturer’s instructions. The concentration of extracted DNA was determined fluorometrically using Qubit® DNA HS (High Sensitivity) Assay Kits (Yeasen Biotechnology, Shanghai, China) on a Qubit® 4 fluorometer (Thermo Scientific™, Waltham, MA, USA). The DNA samples of day 1, day 3, and day 6 were used for isopycnic centrifugation, and the detailed pipeline was described previously with minor modifications6. Briefly, 3 μg DNA were added into 1.85 g mL−1 CsCl gradient buffer (0.1 M Tris-HCl, 0.1 M KCl, 1 mM EDTA, pH = 8.0) with a final buoyant density of 1.718 g mL−1. Approximately 5.1 mL of the solution was transferred to an ultracentrifuge tube (Beckman Coulter QuickSeal, 13 mm × 51 mm) and heat-sealed. All tubes were spun in an Optima XPN-100 ultracentrifuge (Beckman Coulter) using a VTi 65.2 rotor at 177,000 × g at 18 °C for 72 h with minimum acceleration and braking.Immediately after centrifugation, the contents of each ultracentrifuge tube were separated into 20 fractions (~250 μL each fraction) by displacing the gradient medium with sterile water at the top of the tube using a syringe pump (Longer Pump, LSP01‐2A, China). The buoyant density of each fraction was measured using a digital hand-held refractometer (Reichert, Inc., Buffalo, NY, USA) from 10 μL volumes. Fractionated DNA was precipitated from CsCl by adding 500 μL 30% polyethylene glycol (PEG) 6000 and 1.6 M NaCl solution, incubated at 37 °C for 1 h, and then washed twice with 70% ethanol. The DNA of each fraction was then dissolved in 30 μL of Tris‐EDTA buffer. Detailed information (company names and catalog numbers) of the reagents and consumables for qSIP experiment is provided in Supplementary Data 2.Quantitative PCR and sequencingTotal 16S rRNA gene copies for DNA samples of all the fractions and the day-0 soils were quantified using the primers for V4-V5 regions: 515F (5′‐GTG CCA GCM GCC GCG G‐3′) and 907R (5′‐CCG TCA ATT CMT TTR AGT TT‐3′)54. Plasmid standards were prepared by inserting a copy of purified PCR product from soil DNA into Escherichia coli. The E. coli was then cultured, followed by plasmid extraction and purification. The concentration of plasmid was measured using Qubit DNA HS Assay Kits. Standard curves were generated using 10‐fold serial dilutions of the plasmid. Each reaction was performed in a 25-μL volume containing 12.5 μL SYBR Premix Ex Taq (TaKaRa Biotechnology, Otsu, Shiga, Japan), 0.5 μL of forward and reverse primers (10 μM), 0.5 μL of ROX Reference Dye II (50×), 1 μL of template DNA, and 10 μL of sterile water. A two-step thermocycling procedure was performed, which consisted of 30 s at 95 °C, followed by 40 cycles of 5 s at 95 °C, 34 s at 60 °C (at which time the fluorescence signal was collected). Following qPCR cycling, melting curves were obtained to ensure that the results were representative of the target gene. Average PCR efficiency was 96% and the average slope was −3.38, with all standard curves having R2 ≥ 0.99.The DNA of day-0 samples (unfractionated) and the fractionated DNA of fractions with buoyant density between 1.695 and 1.735 g/mL were selected for 16S rRNA amplicon sequencing by using the same primers of qPCR (i.e., 515F/907R). Eleven out of 20 fractions from each ultracentrifuge tube (density between 1.695 and 1.735 g/mL) were selected because they contained more than 99% gene copy numbers of the 20 fractions. A total of 804 DNA samples (12 unfractionated DNA samples + 72 × 11 fractionated DNA samples) were sequenced using the NovaSeq6000 platform (Genesky Biotechnologies, Shanghai, China).The sequences were quality-filtered using the USEARCH v.11.055. In brief, sequences 0.5 were removed. Chimeras were identified and removed. Subsequently, high-quality sequences were clustered into operational taxonomic units (OTUs) using the UPARSE algorithm at a 97% identity threshold, and the most abundant sequence from each OTU was selected as a representative sequence. The taxonomic affiliation of the representative sequence was determined using the RDP classifier (version 16)56. In total, 51,127,459 reads of the bacterial 16S rRNA gene and 11,898 OTUs were obtained. The 16S rRNA amplicon sequences were uploaded to the National Genomics Data Center (NGDC) Genome Sequence Archive (GSA) with accession number CRA006507.Quantitative stable isotope probing calculationsWe used the amount of 18O incorporated into DNA to estimate the growth rates of active taxa24,57. The density shifts of OTUs between 16O and 18O treatments were calculated following the qSIP procedures24. Briefly, the number of 16S rRNA gene copies per OTU in each density fraction was calculated by multiplying the OTU’s relative abundance (acquisition by sequencing) by the total number of 16S rRNA gene copies (acquisition by qPCR). Then, the GC content and molecular weight of a particular taxon were calculated. Further, the change in 18O isotopic composition of 16S rRNA genes for each taxon was estimated. We assumed an exponential growth model over the course of the incubations, and absolute population growth rates were estimated over each of the three time intervals of the incubation: 0–1, 0–3, and 0–6 d corresponding to the samplings at days 1, 3, and 6. The absolute growth rate is a function of the rate of appearance of 18O-labeled 16S rRNA genes. Therefore, the growth rate (g) of taxon i was calculated as:$${g}_{i}={{{{{rm{ln}}}}}}left(frac{{N}_{{{{{{rm{TOTAL}}}}}}{it}}}{{N}_{{{{{{rm{LIGHT}}}}}}{it}}}right)times frac{1}{t}$$
(1)
Where NTOTALit is the number of total gene copies for taxon i and NLIGHTit represents the unlabeled 16S rRNA gene abundances of taxon i at the end of the incubation period (time t). NLIGHTit is calculated by a function with four variables: NTOTALit, molecular weights of DNA (taxon i) in the 16O treatment (MLIGHTi) and in the 18O treatment (MLABi), and the maximum molecular weight of DNA that could result from assimilation of H218O (MHEAVYi)24. We further calculated the average growth rates (represented by the production of new 16S rRNA gene copies of each taxon per g dry soil per day) over each of the three time intervals of the incubation: 0–1, 0–3, and 0–6 d, using the following equation39:$$frac{d{N}_{i}}{{dt}}={N}_{{{{{{rm{TOTAL}}}}}}{it}}left(1-{e}^{-{g}_{i}t}right)times frac{1}{t}$$
(2)
Where t is the incubation time (d). All data calculations were performed using the qSIP package (https://github.com/bramstone/qsip) in R (v. 3.6.2).Grouping of taxa into growth strategiesWe compared the average growth rates of taxa at three time intervals (n = 3 in each time interval) and classified the species into rapid, intermediate, and slow growth strategies based on the timing of the maximum growth rate (Fig. 1d): (1) Rapid responders: Species had the highest growth rates by 1 day of the incubation; (2) Intermediate responders: Species had the highest growth rates at the 3-day incubation; (3) Slow responders: Species had the highest growth rates at the 6-day incubation. The taxa with growth rates significantly greater than zero can be divided into one of three strategies in each treatment.Analyses of phylogenetic conservationPhylogenetic tree analyses were performed in Galaxy/DengLab (http://mem.rcees.ac.cn:8080/) with PyNAST Alignment and FastTree functions58,59. The trees were visualized and edited using iTOL60. To estimate the phylogenetic patterns of growth strategies, phylogenetic dispersion (D) was calculated by the function ‘phylo.d’ of package “caper” in R (v. 3.6.2). The D statistic equal to 1 means the observed trait has a phylogenetically random distribution across the tips of the phylogeny, and 0 means the observed trait is as clumped as if it had evolved by Brownian motion29. Increasing phylogenetic clumping in the binary trait is indicated by values of D decreasing from 1. The p values were obtained by performing 1000 permutations to test D for significant departure from 1 (random distribution). Furthermore, the phylogenetic signal metrics Blomberg’s K and Pagel’s λ, were used to test for significantly nonrandom phylogenetic distributions using the package “phylosignal” in R. The p values of both indices were used to test the significance of phylogenetic signals.The nearest taxon index (NTI) was calculated to determine the degree of phylogenetic clustering as described previously5. The mean nearest taxon distance (MNTD) were calculated in the “picante” package of R61. The values of NTI are equivalent to the negative output of the standardized effect size (SES) of the observed MNTD distances, which test whether the distribution of growth strategies across different phylogenetic groups is random or nonrandom. NTI values > 0 and their p values < 0.05 represent phylogenetic clustering, while NTI values < 0 and p values > 0.95 indicate phylogenetic over-dispersion. The p values of NTI between 0.05 and 0.95 represent random phylogenetic distributions. The data were converted into binary matrices (1 s and 0 s) before all phylogenetic analyses. For the OTUs that were present in more than one treatment and had differing responses (in the analyses when all climate treatments are combined, i.e., in Table 1 and Supplementary Fig. 6), the OTU was classified to the respective growth responder when that taxa exhibited that specific growth strategy in any treatment.Quantification of explained variance in the distribution patterns of strategiesThe phylogenetic distance matrix obtained from the phylogenetic tree and three binary matrices (including four climate treatments) for three growth strategies were decomposed by the package ‘FactoMineR’ in R. The distribution matrices (binary) of three growth responders in the four climate treatments were used to represent the impact of CO2 and temperature on bacterial growth. To predict the microbial growth strategy (i.e., “y”, a 1 or 0 indicating the membership in a response category), variance partitioning was performed by using the first four phylogenetic principal components (PCs) (accounting for over 80% of phylogenetic variance) and two environment PCs (accounting for over 65% of habitat variance). The fraction of the variance explained by phylogenetic constraints and environmental acclimation was calculated by the package ‘car’ in R.Statistical analysesUncertainty of growth rates (95% confidence interval) was estimated using a bootstrapping procedure with 1000 iterations6. The cumulative growth rates at phylum-level were estimated as the sum of taxon-specific growth rates of those OTUs affiliated to the same phylum. The per capita growth rates of each OTU were calculated by dividing absolute growth rates by the total 16S rRNA gene abundance of taxon i. Ecological processes of active populations (e.g., density-dependence) after rewetting dry soil were estimated by correlating species-specific fitness (refers to per capita growth rates of each OTU) with initial population size using linear regression analyses (R v.3.6.2, ‘lm’ function). The beta diversity of the bacterial community at 0 d incubation was visualized by unconstrained principal coordinate analysis (PCoA) and tested by permutational multivariate analysis of variance (PERMANOVA) based on Bray-Curtis distance, using the vegan package in R (v. 3.6.2)62. Comparisons between climate scenarios were tested by two-way ANOVA and LSD post hoc tests (SPSS 19 for Windows, IBM Corp., Armonk, NY, USA).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More
