Site description and experimental design
The study was conducted in a typical semiarid steppe ecosystem (42°06′44′′N, 115°29′10′′E; 1450 m above sea level), which is located in Baochang, Inner Mongolia, China. This region is characterized by a dry, monsoon-influenced continental climate with a mean annual temperature of 1.6°C and annual precipitation of 400 mm, with approximately 86% of the precipitation events occurring during the growing season (May-September). In this region, the soil within the top 20 cm is designated chestnut soil (Chinese Soil Taxonomy Research Group, ISSAS, 2001), with 80.50 ± 2.58% sand, 15.3 ± 2.28% silt, and 4.17 ± 0.40% clay. The soil total C and N contents were 23.06 ± 1.98 g kg− 1 and 2.44 ± 0.14 g kg− 1, respectively; the soil pH was 8.14 ± 0.17. The annual natural N deposition in this area is approximately 3.43 g N m− 2 yr− 1 (He et al. 2016). The plant community in the study area is dominated by Stipa krylovii, Artemisia frigida, and Leymus chinensis.
In 2011, a field experiment with various N addition levels was established to examine the potential effects of N deposition on grassland ecological processes. It comprises 15 plots (each 3 m in width × 6 m in length), consisting of five different levels of N addition (0, 2, 5, 10, and 25 g N m− 2 yr− 1), denoted as N0, N2, N5, N10, and N25, respectively. Adjacent plots were separated by a 1-m buffer zone between plots to avoid N fertilizer transfer. Each plot was randomly assigned one of the five N levels, resulting in three replications for each level. N additions were applied two times in early June and July of each year to the surface of each plot as dilute sodium nitrate (NaNO3) solution using a backpack sprayer. The control plots without N fertilizer (i.e., N0) received equal amounts of water. The water applied to each plot was kept consistent, amounting to approximately 2 mm of local precipitation. Following the estimated annual N deposition rate in northern China (0.042 g N m− 2 yr− 1), N2, N5, and N10 simulate N deposition levels for the next 50, 100, and 200 years, respectively. In turn, N25 simulates the effects of N fertilization in agricultural fields (He et al. 2016).
Field sampling and measurement
In July 2021, the plant community structure (including plant species richness, abundance, height, and coverage) was surveyed in a randomly selected 0.5 m × 0.5 m quadrat within each plot. Meanwhile, another 0.5 m × 0.5 m quadrat was selected in the buffer zone of each plot, and all aboveground plant individuals in the quadrat were harvested after investigating the community structure. We dried the harvested aboveground live (green) material of individual plant species at 65°C for 48 h, weighed them, and then converted them to the plot-level aboveground biomass (AGB, g m− 2) of our corresponding studied plots. Plant diversity (as measured by the Shannon‒Wiener index) for each plot was calculated based on the field investigation (Risser and Rice 1971).
Three randomly located soil cores (8 cm diameter) were taken at depths of 0–10, 10–20, and 20–40 cm for each plot, and subsamples from the same soil layer were composited. In total, 45 soil samples (15 plots × 3 soil layers) were taken in this study. After the removal of stones and visible roots, the soils were sieved with a 2-mm mesh and thoroughly mixed. Each soil sample was then divided into three equal parts. One part was designated fresh soil samples and stored at 4°C for incubation experiments. The second part was air-dried in preparation for subsequent soil property analysis (Table 1), while the third part was freeze-dried and stored at -80°C for assessing soil microbial diversity. All live roots extracted from soil samples were carefully rinsed in deionized water. Fine roots (diameter < 2 mm) were manually picked, dried at 65°C for 48 h, and then weighed to determine the plot-level belowground biomass (BGB, g m− 2). We also measured fine root properties after weighing (Table 1). Further details on the methods of soil and root physiochemical properties determinations are provided in the Text S1.
Microbial community diversity and C-degrading functional groups
To profile soil microbial communities, we amplified the V3-V4 hypervariable region of the 16S rRNA gene for soil bacteria and the internal transcribed spacer region for soil fungi, and then sequenced the PCR products on the Illumina HiSeq 2500 platform (Illumina). To ensure a standardized sampling effort, we rarefied all samples to achieve an even number of sequences per sample for bacteria (52,293) and fungi (38,704). These sequences were clustered into 15,328 bacterial operational taxonomic units (OTUs) and 3,124 fungal OTUs at a 97% similarity level (Wang et al. 2022). We finally calculated the richness (total number of unique OTUs), diversity (as measured by the Shannon‒Wiener index), and evenness (using Pielou’s evenness index) for both soil bacteria and fungi according to the number of OTUs obtained. The sequencing process of this study was performed at Magi Gene Technology (Guangzhou, China).
Microbial community diversity and C-degrading functional groups
To profile soil microbial communities, we amplified the V3-V4 hypervariable region of the 16S rRNA gene for soil bacteria and the internal transcribed spacer region for soil fungi, and then sequenced the PCR products on the Illumina HiSeq 2500 platform (Illumina). To ensure a standardized sampling effort, we rarefied all samples to achieve an even number of sequences per sample for bacteria (52,293) and fungi (38,704). These sequences were clustered into 15,328 bacterial operational taxonomic units (OTUs) and 3,124 fungal OTUs at a 97% similarity level (Wang et al. 2022). We finally calculated the richness (total number of unique OTUs), diversity (as measured by the Shannon‒Wiener index), and evenness (using Pielou’s evenness index) for both soil bacteria and fungi according to the number of OTUs obtained. The sequencing process of this study was performed at Magi Gene Technology (Guangzhou, China).
Table 1
Microbial respiration and relevant abiotic/biotic factors (soil properties, plant attributes, microbial community structure, and C-degrading functional groups) and their abbreviations.
Variables | Abbreviation (unit) |
Microbial respiration | MR (µg C g− 1 soil d− 1) |
Soil properties |
soil water holding capacity | WHC (g g− 1 soil) |
soil gravimetric water content | SWC (g g− 1 soil) |
soil pH | pH |
soil clay content | sand (%) |
soil silt content | silt (%) |
soil sand content | slay (%) |
soil available phosphorus | AP (mg kg− 1 soil) |
soil total carbon content | soil C (%) |
soil total nitrogen content | soil N (%) |
soil carbon and nitrogen ratio | soil C:N |
soil ammonium nitrogen concentration | NH4+-N (mg kg− 1 soil) |
nitrate nitrogen concentration | NO3−-N (mg kg− 1 soil) |
total dissolved carbon concentration | DC (mg C L− 1) |
total dissolved nitrogen concentration | DN (mg C L− 1) |
total dissolved carbon and nitrogen ratio | DC:DN |
total dissolved organic carbon concentration | DOC (mg C L− 1) |
Plant attributes |
fine root total carbon content | root C (%) |
fine root total nitrogen content | root N (%) |
fine root carbon and nitrogen ratio | root C:N |
total amount of carbon in the fine root | FRBC (g m− 2) |
total amount of nitrogen in the fine root | FRBN (g m− 2) |
Microbial community structure |
soil fungal to bacterial ratio | F:B |
Microbial C-degrading functional groups |
decreased carbon-degrading bacteria | DCB |
increased carbon-degrading bacteria | ICB |
decreased carbon-degrading fungi | DCF |
increased carbon-degrading fungi | ICF |
Microbial community diversity and C-degrading functional groups
To profile soil microbial communities, we amplified the V3-V4 hypervariable region of the 16S rRNA gene for soil bacteria and the internal transcribed spacer region for soil fungi, and then sequenced the PCR products on the Illumina HiSeq 2500 platform (Illumina). To ensure a standardized sampling effort, we rarefied all samples to achieve an even number of sequences per sample for bacteria (52,293) and fungi (38,704). These sequences were clustered into 15,328 bacterial operational taxonomic units (OTUs) and 3,124 fungal OTUs at a 97% similarity level (Wang et al. 2022). We finally calculated the richness (total number of unique OTUs), diversity (as measured by the Shannon‒Wiener index), and evenness (using Pielou’s evenness index) for both soil bacteria and fungi according to the number of OTUs obtained. The sequencing process of this study was performed at Magi Gene Technology (Guangzhou, China).
Based on the sequencing results, we identified six bacterial and four fungal phyla reported to be closely involved in C degradation (Fig. S1) (Pankratov 2012; Tláskal et al. 2017). We calculated the relative abundance of OTUs for these phyla in each soil sample, which is the ratio of the number of OTUs for these phyla to the total number of OTUs in the sample. From linear regression of MR against the relative abundance of these C-degrading phyla (Fig. S2), we classified them into four C-degrading functional groups: decreased C-degrading bacteria (DCB), increased C-degrading bacteria (ICB), decreased C-degrading fungi (DCF), and increased C-degrading fungi (ICF). For each group, the relative abundance of all involved phyla was summed to give the group-level relative abundance that was used in the subsequent analysis.
Soil microbial respiration
We performed a laboratory incubation experiment under varying temperature conditions (i.e., 10-15-20-25-30-25-20-15-10°C) to measure the MR for 45 soil samples by a 32-channel MR automatic measurement system (Zhang et al. 2022). This system incorporated a LI-840A CO2/H2O infrared gas analyzer (IRGA; LI-COR Inc., Lincoln, NE, USA), an air pump (LI-COR Inc.), a laptop (Surface Pro 3, Microsoft, USA), a gas exchange system and an artificial weather box (IGS100, Thermo). The MR values (µg C g− 1 soil d− 1) were calculated as the mean daily microbial decomposition rates of soil C across different temperatures.
Statistical analysis
We used two-way analysis of variance (ANOVA) to assess the individual and interaction effects of soil depth and N addition on the MR. This analysis was also applied for soil properties, microbial community structure, and C-degrading functional groups. The Mantel test was then performed to examine the linkage between soil properties, plant attributes and microbial community structure.
We finally constructed structural equation models (SEMs) to quantitatively discern the direct and indirect relationships between N addition and MR through various factors across different soil layers (Grace 2009). To make the SEM more concise and intuitive, we used random forest analysis to identify the variables with the greatest explanatory power on MR (Breiman 2001; Yuan et al. 2023). Random forest is a versatile machine-learning algorithm that is used for making predictions and involving both linear and non-linear main effects, and interaction effects (Strobl et al. 2007). It is relatively insensitive to multicollinearity and overfitting; we incorporated all the previously mentioned variables in this step. The importance of the selected variables on MR was quantified using the percentage increase in mean square error (%lnMSE). We chose the top ten important variables that most decreased the mean square error (MSE) over random permutations (Fig. S3). The ten selected variables were used as starting variables for the subsequent analyses. The random forest algorithm was performed using the R package “randomForest” (Liaw and Wiener 2002).
To avoid overfitting caused by high collinearity among the selected variables, we excluded one variable from any pair of predictor variables with a correlation coefficient (Pearson’s r) larger than 0.70 (Fig. S4) (Dormann et al. 2013; Luo et al. 2019; Zhang et al. 2019). Only the retained variables were then included in the initial SEM. Given that no variable was removed for the medium soil layer (all Pearson’s r values were less than 0.7, Fig. S4b), we proceeded with a principal component analysis (PCA) to reduce the dimensions for soil properties (i.e., WHC, AP, soil N, DN, and silt), simplifying our analyses and facilitating interpretations (Yuan et al. 2023). The first two components of the PCA result (soil_PC1 and soil_PC2) captured 72% of the total variance and were thereafter used to characterize the soil properties of the medium soil layer for subsequent SEM analyses (Fig. S5).
The results from correlation analyses revealed that NH4+-N, fungal diversity, AGB, ICF, soil C:N, and WHC were introduced into the final SEM in the shallow soil layer. Comparatively, ICF, root C, fungal diversity, ICB, soil_PC1, and soil_PC2 were considered in the medium soil layer. In the deep soil layer, we constructed the SEM by including soil pH, sand, fungal diversity, plant diversity, ICB, and bacterial richness.
Our conceptual SEMs were constructed based on the assumption that N addition influences MR indirectly by regulating soil–plant–microbe interactions (Fig. 1). To examine whether the direction and magnitude of the N addition effect on MR changed with soil depth, we performed separate SEM analyses for each soil layer. The fits of SEMs were evaluated using chi-squared (χ2) tests and associated p values (Grace 2009). SEM analyses were carried out using IBM SPSS Amos 26 (Amos Development Corporation). All statistical analyses, except for the SEM, were performed by R version 4.2.2 (R Core Team, 2021).