Introduction

Global change can greatly affect plant productivity and C allocation between above- and below-ground components, consequently driving uncertain ecosystem C cycling in the terrestrial biosphere1,2. Variations in the quantity and quality of plant litter and roots under global change, the primary sources of soil C pools, have profound influences on soil C dynamics1,2. Considerable efforts have been focused on estimating the impact of litter manipulations on the formation, stabilization, and turnover rates of soil organic C (SOC) using Detritus Input and Removal Treatment, and the results suggest that the increase in SOC under litter addition (LA) may be lower than the decrease in SOC under litter removal (LR)2,3. Moreover, a potentially stronger contribution of below-ground C inputs compared to aboveground C inputs to SOC sequestration has been recognized, primarily because root-derived C can be more efficiently retained in soils4,5. Therefore, better understanding on the effects of above- and below-ground C inputs on soil biogeochemical processes will help to reveal the feedback mechanisms of terrestrial biosphere in response to global change.

Microbes are crucial for soil C and nutrient cycling in terrestrial ecosystems6,7. Given the strong dependences of microbial growth and activity on C substrates, changes in soil C sources can lead to alterations in microbial biomass, diversity, composition, and functions1,8. Plant litter may have equal or greater contributions than roots to soil biota and ecological processes9,10,11, but the effects of litter manipulations vary with time and ecosystem types12,13. By contrast, growing evidence shows that plant roots have a greater impact on soil C budget than litter8,14 by providing the primary source of C for soil microbes that rely closely on labile C in the rhizosphere15. Consequently, excluding soil C inputs from plant roots may exert stronger negative influences on soil microbial community than LR8,14. However, on the one hand, some short-term field experiments could not fully reflect the different responses of soil microbial communities to plant litter and root C inputs due to the temporal lag of microbial responses to LA16,17,18. On the other hand, most previous studies have focused on the effects of C inputs on microbial biomass and their contributions to SOC decomposition, the relative roles of above- and below-ground C inputs in regulating microbial diversity, community structure, and co-occurrence interactions remain largely unclear.

Different soil microbial groups may have distinct preferences for C sources derived from plant litter or roots. Enhanced C inputs from aboveground litter can promote fungal biomass and shift the soil microbial community from bacteria dominance to fungal dominance, whereas LR can have positive, negative, or neutral effects on soil microbial composition due to the complex soil properties2,19,20. Root exclusion has a significant impact on fungal communities, especially on symbiotic fungi and rhizosphere-associated microbes1,14. In addition, protists, representing most of the eukaryotic diversity in the soil, occupy critical niches in soil food webs by preying on bacteria and fungi, leading to changes in microbial biomass, diversity, and community structure, which may have essential functional consequences21. Previous studies have indicated that C sources from both litter composition and root secretion might vary the assembly process of soil protists, and then reshape bacterial and fungal communities through the top-down effect21,22,23. However, due to the higher trophic level and a wider range of ecological niches, soil protists are likely to show stronger redundancy features than bacteria and fungi facing environmental changes23,24,25. Therefore, considering the complex interactions between various microbial groups7,26, multitrophic level interactions under changing C inputs from plant litter and roots may potentially affect soil functions.

Ecosystem multifunctionality (EMF), which refers to the simultaneous provision of multiple ecosystem functions, can be used to assess ecological services associated with soil microbial diversity under changing C input pathways6,7,26,27. Due to their differential substrate utilization strategies, bacteria and fungi may have contrasting associations with EMF due to the niche complementarity and species interactions28,29. In terrestrial ecosystems, bacterial communities may dominate in regulating soil functions, because they possess greater metabolic versatility and are better adapted to relatively barren conditions4,29,30,31. For example, soil bacterial communities have been documented to dominantly regulate the decomposition and formation of soil organic matter, supply soluble nutrients promote plant growth, and reduce phytopathogens32,33. In contrast, a growing body of research suggests the more critical role of soil fungal diversity in driving EMF in drylands, grasslands, and subtropical forests7,34,35,36, which can be attributed to their stronger resistance to environmental changes, greater sensitivity to SOC dynamics, and higher C use efficiency37,38. Such superiorities enable soil fungi to exert greater impacts on maintaining soil functional stability, accelerating the migration and transformation of environmental pollutants, and thus increasing ecosystem productivity1,7. However, our knowledge of how protists affect microbial associations and, subsequently, EMF under changing C inputs is extremely limited.

Grasslands play a crucial role in regulating complex ecosystem functions due to their extensive coverage (over one-third of Earth’s land surface) and huge soil C stocks29,39. As part of a field experiment with 13-year soil C input manipulations in a semi-arid grassland on the Mongolian Plateau, this study was conducted to investigate the responses of soil microbiomes, including bacterial, fungal, and protistan diversity, community composition, and co-occurrences networks, to different C input ways, and their contributions to EMF. Three hypotheses were proposed: (1) LA enhances microbial diversity due to extra C substrate supply whereas LR reduces it due to decreased C substrate inputs, (2) The effect of plant removal (PR) on microbial community is greater than that of litter manipulations because root-derived C is more efficient for microbial growth, and (3) Soil protistan community is not sensitive to changes in C inputs and thus has a minor impact on EMF through bottom-up and top-down effects on soil trophic cascades because predatory protists tend to be generalists.

Results

Soil properties and microbial activity

Soil moisture, temperature, NO3−-N, dissolved organic C (DOC), and N (DON) exhibited significant fluctuations between the early and middle growing seasons (Supplementary Fig. 1; all P < 0.001, Repeated-measures ANOVA, Supplementary Table 1). Averaged over the two sampling stages, PR increased soil temperature (0.35 °C, absolute change, P < 0.01) and NO3−-N (39.0%, P < 0.001), reduced NH4+-N (30.2%, P < 0.01), DOC (40.4%, P < 0.001), DON (39.2%, P < 0.001), soil organic C (SOC; 0.16%, absolute change, P < 0.001), TC (28.4%, P < 0.001) and TN (27.5%, P < 0.001), but had no effect on soil moisture (P = 0.794) or C/N ratio (P = 0.219). LA enhanced DOC (26.1%, P < 0.001), DON (30.0%, P < 0.01), and TC (8.97%, P < 0.05), but decreased soil temperature by 0.39 °C (P < 0.05). LR reduced soil moisture only (0.01%, absolute change, P < 0.01). Neither LA nor LR affected SOC, NH4+-N, NO3−-N, TN, or C/N ratio (all P > 0.05). PR interacted with LA to impact DOC (P < 0.05), and with LR to influence soil moisture (P < 0.01), and SOC (P < 0.05).

Across the six treatments, potential nitrification rate (PNR), nitrogen mineralization rate (Nmin), nitrogenase activity (NA), and soil basal respiration (SBR) exhibited significant variations between the early and middle growing season (Supplementary Fig. 2; all P < 0.01, Repeated-measures ANOVA, Supplementary Table 2). PR significantly reduced PNR, Nmin, NA, and SBR by 27.2%, 50.3%, 50.4%, and 68.7% averaged over the two stages, respectively (all P < 0.001). LA increased SBR by 19.7% (P < 0.05), but had no effect on other microbial variables. There was no interaction between PR and LA on PNR, Nmin, NA, and SBR (all P > 0.05). Neither LR nor its interaction with PR affected microbial activities (all P > 0.05).

Microbial diversity and composition

PR significantly decreased the α-diversity of both bacterial and fungal communities by 10.6% and 15.5%, respectively (Fig. 1a, b; both P < 0.001, Repeated-measures ANOVA, Table 1; Supplementary Fig. 3). LA did not affect bacterial and fungal α-diversity, but LR reduced the α-diversity of bacteria and fungi by 4.15% (P < 0.001) and 20.0% (P < 0.01), respectively. PR interacted with LA and LR to affect bacterial α-diversity, whereas fungal α-diversity was influenced only by the interaction of PR and LR. Protistan α-diversity remained unchanged under PR, LA, or LR (Fig. 1c; all P > 0.05). In addition, the β-diversity of all the bacterial, fungal, and protistan communities exhibited significant variations between the early and middle growing season (Fig. 1d–f; all P < 0.05, PERMANOVA, Supplementary Table 3), and were significantly influenced by PR (P < 0.01), LA (P < 0.01) and LR (P < 0.05), with PR having stronger effects on bacterial β-diversity than LA and removal.

Fig. 1: Microbial diversity under the six treatments during two sampling stages.
figure 1

Box plots of Chao1 richness of soil a bacterial, b fungal, and c protistan communities. Nonmetric multidimensional scaling (NMDS) analyses based on Bray-Cutis distance metrics of soil d bacterial, e fungal, and f protistan β-diversity, and the 95% confidence ellipses are shown around the samples. C control, LA litter addition, LR litter removal, PR plant removal, LAPR litter addition plus plant removal, LRPR litter removal plus plant removal. With plant: C, LA, and LR; without plant: PR, LAPR, and LRPR.

Full size image
Table 1 Results (P values) of two-way repeated-measures ANOVAs on the effects of plant removal (PR), liter manipulations (L) (litter addition (LA) and litter removal (LR)), sampling stage (S), and their interactions on the α-diversity of bacteria, fungi, and protists, as well as functions of microbial activity (EMFMA), carbon storage (EMFC), nutrient provisioning (EMFN), and averaging ecosystem multifunctionality (EMF)
Full size table

Across the six treatments, bacterial community was mainly composed of Actinobacteria (42.0%), Proteobacteria (22.5%), and Acidobacteria (12.1%). Plant removal significantly increased the relative abundance of Actinobacteria and Chloroflexi by 2.15% and 20.2%, respectively, but decreased that of Proteobacteria by 24.5%. LA elevated the relative abundance of Proteobacteria by 4.65%, whereas LR reduced it by 7.24% (Fig. 2a, P < 0.05, one-way ANOVA). Ascomycota (71.5%) was the most dominant phylum of fungi, but remained largely unchanged under all the treatments. PR stimulated the relative abundance of Mortierellomycota by 3.5-fold (Fig. 2b, P < 0.05). Rhizaria (36.5%), Alveolata (24.1%), Archaeplastida (14.4%), and Stramenopiles (12.3%) were the main protistan taxa (Fig. 2c). Archaeplastida were stimulated whereas Alveolata was reduced under the PR treatment (Fig. 2c, P < 0.05). At the genus level, LA did not change, whereas LR enriched 11 genera mostly belonging to fungi, but depleted three protistan genera (Fig. 2d, e). PR enriched 221 genera primarily belonging to bacteria and fungi, but depleted a total number of 183 genera, including 15 protistan genera only (Fig. 2f, P < 0.05).

Fig. 2: Microbial composition and differential genera among treatments.
figure 2

The relative abundance of soil a bacterial, b fungal, and c protistan communities under the six treatments. Volcano plots of the differential genera under the main effects of d litter addition, e litter removal, and f plant removal. The blue dots represent the up-regulated genera and the pink dots represent the down-regulated genera under each treatment. The square, circle, and triangle shapes represent genera belonging to the bacterial, fungal, and protistan communities.

Full size image

Co-occurrence networks

Co-occurrence patterns of soil microbiomes between bacterial-fungal-protistan interkingdom were assessed by combining data from the early and middle growing seasons. Compared to the control and PR treatments, all the networks with LA and LR had higher proportions of negative correlations, higher average path distance, and lower average clustering coefficient (Supplementary Fig. 4). There were significant differences between the networks with and without plants. PR decreased the average degree (Fig. 3a; 6.96 vs. 2.61, Supplementary Table 4), path distance (4.28 vs. 2.28), the proportions of fungi nodes (22.6 vs. 4.1), but enhanced the proportions of bacteria nodes (67.1 vs. 87.8). Specifically, bacteria taxa had higher network connectivity (i.e., network degree) than fungal and protistan taxa (Fig. 3b). The network with plants had more stable stability than that without plants (R2 = 0.26 vs. 0.12; Fig. 3d). The number of ‘hub nodes’ (nodes with high values of degree (>30) and closeness centrality (>0.3)) was also higher in the network with than without plants, and all these hub nodes belonged to bacteria (Fig. 3e).

Fig. 3: Microbial co-occurrence networks and topological features.
figure 3

a Co-occurrence patterns of soil bacterial-fungal-protistan interkingdom with and without plant. The topological features of b average degree, c path length, and d robustness between these two co-occurrence networks. e The comparison of node-level topological features (degree and closeness centrality) between bacterial, fungal, and protistan taxa in the networks with and without plant.

Full size image

EMF and influencing factors

Averaged over the two stages, PR reduced the multifunctionality of microbial activity, C storage, nutrient provisioning, and thus the averaging EMF (Fig. 4a–d; all P < 0.001, Repeated-measures ANOVA, Table 1). LA enhanced C storage function (P < 0.01) and averaging EMF (P < 0.05), whereas LR did not affect all the induvial functions and EMF (P > 0.05).

Fig. 4: Different components of ecosystem multifunctionality (EMF).
figure 4

a EMF of microbial activity (EMFMA) including potential nitrification rate, nitrogen mineralization rate, nitrogenase activity, and soil basal respiration. b EMF of C storage (EMFC) including soil-dissolved organic carbon, organic carbon, total carbon, heterotrophic respiration, and soil respiration. c EMF of nutrient provisioning (EMFN) including soil NH4+-N, NO3−-N, dissolved organic nitrogen, total nitrogen, and C/N ratio. d Averaging EMF including all the above variables under the six treatments.

Full size image

Soil moisture and DOC were the most important factors influencing the α-diversity of soil bacteria and fungi, whereas soil protistan diversity was not affected by the changes in environmental variables (Supplementary Fig. 5). Mantel test and heatmap analysis showed that bacterial communities had stronger correlations with soil properties and microbial activities, especially with DOC, TC, and TN, than fungal and protistan communities. By contrast, EMF (including microbial activity, soil carbon storage, and nutrient provisioning) was more closely related to fungal diversity than to bacterial and protistan diversity (Supplementary Fig. 6). A structural equation model (SEM) was further used to quantify the effects of PR and litter manipulations on soil microbiomes and EMF. The results showed that EMF was not affected by bacterial and protistan diversity under PR or litter manipulations, but was suppressed by reducing fungal diversity under PR (Fig. 5a). LR, rather than LA, interacted with PR to affect bacterial and fungal diversity and then regulated EMF (Fig. 5b).

Fig. 5: The influence pathways of plant removal and litter manipulations on EMF.
figure 5

Structural equation model (SEM) analyses of the effects of a plant removal and litter addition, and b plant removal and litter removal on soil microbial diversity and EMF. The final model fits well based on the chi-square and RMSEA tests (a χ2= 0.230, P = 0.973, RMSEA = 0.000, df = 3; b χ2 = 0.758, P = 0.859, RMSEA = 0.000, df = 3). The black dashed arrow represents a non-significant pathway. Orange and black solid arrows indicate significant positive and negative pathways, respectively. Numbers adjacent to the arrows are standardized path coefficients with different significance levels: ^P < 0.100, *P < 0.050, **P < 0.010, and ***P < 0.001. R2 value for each dependent variable refers to the proportion of variance explained by the model.

Full size image

Discussion

Our results do not fully support Hypothesis 1, that 13-year LA would increase soil microbial diversity, whereas LR would reduce it. In this study, LR decreased bacterial and fungal diversity but did not affect protistan diversity, whereas LA had no effect on any type of microbial diversity. These findings are not consistent with those of previous research, which has shown that LA typically increases soil microbial diversity by providing additional resources, whereas LR has the opposite effect20,40,41. However, our study area is a semi-arid grassland that is largely limited by water, and the changes in active C and N levels resulted from our litter manipulations may not have been significant enough to cause substantial variation in microbial diversity20,42. Our data do support the proposition that reduced soil water availability resulted from LR can decrease bacterial and fungal diversity, because water stress caused by the removal of aboveground litter can exacerbate the negative effect of decreased C substrate supply on microbial diversity20. Although it has well been documented that soil water availability is the primary factor influencing protistan communities at large spatial scales43,44, our observations suggest that decreased soil moisture, bacterial and fungal diversity under LR do not significantly alter protistan diversity in our temperate steppe study site. Nevertheless, it is possible that specific abiotic and biotic factors have differential impacts on protistan lineages, resulting in shifts in protistan composition or microbial associations.

Our study revealed significant variations of soil microbial β-diversity and compositions under LA and removal. These findings agree with those of previous studies in both laboratory and field experiments that litter manipulations have a significant impact on bacterial and fungal community structure2,27,45. This is likely due to the expansion of fast-growing bacteria and fungi, which can respond quickly to changes in labile C substrates45,46. The predominant phylum Proteobacteria increased in abundance under LA and decreased under LR because many members of this group primarily utilize complex C substrates derived from plant litter45,46,47. In contrast, fungal abundance at the phylum level remains stable under both LA and removal, which is consistent with results of previous studies that fungal community has great resistance to environmental disturbances and thus maintains the stability of soil microbiomes2,48. Moreover, several possible reasons may explain the discrepant responses between bacteria and fungi to litter manipulations. First, fresh litter in grasslands contains less recalcitrant components which can be more easily used by dominant bacteria during the decomposition process, hence leading to the shifts in bacteria community27. Second, the enhanced soil C/N ratio caused by litter input (with relative higher C/N ratio of approximately 61.2 globally) may induce a lower N availability and favor fungal growth49,50, neither LA nor removal affected soil C/N ratio, thus exerting neutral effects on fungal community. Third, the loss of soil labile C under LR has less impacts on fungi because they can produce oxidative enzymes more efficiently than bacteria to decompose recalcitrant C to maintain their growth19. Therefore, bacteria may show more sensitive responses than fungi to litter manipulations in semi-arid grasslands.

Protists, the most taxonomically diverse eukaryotes, play critical roles in soil food webs21,43. Most soil protists are phagotrophic and major consumers of other soil microorganisms, likely resulting in important functional consequences on litter decomposition processes43,44. Despite the fact that the α-diversity of the protistan community does not change under litter manipulations, both the β-diversity and the community composition of soil protists have varied under litter manipulations. In specific, the abundance of phagotrophic Amoebozoa, which generally participates in the SOC decomposition and nutrient mineralization processes23,51, and is closely related to lower-tropical microbial diversity and ecological functions, has been enriched by LA, suggesting that litter manipulations may affect the niches of specific protists in soil microbial communities23,51,52.

Given that root-derived C, rather than litter-derived C, is the primary source of C for soil microbes, especially those associated with the rhizosphere in most terrestrial ecosystems1,15, PR may have a greater effect on soil microbial communities than litter manipulations. In our study, we observed more dramatic decreases in bacterial and fungal α-diversity, as well as variations in microbial β-diversity and compositions under PR than under litter manipulations. This is likely due to the strong C limitation and large decline in soil nutrient pools under decreased fresh C substrates (e.g., exudates) from living roots2. Plant root exudates generally contain amino acids, sugars, and solubilize phosphate, which can not only supply the C sources for the growth of soil microorganisms, but also promote C and nutrient cycles that are regulated by microorganisms in the soil53. Consequently, soil microbes may be more restrained by the reduced availability of soil C and nutrients under PR compared to LR. Although there is a lack of C component estimation of root secretion in this study, the reduced soil DOC, SOC and TC can partially support our above speculations. However, due to their ability to efficiently use litter-derived C sources, Actinobacteria can live in soils with limited C sources and maintain a relative higher abundance under PR37. PR can also directly decrease the biomass and diversity of symbiotic fungi (e.g., arbuscular mycorrhizal and ectomycorrhizal fungi) associated with plant roots1,14, which is supported by the sharper decline in fungal diversity under PR than under LR in our experiment. Although PR may alleviate the adverse effects of some exudates on soil environment and microorganisms, for example, the organic acids secreted by the roots may reduce soil pH, driving the abundance changes in Burkholderia species which dominantly metabolize citrate and oxalate; and the bacterial and fungal diversity are likely to be inhibited by the isoflavonoids and fungal cell wall-degrading enzymes53, our findings suggest that PR does have much stronger negative effects on soil bacterial and fungal community than litter manipulations in this semi-arid grassland, which is consistent with our Hypothesis 2.

The protistan richness remains unchanged after PR, suggesting that protistan diversity might be not affected by changing soil C inputs or their influence on soil microclimate, C and nutrient availability, as well as the bottom-up effects of soil bacterial and fungal diversity. All the driving forces on protists can be explained by the differences in soil moisture because higher water availability facilitates protists’ feeding on bacteria and fungi in water-filled soil pores43,44. Correspondingly, soil moisture does not affect protistan diversity. However, PR reduces the abundance of the phagotrophic taxa Alveolata, but increases the abundance of the phototrophic taxa Archaeplastida, suggesting that the availability of C resource and light due to PR had substantial effects on protistan communities. In contrast, our results suggest that PR, rather than litter manipulations, predominantly determines the complexity and stability of soil microbial networks, which can be supported by the decreased topological features. Therefore, PR exerts strong effects on soil microbiomes in this semi-arid grassland.

Several studies from different ecosystems2,44,54 partly supported the observations in this study that some soil functions are enhanced by LA but remain unchanged under LR, whereas PR reduces not only microbial metabolism but also the ecosystem functions of soil C storage and nutrient provisioning. The input of litter C into the soil, where microbes face C and nutrient-poor conditions, may stimulate soil respiration (SR) through positive priming effects37 and nutrient absorption55. However, declines in these key C and N cycling processes are likely due to the depletion of C sources caused by excluding plant roots56. All the individual functions and EMF show similar responses to plant C inputs, suggesting that PR rather than litter manipulations can significantly affect EMF in this semi-arid grassland.

The relationship between EMF and microbial diversity has been extensively studied, but the contributions of different functional groups have not reached a consensus57,58. Soil microbiomes with greater taxonomical diversity usually contain greater redundancy effects and more stable microbial associations that underpin complex ecosystem functions. However, fungal diversity regulates soil EMF more than bacterial diversity in temperate cropland, a boreal forest, and semi-arid grassland31,59. In contrast, bacterial community play dominant roles in alpine meadows and forest ecosystems4,29. This discrepancy may be explained by the fact that bacteria are primarily responsible for decomposing easily degradable substrates, whereas fungi contribute significantly to the breakdown of recalcitrant C compounds29,42. However, a stronger positive correlation of EMF with fungal diversity rather than bacteria and protists under changing C inputs observed in our study indicates that fungi may be more pivotal in driving ecological functions in soils.

Different microbial groups simultaneously play crucial roles in soil EMF, a division of metabolic labor among soil microbes can result in complementarity among various ecological processes60,61. For instance, rhizosphere-associated microbes may complement each other by spanning soil pores to facilitate the movement of other microbes to acquire new resources and maintain ecosystem functions62,63. In our study, both PR and litter manipulations can alter microbial associations because of the depletion of root- and litter-derived C sources. In accordance with our Hypothesis 3, protistan community play a little role in regulating EMF under changing C input pathways, likely due to the small changes in their diversity, composition, and microbial associations with bacteria and fungi. On the one hand, the decomposition of SOC and the mineralization of nutrients require the cooperation of diverse microbes including bacteria, fungi, and protists with specific enzymes23,53. Although soil protists have been reported to stimulate litter breakdown and soil CO2 release in different temperature regimes21, these processes (both C cycling and N provisioning) may be predominantly mediated by bacterial and fungal metabolisms, and thus have no effect on protistan community through the bottom-up regulatory pathway. On the other hand, it has been suggested that the tighter fungal-bacterial interactions and more complex ecological associations can enhance ecosystem functioning in grasslands61, therefore, the role of protists in the microbial communities might be limited due to the relatively lower proportions in the co-occurrence networks. However, we must realize that the unchanged diversity of soil protists and the weaker correlation with EMF do not necessarily mean the decoupling of protists and soil functions. In contrast, more efforts need to be made to estimate the potential linkage between them using more accurate manipulations such as the isotope tracing method.

Over all, these findings demonstrate that changes in soil microbial diversity under altered C inputs from plant roots and litter have positive associations with EMF, but protists may not be as sensitive to changing C inputs because of their diverse metabolic genotypes. Therefore, fungi and bacteria may play more significant roles than protists in driving ecological processes. Under the scenarios of global change, this work highlights the importance of plant C inputs on the assembly and functions of soil microbiomes, which greatly mediate the complex ecological processes.

Methods

Site description

This study site is located in a semi-arid grassland on the Mongolian Plateau, in Duolun County (42°02′N, 116°17′E, 1324 m a.s.l.), China. Long-term (1953-2019) mean annual precipitation is 382 mm with 90% falling from May to October. The mean annual temperature is 2.4 °C ranging from −17.5 °C to 18.9 °C. The sandy soil at this site is classified as Chestnut (Chinese classification), or Haplic Calcisols (Food and Agricultural Organization of the United Nations classification). Plant species are dominated by grass species (Agropyron cristatum, Stipa krylovii, Leymus chinensis, and Cleistogenes squarrosa) and forb species (Artemisia frigida, Heteropappus altaicus, Potentilla tanacetifolia, and Potentilla acaulis).

Experimental design

A randomized complete block design experiment was established in 2006 with six treatments including a control (C), LA, LR, PR, LA plus PR (LAPR), and LR plus PR (LRPR) (Supplementary Fig. 7). Each treatment has six replicates. Thirty-six 3 m × 4 m plots were arranged in 6 × 6 matrix, with a 1-m buffer zone between any two adjacent plots. For the LR plots, dead plant materials were gently removed using a rake in late October each year, avoiding disturbing surface-soil structure as little as possible. For the LA plots, the added plant litter was collected from the adjacent LR plots, and the average added litter mass was 0.93 kg plot−1 (Supplementary Fig. 8; 77.5 g m−2) averaged over the 14 years and all the LR plots. For the PR plots, native vegetation was removed prior to the establishment of this experiment, and germinating plants were continuously uprooted monthly by hand. Averaged over the experimental years, the mean plant aboveground biomass was 2.89 kg plot−1 (240.6 g m−2). After finishing litter manipulation and PR treatments, all the plots were covered with a 5-cm nylon net to hold plant litter on the ground.

Soil sampling and analysis

Surface-soil (0–10 cm) surveys were carried out in May (early growing season) and August (middle growing season) after 13 years of experiment in 2020. Each sample was a composite of two cores (7.5 cm diameter) randomly collected from each plot. After sieving (<2 mm) and removing organic debris, soil samples were transported to the laboratory immediately and divided into two parts: one was stored at 4 °C for analyzing soil chemical properties, and the other was stored at −80 °C for DNA extraction.

Soil-filled tins were placed in a drying oven at 105 °C for over 12 h, and then soil moisture was calculated. Soil temperature at the depth of 10 cm was detected in situ using a thermocouple probe (Li-8100-201, Li-Cor Inc., Lincoln, NE, USA) three times a month. Soil NO3−-N and NH4+-N were analyzed with a continuous flow injection analyzer (TRAACS2000, Hamburg, Germany) after extraction with 2 mol L−1 KCl. Concentrations of soil DOC and N (DON) were extracted with Milli-Q water, and detected using an automated total organic C analyzer (Elementar vario TOC, Elementar Co., Germany). SOC was also measured by TOC analyzer. Soil total C (TC) and N (TN) were analyzed on a Vario EL III Element Analyser (Elementar, Hanau, Germany), and the ratio of TC to TN (C/N ratio) was calculated.

Microbial activities relating to soil N and C cycling processes of PNR, Nmin, NA, and SBR were analyzed by the chlorate inhibition method, acetylene reduction method, and microcosm incubation methods, respectively64,65,66. In situ SR and heterotrophic respiration (HR) were detected using a CO2 flux chamber attached to the LI-8100 Soil CO2 Flux System (Li-Cor Inc., Lincoln, NE, USA) three times a month during the growth period from May to October each year.

DNA extraction, qPCR, and high-throughput sequencing

Soil DNA was extracted from 0.5 g soil using MoBio UltraClean™ soil DNA isolation kit (MoBio Laboratories, Inc., San Diego, CA, USA) according to the specification. The quality of DNA extraction was examined with a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). Extracted DNA samples were subjected to high-throughput sequencing to profile microbial communities targeting bacterial 16S rRNA genes V4 region (F-515/R-907)67, fungal ITS2 region (F-ITS7/R-ITS4)68, and protistan 18S rRNA genes V9 region (F-1389/R-1510), respectively69. Polymerase chain reactions (PCR) were performed at four parallels in 25 μl mixtures (12.5 μl of Premix TaqTM, 0.5 μl of each primer (10 μM), 2 μl of template diluted 5 folds of extracted DNA, and 9.5 μl of sterilized ddH2O). Negative control samples were also included in the PCR assay to ensure that the reaction system was not contaminated. Purified amplicons were sequenced by Illumina MiSeq PE 300 × 2 sequencer (Magigene Biotechnology Co., Ltd., Guangzhou, China).

Bioinformatics

Amplicon sequencing data were processed using the standard operation procedure of Quantitative Insights into Microbial Ecology (QIIME) 2. Paired ends of raw sequences were merged and quality-filtered using USEARCH (v10.0). The 100% similarity of correct biological reads (i.e., zero-radius operational taxonomic units, ZOTUs) was picked using the unoise3 command with default parameters68. The sequencing raw data of bacteria, fungi, and protists were assigned according to the SILVA (v132), UNITE (v8.0), and Protist Ribosomal Reference (PR2 V12.0) databases, respectively. After strict quality filtering and removing ZOTUs belonging to plant, animal, and unidentified eukaryotes, the total count of 3,785,527, 5,782,458, and 1,711,125 sequences were clustered into 25,251 bacterial ZOTUs, 3330 fungal ZOTUs and 2869 protistan ZOTUs. However, the sequencing depth of protistan dataset is relatively low, probably due to the barren soils in the semi-arid grassland or the sequencing method. In order to reduce the influence of such issue, the rarefaction curves of soil bacterial, fungal, and protistan communities under different treatments were presented in the supporting information (Supplementary Fig. 9). The raw sequencing reads of bacteria, fungi, and protists were deposited in the NCBI Sequence Read Archive (SRA) under the accession numbers PRJNA938597, PRJNA938601, and PRJNA938602.

Statistics analyses

Two-way repeated-measures analysis of variance (ANOVA) was employed to explore the effects of PR, litter manipulations (L) (LA and LR), years, and their interactions on soil microclimate, chemical properties, microbial activity, diversity, and compositions by SAS 9.0 (SAS Institute). The main effects of LA and LR on the environmental variables were calculated as \(\frac{{LA}-C+{LAPR}-{PR}}{C+{PR}}\times 100 \%\), and \(\frac{{LR}-C+{LRPR}-{PR}}{C+{PR}}\times 100 \%\), respectively, whereas the main effects of PR on those were determined as : \(\frac{{PR}-C+{LAPR}-{LA}+{LRPR}-{LR}}{C+{LA}+{LR}}\times 100 \%\).

The alpha (α)-diversity as revealed by Chao1 richness of microbial communities was calculated at ZOTU level after rarefying to a minimum number of reads per sample. The beta (β)-diversity of bacterial, fungal, and protistan communities was assessed by nonmetric multidimensional scaling (NMDS) based on Bray-Cutis distance matrices and permutational multivariate analyses of variance (PERMANOVA) analysis with the Vegan package in R (v 3.4.2). Differential abundant genera of bacteria, fungi, and protists under the main effects of LA, LR, and PR were determined by the volcano plots using the DESeq2 package in R70. The co-occurrence networks were constructed with the Molecular Ecological Networks Analysis Pipeline platform and presented using Gephi (0.9.2) software71. The topology parameters of the soil microbiome networks were calculated in the Gephi platform. The robustness was used to reflect microbial stability and was evaluated with ggClusterNet and igraph packages in R. In order to determine the general effects of treatments on soil microbial community, all networks were constructed by combining microbiomes originating from the two sampling stages.

Both biotic and abiotic factors were analyzed by the RandomForest modeling in R (v 3.4.2) with the Vegan and RandomForest packages to assess the significant predictors of changes in microbial activity and diversity. Fourteen soil parameters supporting EMF in this study were separated into three components including (i) Microbial activities relating to soil N and C cycling processes: PNR, Nmin, NA, and SBR; (ii) Soil C storage: DOC, SOC, TC, HR, and SR; and (iii) Nutrient provisioning: NO3−-N, NH4+-N, DON, TN, and C/N ratio. All the functional parameters were standardized using the Z-score transformation to obtain the averaging multifunctionality in SPSS (IBM SPSS Statistics 22)72. Mantal test of soil bacterial, fungal, and protistan communities with soil properties and microbial activity, as well as the heatmap analysis of the Spearman correlations between environmental functions and microbial diversity and community structure, were conducted with vegan and ggplot2 packages in R (v 3.4.2). SEM was constructed to analyze the effects of PR, LA, and LR on bacterial, fungal, and protistan diversity as well as EMF in AMOS 17.0. To ensure the fit of SEM Chi-square (P > 0.050), Akaike information criterion, goodness-of-fit index (>0.90), and root mean square error of approximation (<0.050) were calculated.