The study was conducted in the central area of the Loess Plateau, Shaanxi province, China. This area is characterized by a typical hilly and gully landform with a semi-arid continental climate. The main soil type is a loessial soil, which is classified as a Calcaric Cambisol according to the FAO system (Zhong et al., 2021). We selected four sites of restored grasslands in the study area (Supplementary Figure S1), named Yan’gou (YG), Ansai (AS), Suide (SD), and Shenmu (SM). The mean annual precipitation (MAP) is around 523 mm at YG and 405 mm at SM, and the mean annual temperature ranges from 8.4°C to 9.8°C (Table 1).

We used the chronosequence approach (Dana and Mota, 2006) to investigate changes in community diversity and EMF of semi-arid grasslands during their secondary succession. Through interviews with local landowners, we selected parcels of land differing in elapsed time since their abandonment but with a similar land-use history. No other activities have been conducted and no disturbances have occurred since abandonment. In this way, we selected 21 plots of different stand ages from four sites, including six plots at AS and five plots each at YG, SD, and SM (Table 1). These plots had been abandoned for 3–40 years. The dominant species composition at different stand ages were consistent with the succession sequence previously reported in this region (Du et al., 2005). Hence, these plots were appropriate for analyses of grassland secondary succession. The stand conditions, including the direction, degree, and position of the slopes, were similar and consistent among the plots. Five randomly distributed quadrats (1 × 1 m) were established in each plot (i.e., 21 plots × 5 quadrats = 105 sampling quadrats in total). Sampling was carried out in August 2017 at YG, July 2017 at AS, August 2018 at SM, and August 2020 at SD. We found no significant difference in precipitation during growing seasons between different sampling years. Thus, we supposed that the different sampling years were acceptable in our study.

We recorded the name, number, and relative coverage of each plant species with one or more extant individuals present in a given quadrat. We hand-clipped all the standing vegetation in each quadrat and obtained the roots from the middle of the quadrat in three soil layers: 0–10, 10–20, and 20–30 cm, using a 10 cm-diameter soil auger. Root samples were placed in a sieve (mesh size: 0.18 mm) and washed under flowing water to remove any adhering soil. To obtain the dry weight of aboveground biomass (AGB) and root biomass (RB), the aboveground vegetation and roots were heated in an oven at 105°C for 20 min and then dried again at 70°C for 24 h.

Three soil samples were collected from random locations within each quadrat, using a cutting ring (volume of 100 cm³) inserted at the same three depths (i.e., 0–10, 10–20, and 20–30 cm). The soil bulk density and SWC were determined by oven-drying the soil at 105°C for 48 h (Fornara and Tilman, 2008). For each quadrat, the other soil samples from the same depth were thoroughly mixed to form a composite sample. The soil total carbon (TC) and nitrogen (TN) content and SOC were measured with a carbon–nitrogen analyzer (Primacs SNC100-IC-E). The available nitrogen (AN), phosphorus (AP), and total phosphorus content (TP) were determined using a continuous flow analyzer (SEAL Auto Analyzer 3).

Those species contributing more than 80% of the total aboveground biomass in each plot were designated the dominant species. For each dominant species, we measured five plant traits: plant height (Ht), SLA, LDMC, leaf nitrogen (LN), and leaf phosphorus (LP) content. Full details of functional trait measurements are provided in the Supplementary Material.

The eight soil-related function variables we measured are closely related to carbon (AGB, RB, TC, and SOC), nitrogen (TN and AN), and phosphorous (TP and AP) storage and cycling. These variables were reliable proxies for carbon sequestration and soil nutrient cycling in previous studies (Valencia et al., 2018; Wang et al., 2019). Here we used the averaging approach to estimate the EMF (Xu Y. et al., 2021b; Guo et al., 2021), as it can provide a straightforward way to assess changes in several functions at once (Byrnes et al., 2014; Wang et al., 2019).

We performed a pairwise Pearson correlation analysis of the eight function variables to avoid redundancy in the EMF calculation. Only one of the 28 correlation coefficients was higher than 0.7 (Supplementary Figure S3), which suggested all eight variables were sufficiently independent and acceptable for deriving the EMF (Valencia et al., 2018). We then standardized all individual functions by min–max normalization (Eq. (7)), such that each transformed function had values that ranged from 0 to 1. The eight standardized functions were finally averaged (Eq. 8) to obtain the EMF of each quadrat (Wang et al., 2019). F i = F r a w i − F m i n i F m a x i − F m i n i (7) EMF = ∑ i = 1 n F i n (8) where F _i is the transformed value of the ith function, F _rawi is the untransformed value of the ith function, and the F _mini and F _maxi are respectively the minimum and maximum untransformed value of the ith function.

All variables were respectively checked for a normal distribution, using the Shapiro–Wilk test. Data that were not normally distributed were log-transformed for the statistical analyses. Pearson correlation analyses were performed to determine the changes of diversity indices and EMF in relation to stand age. Two-way analysis of variance (ANOVA) was used to test the influence of site, stand age, and their interaction on the EMF and likewise for each diversity index. After a significant F-value, pairwise differences in the means of the EMF and diversity indices among the four sites were assessed using the Scheffe test in a post-hoc comparison. These analyses were carried out in SPSS Statistics 20 software.

We used SEM to quantify the plausible impacts of stand age, SWC, and diversity indices upon EMF. We chose the model with the lowest AIC, and assessed model fitting by the chi-square (χ²) test, root mean square error of approximation (RMSEA), and comparative fit index (CFI). A good model fit was indicated by χ² values associated with a p-value > 0.05 (suggesting that observed and expected covariance matrices are statistically different), χ²/degree of freedom (df) < 2, an RMSEA <0.05, and a CFI > 0.90. We also calculated the standardized total effect of stand age, SWC, and diversity indices to infer each variable’ total contribution to EMF. The SEM was implemented using AMOS 21.0 software.

The diversity indices were significantly correlated with stand age except for FDvar of LP and Ht (Table 2). Both site-dependent effects and interactions between stand age and sites were also detected (Table 2). The diversity indices showed different trends during succession among the four sites (Figure 2). In particular, species evenness significantly increased with stand age at SD and SM while the H and D decreased at AS, whereas all taxonomic diversity indices were insignificantly correlated with stand age at YG. The PD was negatively correlated with stand age at AS and SM (p < 0.001), but it had no obvious linear trends at the other two sites (Figure 2). The CWMs of LN and LP were each significantly correlated with stand age, but the dynamic trends were inconsistent among the four sites. The CWM_LN decreased significantly with stand age at YG and SM but it increased at AS; CWM_LDMC increased significantly with stand age at YG, AS, and SD (p < 0.05); CWM_SLA decreased significantly with stand age at YG and AS, but it showed a non-linear trend at SM; CWM_Ht was significantly changed with stand age only at SD. By way of comparison, the FDvar indices were less influenced by stand age. Only the FDvar of LN and LP decreased significantly with stand age at SD and LDMC decreased at SM (Figure 2).

There were significant differences in the diversity indices among the four sites (Table 2; Figure 3). The H, D, and PD were significantly lower at SD (p < 0.05), yet the evenness was similar among the four sites (Figure 3). The CWMs of SLA and LN were both highest, while the CWM_LDMC was lower at SD (p < 0.05) and CWM_Ht was lowest at AS (p < 0.05). In stark contrast, the CWM_LP did not differ significantly among the four sites (Figure 3). FDvar values of leaf traits were significantly lower at SD than the other three sites, and FDvar_Ht was greatest at AS (p < 0.05). The two-way ANOVAs revealed that stand age and site interacted significantly for H, E, and the CWMs of LN, LP, and SLA, in addition to the FDvars of LN and SLA (p < 0.05), but not so for PD and the other functional diversity indices (Table 2).

The EMF was significantly increased with stand age at YG, AS, and SD (p < 0.05) (Figure 4). The EMF values fluctuated with stand age at SM, precluding a clear trend of increase discernible during succession. The EMF increased mainly at the mid-later successional stages, indicating the effects of grassland restoration upon EMF may emerge after several years of ecosystem recovery. The mean value of EMF was significantly lower at SM than at either AS or YG (p < 0.05) (Figure 4).

We combined all the data to quantify the direct and indirect effects of stand age, SWC, and diversity indices on EMF using the SEM approach. The final model showed that these factors collectively explained 56.0% of the variation in EMF (Figure 5A), and the model (χ² = 7.56, df = 9, p = 0.579; RMSEA <0.001, CFI = 1.00) fitted well with the data. Notably, we found a strong direct effect of stand age upon EMF (p < 0.001) (Figures 5A,B), yet the direct impact of SWC on EMF was not significant (p > 0.05) (Figure 5A). Beyond direct effects, the functional and phylogenetic diversity indices had significant mediating effects for the ways by which stand age and SWC influenced EMF, whereas taxonomic diversity did not (p > 0.05). More specifically, stand age had significant indirect effect on EMF by regulating the values of CWM_LDMC. We detected an indirect effect of SWC upon EMF by influencing both PD and CWM_Ht; meanwhile, PD had a negative indirect influence on EMF by altering the FDvar_LN (p < 0.01). Overall, functional and phylogenetic diversity indices outperformed taxonomic diversity in predicting EMF in terms of the standardized total effects (Figure 5B).

Biodiversity can be used as a surrogate to evaluate the response of a community to restoration succession (Jing et al., 2019; Guo et al., 2021). Nevertheless, the trend for taxonomic diversity’s response along the successional gradient is not always consistent among studies. For instance, Shannon diversity of plants increased along the secondary succession from abandoned agricultural lands to forest in Mediterranean semi-arid regions (Heydari et al., 2020). Yet a decline in taxonomic diversity happened during grassland restoration on the Qinghai-Tibet Plateau, due to the inhibition of litter coverage time on seedling emergence and establishment (Liu et al., 2016). We found that the species evenness strengthened during succession, but the diversity and dominance indices changed negligibly with stand age (Figure 2). The Shannon diversity indices increased at first and then decreased with stand age, peaking at around the 20th year at the four sites (Supplementary Figure S4). This reduced taxonomic diversity in the late successional stage may arise from greater coverage of vegetation and accumulation of litter in the early stage, which would constrain new seedlings’ access to light and slow the recruitment of other species (Lamb, 2008; Jing et al., 2013).

Shifts in functional traits are usually indicative of responses to changes in local resource availability and the consequent trade-offs between acquiring vs. conserving resource-use strategies of plants (Kahmen and Poschlod, 2004; An and Shangguan, 2014). Those species predominant in the early successional stage are fast-growing ones characterized by an acquiring resources-use strategy (Buzzard et al., 2016). They usually have a high SLA, high nitrogen concentration in their tissues, and low LDMC to efficiently acquire more resources and grow more quickly (Lavorel and Garnier, 2002; Wright et al., 2004). However, the opposite plant strategy and trait characteristics can be found in later successional stages (Kelemen et al., 2017). We found inconsistent temporal dynamics among the different functional diversity indices during the succession of semi-arid grasslands on the Loess Plateau, with only the CWMs of LDMC and SLA showing significantly increasing and decreasing trends across the four sites, respectively (Figure 2). The shifts in CWMs of LDMC and SLA provide compelling evidence for the resource-use strategy shift from acquisition to conservation at the community level during succession, which is consistent with some previous studies (Buzzard et al., 2016; Kelemen et al., 2017).

The overall trend of CWM_SLA was decreasing as succession progressed, but its dynamics varied depending on the site: it decreased with stand age at YG (p < 0.05) and SD (p < 0.001), but did not respond significantly to succession at AS or SM (Figure 2). Both YG and SD had higher values of CWM_SLA in the early successional stage, indicating the efficient photosynthesis by fast-growing plants at the beginning of succession. During the successional processes, fast-growing species were replaced by the slow-growing ones, resulting in the lower values of CWM_SLA at the later successional stages (Kelemen et al., 2017). The insignificant response of SLA to stand age at AS and SM was consistent with findings reported by Kahmen and Poschlod (2004). Further, the dynamics of CWMs exhibited an opposite trend at different sites; e.g., decreasing CWM_LN at YG and SM but increasing at AS (Figure 2). The differential response of functional traits, especially SLA and LN, might due to the initial site conditions in soil water and nutrient levels (Douma et al., 2012).

Our study confirmed the enhancement of EMF during grassland restoration, which is consistent with the previous findings in grassland and forest succession (Lucas-Borja and Delgado-Baquerizo, 2019; Guo et al., 2021), suggesting the improvement of carbon sequestration and soil nutrient retention by grassland restoration. We also observed that most individual functions increased with stand age during succession, but those variables related to soil phosphorus seemed unresponsive to succession (Supplementary Table S1). Similar results were found in grasslands of the Loess Plateau, in which their soil total phosphorus and available phosphorus content lacked significant correlations with stand age (Wang et al., 2011; Guo et al., 2021). Compared with the research which mainly focused on the dynamics of individual ecosystem functions during vegetation succession (Garnier et al., 2004; Ohtsuka et al., 2010), the present evaluation based on EMF in our study offers a more reliable way to assess the restoration of semi-arid grasslands.

We found that EMF was significantly enhanced at YG, AS, and SD (p < 0.001), but it had no obvious trend at SM site (Figures 2, 4). Similarly, research on forest succession has found that multifunctionality can increase in some but not all sites (Cruz-Alonso et al., 2019). Results of our correlation analysis showed that, with stand age, the AN and AP are significantly decreased at SM site while AGB and RB are significantly increased, though other functions had no significant dynamics (Supplementary Table S1). The conflicting trends in individual functions at SM likely drove the insignificant response of multifunctionality to succession there. The significantly increased AN at AS and SD is consistent with previous findings in the Loess Plateau (Wang et al., 2011; Guo et al., 2021). The trend of AN decline at the SM site may be due to the slow decomposition rate in the late successional stage (Weintraub and Schimel, 2005). The SM site featured a higher CWM_LDMC and lower CWM_SLA (Figure 3), which may provide low quality litter and accordingly a slow decomposition rate (Rosenfield and Muller, 2020). In addition, the litter decomposition rate can change along a long-term precipitation gradient (Schwartz et al., 2007). In another study, soil AN content was significantly positively related to precipitation along the Northeast China Transect (Wang et al., 2005). Therefore, the accumulation of AN during succession might have been constrained by less precipitation at SM (a drier site).

Previous studies have found that species richness, trait distribution, and phylogenetic diversity can improve multiple ecosystem functions and services, which are closely related to nutrient cycling, biological productivity, and carbon sequestration (Gross et al., 2017; Xie et al., 2018; Guo et al., 2021). In our study, the EMF increased with stand age and it was associated with the functional diversity indices and PD, whereas taxonomic diversity had no significant impact on EMF. Our findings demonstrate that functional diversity and PD can play more important roles in predicting EMF during succession in semi-arid grasslands.

The SEM results showed that stand age exerted a significant, direct effect on EMF and an indirect effect by altering the CWM_LDMC (Figure 5A). LDMC, which is strongly related to nutrient availability and plants’ relative growth rate, has been used before as a functional effect trait to predict ecosystem functions (Garnier et al., 2004; Kazakou et al., 2006). For example, the communities dominated by fast-growing species (i.e., these having a lower LDMC and higher SLA) in early successional stages tend to have higher annual aboveground net primary productivity (Garnier et al., 2004; Kohler et al., 2020). However, the AGB significantly increased with increasing CWM_LDMC in the understory and at the whole-community level in a subtropical forest (Ali and Yan, 2017), which implied that conservative species with high LDMC drive high AGB in resource-limited environments (e.g., where light resources are limiting). In addition, LDMC can influence the quality of litter produced, and further affect litter decomposition and soil nutrient retention in grasslands (Pakeman et al., 2011). Increasing trends of soil organic matter content and soil nitrogen content have been observed with increasing CWM_LDMC in grasslands (Garnier et al., 2004; Grigulis et al., 2013). Thus, we suppose CWM_LDMC might have promoted the EMF related to carbon sequestration and nutrient cycling in our study. In addition, CWM_Ht was found positively controlled by SWC, and it further improved the EMF (Figure 5A). Plant height is generally related to light resource acquisition, and so a higher CWM_Ht can increase the accumulation of aboveground biomass, and accordingly soil carbon and nitrogen storage (De Deyn et al., 2008; Zuo et al., 2016; Wang et al., 2021). The positive effect of CWM indices on multifunctionality supports the mass ratio hypothesis, which posits that ecosystem functions are mainly determined by the most abundant species in a given community (Grime, 1998).

Further, only FDvar_LN remained in the final SEM, having a negative effect on EMF (Figure 5A). Negative relations between FDvar and ecosystem functions have been observed in previous studies (Conti and Díaz, 2013; Zuo et al., 2016). In our study, the Spearman correlation analysis showed that FDvar_LN was negatively correlated with CWM_LN (r = –0.241, p = 0.013); hence, the leaf nitrogen distribution might be concentrated towards high values with low divergence in the community. The effect of the FDvar_LN on EMF ultimately depended on the associated CWM, and this is consistent with the mass ratio hypothesis (Conti and Díaz, 2013).

Phylogenetic diversity has received increasing interest by those studying the relationships between biodiversity and ecosystem functions, but there is currently no consensus of the effect of PD upon ecosystem functioning (Srivastava et al., 2012). The functional traits did not show any detectable phylogenetic signals (Supplementary Table S2), which means there were hidden and unmeasured phylogenetically conserved traits that are related to carbon sequestration and nutrient cycling in our study (Cadotte et al., 2009; Le Bagousse-Pinguet et al., 2019). Our study did find that PD has a direct positive effect on EMF and an indirect effect via increasing the FDvar_LN (Figure 5A), which suggests that PD might capture some important unmeasured functional traits, thereby influencing the ecosystem functions (Zirbel et al., 2019). Recent research showed that phylogenetic diversity is a better predictor of aboveground biomass and decomposition rate than taxonomic diversity in restored grasslands (Zirbel et al., 2019). The positive effect of PD on ecosystem functions may due to the niche complementary effect, whereby an increasing evolutionary distance among species promotes increasing niche space in the community (Srivastava et al., 2012). This evolutionary divergence inevitably generates ecological differentiation among species, which can further reduce the overlap of resource use among co-occurring (sympatric) species and improve their resource-use efficiency as well as ecosystem functions (Cadotte et al., 2008; Cadotte et al., 2009; Xiao et al., 2020). For example, PD indirectly promoted multifunctionality by augmenting the structural complexity of stands and light capture efficiency in a temperate forest (Yuan et al., 2020). The positive relationship between PD and multifunctionality suggests that a higher diversity of phylogenetic lineages can enhance multiple ecosystem functions (Le Bagousse-Pinguet et al., 2019; Li et al., 2020).

Previous studies have also revealed that different ecosystem processes are affected by different species, and maintaining ecosystem multifunctionality requires greater biodiversity than studies focusing on individual function (Hector and Bagchi, 2007; Zavaleta et al., 2010). The bivariate correlation showed that Shannon index was significantly positively correlated to EMF (r = 0.227, p = 0.020). The insignificant pathway between Shannon diversity and EMF in our final model indicated that plants’ taxonomic diversity had no significant effect on ecosystem functions when concurrently controlling for functional and phylogenetic diversity indices. This implied that an increasing taxonomic diversity may affect the underlying functional diversity and thereby promote ecosystem functions (Mouillot et al., 2011). Thus, functional and phylogenetic diversity indices played more important roles than taxonomic diversity in predicting EMF in our study. Future restoration projects in semi-arid regions should consider the phylogenetic distance of species or the functional composition in the community, and not just manipulate the species richness in the grassland (Zhu et al., 2015).

Water availability is the primary limiting factor for determining community structure, ecosystem processes, and the course of succession in arid and semi-arid regions (Gabarrón-Galeote et al., 2015; Zuo et al., 2018). In our study, the significant differences in multiple diversity indices and EMF uncovered among the four sites (Figures 3, 4) were partly related to soil water availability. We found that the mean values of Shannon diversity were significantly higher at YG and AS (p < 0.05) but Pielou evenness was similar among the four sites (Figure 3). The pattern was consistent with the findings for other grasslands of the Loess Plateau whose Shannon–Wiener diversity index was positively correlated with MAP (p < 0.05) but not the evenness index (Zhang et al., 2018). In our study, both CWM_LN and SLA were significantly higher at SD, which is the drier site in terms of MAP and precipitation during the growing season in the sampling year (Precip) (Table 1; Figure 3). Our results are thus consistent with previous findings that the LN is negatively correlative with SWC in a semi-arid region (Zuo et al., 2018). We also found that CWM_Ht was lowest at AS which had the lowest SWC among the four sites (Figure 3), suggesting that plant height responds rapidly to soil water availability in semi-arid regions.

Ecosystem functions related to carbon sequestration and nutrient cycling are strongly affected by SWC and precipitation in arid and semi-arid regions. For example, NPP was positively correlated with MAP in arid and semi-arid grasslands in China and globally (Chen et al., 2016; Liu et al., 2019). More precipitation promotes the growth of aboveground vegetation and further affects the organic matter inputs and accumulation of soil carbon (Gabarrón-Galeote et al., 2015). The SOC stock markedly increased along the precipitation gradient on the Loess Plateau, suggesting the stronger ability of carbon sequestration under wetter environmental conditions (Han et al., 2018). More plant litter input to the soil carbon and nitrogen pools might increase the soil nitrogen transformation rate along the precipitation gradient (Feyissa et al., 2021). Further, we found that the mean value of EMF was significantly higher at YG than either SD or SM (Figure 4), a result consistent with another reported finding that MAP positively affected carbon and nitrogen cycling and ecosystem multifunctionality in a dryland (Le Bagousse-Pinguet et al., 2019). The positive effect of MAP on ecosystem multifunctionality, which is closely related to dynamics of plant growth, the soil carbon stock, and nutrient cycling (Xiong et al., 2016; Xu Y, et al., 2021a), would suggest that carbon sequestration and nutrient cycling perform better at wetter sites. In this way, precipitation can simultaneously enhance multiple ecosystem functions in arid and semi-arid regions.