Site characteristics and experimental design
In this study, agricultural soils with five SOM contents were collected in 2015 from the following three different locations with the same climate type (the moderate temperate continental climate) in Northeast China (Table S3 and Fig. 1): Bei’an (BA), Hailun (HL), and Dehui (DH). Their MAT and MAP range from 1.0 to 4.4 and 520 to 550, respectively. After collection, the samples were transported to the Hailun Agricultural Ecological Experimental Station (HL), where the samples were packed into the same PVC tubes. Moving the soil from these three initial sampling points to the HL may have had some influence on the microbes, but compared with longer-distance soil translocation across different climatic zones, the HL site can be regarded as an in situ site that reflects the original climatic conditions. The SOM contents were 2%, 3%, 5%, 7%, and 9% (equivalent to 10, 18, 28, 36, and 56 g C kg−1 soil−1, respectively), and all the soils were classified as Mollisols according to the FAO classification. Here, we designed a unique latitudinal soil translocation experiment to investigate the relationship between the bacterial and fungal community stability and the responses of soil C molecular structure to climate warming. The detailed protocol for the experiment was the following: (1) Forty kilograms of topsoil (0–25 cm) was collected for each SOM. The latitude and longitude of the sampling sites and soil geochemical characteristics are shown in Tables S3 and S4. Detailed data can be found in Supplementary Data 1. (2) The soil was homogenized using a 2 mm sieve and filled with sterilized PVC tubes. The PVC tube was 5 cm in diameter at the bottom and 31 cm in height. Each tube was filled with a 25 cm-high soil column, which corresponded to approximately 1 kg of soil. The bottom of the pipe was filled with 1 cm quartz sand, and a 5 cm space was left at the top. (3) From October to November 2015, 90 PVC pipes containing soil (5 SOM gradients × 3 replicates × 6 climatic conditions) were transported to six ecological research stations with different geoclimatic conditions and SOM contents, and 15 PVC pipes were placed in each station. Once the experiment was set up, the weeds growing in each PVC pipe were manually removed every 2–3 weeks to avoid the impact of plants.
The six ecological research stations were the Hailun Agricultural Ecological Experimental Station (HL, N 47°27′, E 126°55′) in Heilongjiang Province, Shenyang Agriculture Ecological Experimental Station (SY, N 41°49′, E 123°33′) in Liaoning Province, Fengqiu Agricultural Ecological Experimental Station (FQ, N 35°03′, E 114°23′) in Henan Province, Changshu Agricultural Ecological Experimental Station (CS, N 31°41′, E 120°41′) in Jiangsu Province, Yingtan Red Soil Ecological Experiment Station (YT, N 28°12′, E 116°55′) in Jiangxi Province and Guangzhou National Agricultural Science and Technology Park (GZ, N 23°23′, E 113°27′) in Guangdong Province. The MAT and MAP at the six ecological research stations ranged from 1.5 to 21.9 °C and from 550 to 1750 mm from north to south, respectively. Details of their climatic conditions (e.g., climatic types) are shown in Table S5. All tubes were removed from each station after 1 year.
The soil samples were stored on dry ice and rapidly transported back to the laboratory. The soil pH was measured by the potentiometric method. Nitrate (NO3−-N) and ammonium nitrogen (NH4+-N) were measured by the Kjeldahl method. DOC was measured using a total organic carbon analyzer (Shimadzu Corporation, Kyoto, Japan). SOC was determined by wet digestion using the potassium dichromate method53. Microbial biomass C (MBC) was measured by the chloroform fumigation-incubation method54. All geochemical attributes are shown in Table S4.
Solid-state 13C NMR analysis of soil C molecular groups
Solid-state 13C NMR spectroscopy analysis was performed to determine the molecular structure of SOC. A Bruker-Avance-iii-300 spectrometer was used at a frequency of 75 MHz (300 MHz 1H). Before the examination, the soil samples were pretreated with hydrofluoric acid to eliminate the interference of Fe3+ and Mn2+ ions in the soil. Specifically, 5 g of air-dried soil was weighed in a 100 ml centrifuge tube with 50 ml of hydrofluoric acid solution (10% v/v) and shaken for 1 h. The supernatant was then removed by centrifugation at 3000 rpm for 10 min. The residues were washed eight times with a hydrofluoric acid solution (10%) with ultrasonication. The oscillation program consisted of the following: four × 1 h, three × 12 h, and one × 24 h. The soil samples were washed with distilled water four times to remove the residual hydrofluoric acid. The above-mentioned treated soil samples were dried in an oven at 40 °C, ground and passed through a 60-mesh sieve for NMR measurements.
The soil samples were then subjected to solid-state magic-angle rotation-NMR measurements (AVANCE II 300 MH) using a 7 mm CPMAS probe with an observed frequency of 100.5 MHz, an MAS rotation frequency of 5000 Hz, a contact time of 2 s, and a cycle delay time of 2.5 s. The external standard material for the chemical shift was hexamethyl benzene (HMB, methyl 17.33 mg kg−1). The spectra were quantified by subdividing them into the following chemical shift regions55: 0–45 ppm (alkyl), 45–60 ppm (N-alkyl and methoxyl), 60–110 ppm (O-alkyl), 110–140 ppm (aryl), 140–160 ppm (O-aryl), 160–185 ppm (carboxy), and 185–230 ppm (carbonyl) (Fig. 3a). We classified O-alkyl, O-aryl, and carboxy C as labile C and alkyl, N-alkyl/methoxyl, and aryl C were classified as recalcitrant C.
Soil microbial C metabolic profiles
The soil microbial C metabolic capacities were measured with BIOLOG 96-well Eco-Microplates (Biolog Inc., USA) using 31 different C sources and three replicates in each microplate. These C sources included carbohydrates, carboxylic acids, polymers, amino acids, amines, and phenolic acids (Table S2). Carbohydrates, amino acids, and carboxylic acids are generally considered labile C sources, amines and phenolic acid compounds are relatively resistant C sources, and polymers are recalcitrant C. The diverse nature of these C sources allowed us to identify differences in the capacity of microbes to degrade different C sources56. Soil microbes were extracted as follows: (1) Five grams of soil (dry weight equivalent) was incubated at 25 °C for 24 h, and 45 ml of sterile 0.85% (w/v) sodium chloride solution was added57. (2) At room temperature (25 °C), the mixture was shaken at 200 rpm for 30 min and allowed to stand for 15 min. (3) Subsequently, 0.1 ml of the supernatant was collected and diluted to 100 ml with sterile sodium chloride solution. (4) Soil suspensions were dispensed into each of the 93 wells (150 μl per well), and the plates were then incubated at 25 °C in the dark for 14 days. The optical density (OD, reflecting C utilization) of each well was read at 590 nm (color development) every 12 h. The normalized OD of different C sources was calculated as the OD of the well that contained the C source minus the OD of the well that contained sterile sodium chloride solution (control well). The normalized OD at a single time point (228 h) was used for the posterior analysis when it reached the asymptote.
DNA extraction, PCR amplification, and sequencing
DNA was extracted from all 90 soil samples. Briefly, well-mixed soil samples (0.6 g) were analyzed using the Power Soil DNA Isolation Kit (MoBio Laboratories, Inc., Carlsbad, CA, USA) following the manufacturer’s instructions. The quality of the DNA extracts was determined by spectrophotometry (OD-1000+, OneDrop Technologies, China). The DNA extracts were considered of sufficient quality if the ratio of OD260 to OD280 (optical density, OD) and the ratio of OD260 to OD230 were approximately 1.8. All eligible DNA samples were stored at −80 °C.
Taxonomic profiling of the soil bacterial and fungal communities was performed using an Illumina® HiSeq Benchtop Sequencer. PCR amplification was performed using an ABI GeneAmp® 9700 (ABI, Foster City, CA, USA) with a 20 μl reaction system containing 4 μl of 5× FastPfu Buffer, 0.8 μl of each primer (5 μM), 2 μl of 2.5 mM dNTPs, 2 μl of template DNA, and 0.4 μl of FastPfu Polymerase. For bacterial analysis, the forward the primer 515F (GTGCCAGCMGCCGCGG) and the reverse primer 907R (CCGTCAATTCMTTTRAGTTT) were used to amplify the bacteria-specific V4-V5 hypervariable region of the 16S rRNA gene58. For fungal analysis, the internal transcribed spacer 1 (ITS1) region of the ribosomal RNA gene was amplified with primers ITS1-1737F (GGAAGTAAAAGTCGTAACAAGG) and ITS2-2043R (GCTGCGTTCTTCATCGATGC)59. The PCR protocol for bacteria consisted of an initial predenaturation step of 95 °C for 2 min, 35 cycles of 20 s at 94 °C, 40 s at 55 °C and 1 min at 72 °C, and a final 10 min extension at 72 °C. The PCR protocol for fungi consisted of an initial predenaturation step of 95 °C for 3 min, 35 cycles of 30 s at 95 °C, 30 s at 59.3 °C, and 45 s at 72 °C and a final 10 min extension at 72 °C.
Each sample was independently amplified three times. Following amplification, 2 μl of each of the PCR products was checked by agarose gel (2.0%) electrophoresis, and all the PCR products from the same sample were then pooled together. The pooled mixture was purified using the Agencourt AMPure XP Kit (Beckman Coulter, CA, USA). The purified products were indexed in the 16S and ITS libraries. The quality of these libraries was assessed using Qubit@2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 systems. These pooled libraries (16S and ITS) were subsequently sequenced with an Illumina HiSeq 2500 Sequencer to generate 2 × 250 bp paired-end reads at the Center for Genetic & Genomic Analysis, Genesky Biotechnologies Inc., Shanghai, China.
The raw reads were quality filtered and merged as follows: (1) TrimGalore was used for truncation of the raw reads at any site with an average quality score <20, removal of reads contaminated by the adapter and further removal of reads with less than 100 bp; (2) the paired-end reads were merged to tags by Fast Length Adjustment of Short reads (FLASH, v1.2.11); (3) the removal of reads with ambiguous base (N base) and reads with more than 6 bp of homopolymer was performed with Mothur; (4) reads with low complexity were then removed to obtain clean reads for further bioinformatics analysis. The remaining unique reads were chimera checked by comparison with the gold.fa database (http://drive5.com/uchime/gold.fa) and clustered into operational taxonomic units (OTUs) using QIIME2 with a 97% similarity cutoff. The Ribosomal Database Project classifier performed the taxonomic assignment of OTUs with a minimal 70% confidence score60. For 16S data, taxonomic assignment was performed using the SILVA database (version 138, https://www.arb-silva.de/)61; for the ITS, the UNITE database was used (version 8.2, https://unite.ut.ee/)62.
To obtain the relationship between the number of detected species and the sequencing depth, we performed a linear fitting analysis between the number of OTUs and the number of reads (Fig. S11). The results showed that the number of species (OTUs) detected in the samples tended to level off once the number of bacterial and fungal reads was higher than 42,254 and 54,932, respectively (Fig. S11a, b). In this study, from all 90 samples, 14,414 OTUs were acquired for the bacterial community, and their read numbers ranged from 42,254 to 150,465; for the fungal community, a total of 9811 OTUs were obtained, and their read numbers ranged from 54,932 to 492,981. In the subsequent bioinformatics analysis, to minimize the impact of changes in the reading count among different samples, we rarefied all the samples based on the smallest read numbers (42,254 bacterial reads and 54,932 fungal reads per sample).
Microbial co-occurrence network construction
CoNet was used to generate interaction networks to determine the effects of climate change on the connections of bacteria and fungi in C-poor and C-rich soils. Zero-rich data were filtered before network construction. Briefly, the construction of a network graph was divided into four steps: basic configuration, permutation, bootstrapping, and restoration of the network from random files. Pairwise associations among OTUs were calculated using the Pearson, Spearman, Bray-Curtis, and Kullback-Leibler methods simultaneously. The initial top and bottom edge numbers were set to 1000. An edge- and measure-specific p-value was obtained as the area under the bootstrap distribution limited by the mean of the permutation distribution63. The edges were retained when supported by at least two correlation methods. The edges were discarded when the 95% confidence interval defined limits by the bootstrap distribution or the adjusted p-values were higher than 0.05. A final network was restored from the permutation and bootstrap files64. Graphs of the interaction network were built and visualized using Cytoscape (version 3.8.1). Each node in the network represents a bacteria/fungus, and the edges represent the correlation among different nodes.
Statistics and reproducibility
Analysis of changes in the SOC content and microbial C metabolism capacity
The relationships of SOC and DOC with latitude were evaluated using a fitting analysis. The fitting model corresponds to the minimum Akaike information criterion (AIC) value as the final model (16 samples under each OM gradient). In soils with the same OM, the relative abundance of C molecules in the translocated sites was divided by the relative abundance of C molecules in the in situ HL site to evaluate the stability of different C components (i.e., alkyl, O-alkyl, N-alkyl/methoxyl, aryl, O-aryl, carboxy C) under elevated temperatures. A total of 45 ratios were obtained for each OM content. A ratio greater than 1 represents increased stability, and a ratio less than 1 represents decreased stability. Similarly, in C-poor (OM ≤ 5%) and C-rich (OM > 5%) soils, changes in the C metabolic capacity of microbes under elevated temperatures were characterized using the ratio of the OD of microbes measured in the translocated soils to the OD of microbes in the in situ HL soil. A ratio greater than 1 indicates that translocation warming increases the C metabolism of microbes.
Mantel and partial Mantel analysis
A previous study showed that partial Mantel analysis is a robust method for evaluating the relationship among three variables65. This approach can control the z-axis and assess only the relationship between the x– and y-axes, avoiding the interaction between the z– and x-axes on the y-axis. In this study, Mantel analysis was employed to assess the relationships between the stability of the bacterial and fungal communities and C metabolic capacity. Stability refers primarily to the ability of the microbial community to resist translocation warming66. A higher similarity between the microbial communities in translocated soil compared with that in the in situ HL area indicates that the community is more resistant to translocation-related warming and that the microbial community is more stable.
Calculation of the microbial β-diversity
Bray-Curtis and Euclidean dissimilarity metrics were calculated to estimate the bacterial and fungal taxonomic dissimilarity (β-diversity) and environmental dissimilarity (e.g., latitude, MAT, and MAP), respectively, using the vegan package (version 2.5–6) in the R statistical program (version 4.0.2, https://www.r-project.org/)67. Corresponding to the 45 C metabolism ratios in soils with the same OM content, the β-diversity values of bacteria and fungi were selected to analyze the relationship between the community similarity (1-β-diversity) of bacteria and fungi and changes in microbial C metabolism.
Impact of the SOM content and climate change on changes in microbial communities
The distribution patterns of the bacterial and fungal communities under different SOM gradients and climatic regimes were determined through nonmetric multidimensional scaling (NMDS)68. To quantitatively compare the effects of the SOM gradient and climatic regimes on the bacterial and fungal community composition, three nonparametric multivariate statistical analyses were used in this study: nonparametric multivariate analysis of variance (Adonis), analysis of similarity (ANOSIM), and multiple response permutation procedure (MRPP)69. The linear fit between environmental dissimilarity and microbial β-diversity was analyzed using the lm function in R. A significant difference in the bacterial and fungal β-diversity among different SOM contents was evaluated by Student’s paired t-test using the ggpubr (version 0.4.0) package70. RDA was performed to analyze the relationships of bacterial and fungal communities with various environmental factors (soil geochemical attributes and climatic conditions, such as MAP and MAT). In parallel, the Monte Carlo permutation test (999 permutations) was employed to determine whether the explanation of the microbial distribution by individual factors (e.g., pH, SOC, and TN) was significant71.
Construction of the structural equation model and random forest model
A SEM was fitted to illustrate the direct or indirect effects of soil properties (e.g., pH, moisture, ammonia, and nitrate nitrogen), climate change (e.g., MAT and MAP), and bacterial and fungal β-diversity on soil C metabolic capacity72. Based on the Euclidean method, the changes in soil properties and climatic conditions of five translocated sites compared with those in the in situ HL site were calculated. A total of 45 ratios were obtained for each OM content. Corresponding to the 45 ratios in soils with the same OM content, the β-diversity values of bacteria and fungi were selected. The model construction process was mainly divided into three steps. In brief, these steps include the establishment of an a priori model, data normality detection, and an overall goodness-of-fit test. The prior model was constructed based on a literature review and our knowledge. For the variables that did not conform to the normal distribution, we performed logarithmic transformation. Here, we used the χ2 test (the model was assumed to exhibit a good fit if p > 0.05), the goodness-of-fit index (GFI; the model was assumed to show a good fit if GFI > 0.9), the root mean square error of approximation (RMSEA; the model was assumed to exhibit a good fit if RMSEA < 0.05 and p > 0.05)73 and the Bollen-Stine bootstrap test (the model was assumed to show a good fit if the bootstrap p > 0.10) to test the overall goodness of fit of the SEM. All SEM analyses were conducted using IBM® SPSS® Amos 21.0 (AMOS, IBM, USA). Additionally, the importance of the metabolic capacity of different types of C on labile and recalcitrant C was assessed by random forest models using the randomForest package (version 4.6-14) in R74, and the model significance and amount of interpretation were evaluated using the rfUtilities package (version 2.1–5)75.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
