Plant Functional and Phylogenetic Diversity Regulate Ecosystem Multifunctionality in Semi-Arid Grassland During Succession

The relationship between biodiversity and ecosystem multifunctionality (EMF) is crucial for understanding the processes of ecological restoration in semi-arid regions. However, partitioning the relative influence of various biodiversity attributes, namely taxonomic, functional, and phylogenetic diversity, on EMF during secondary succession is still unclear. This study aimed to bridge the gap by employing field measurements and the chronosequence approach at 21 plots with different stand ages and precipitation conditions on the Loess Plateau of China. For diversity indices, we calculated the Shannon–Wiener diversity index, Simpson’s dominance index, Pielou evenness index, community weighted mean (CWM), functional variance (FDvar), and Faith’s phylogenetic diversity (PD) based on the empirically measured composition and traits of plant species. The EMF was expressed as the averaged value of eight function variables (including aboveground biomass, root biomass, soil total carbon, total nitrogen, and total phosphorus content, soil organic carbon, available nitrogen and available phosphorus content). The results showed that species evenness and CWM of leaf dry matter content (LDMC) significantly increased yet the CWM of specific leaf area (SLA) decreased with stand age, indicating the resource-use strategy of the plants became more conservative through succession into its later stages. The EMF increased with both stand age and mean annual precipitation. The structural equation model revealed that stand age, soil water content (SWC), and the multiple diversity indices altogether accounted for 56.0% of the variation in the EMF. PD and the CWMs of plant height and LDMC had positive effects on the EMF, and the FDvar of leaf nitrogen had negative effects on EMF. However, the Shannon Wiener diversity had no significant effect on the EMF. Our results suggest that functional and phylogenetic diversity are more important than taxonomic diversity in predicting EMF, and that multidimensional biodiversity indices should be jointly considered to better predict EMF during the succession of semiarid grasslands.


INTRODUCTION
It is increasingly recognized that ecological restoration is vital for mitigating the loss of biodiversity and improving ecosystem functioning in semiarid regions (He et al., 2016;Zuo et al., 2017;Jing et al., 2019). The shifts in community composition and biodiversity during secondary succession were usually associated with changes in ecosystem processes and functions (Purschke et al., 2013;Kelemen et al., 2017). Previous studies of the biodiversity-ecosystem functioning relationship have mainly focused on some aspects of diversity and individual ecosystem function, such as net ecosystem productivity, soil carbon stock, and soil erosion (Everwand et al., 2014;Zhu et al., 2015;Huang et al., 2020). The single factor studies usually neglect the tradeoffs between multiple ecosystem functions (Allan et al., 2015) and limit or even bias our understanding of their relationship (Byrnes et al., 2014). The concept of ecosystem multifunctionality (EMF) has provided a way to tackle the complex issue. Multifunctionality is defined as the ability of ecosystems to simultaneously provide multiple ecosystem functions or services (Nelson et al., 2009;Xu Y. et al., 2021b). EMF is essentially an integrative metric that embodies multiple ecosystem functions, providing a more comprehensive understanding of the overall functioning of an ecosystem (Manning et al., 2018).
Biodiversity is a complex multifaceted concept that entails different scales and entities (Pavoine and Bonsall, 2011). During the past decades, the mostly mentioned diversity indices contained taxonomic, functional and phylogenetic diversity (Purschke et al., 2013;Yuan et al., 2016). Species richness was frequently used as the taxonomic diversity to predict ecosystem productivity in grasslands (Aarssen et al., 2003;Gross et al., 2014). Yet, many studies have suggested to consider functional and phylogenetic diversity in predicting ecosystem functions, which complements the relational pattern through functional traits and evolutionary history (Flynn et al., 2011;Steudel et al., 2016). Thus, multiple aspects of diversity should be jointly considered in predicting the biodiversity-ecosystem function relationships (Yuan et al., 2016).
Numerous studies have discussed the relationships between biodiversity and ecosystem multifunctionality (Hector and Bagchi, 2007;Meyer et al., 2018). According to the analysis of published data from grassland biodiversity experiments in Europe, greater numbers of species were required to maintain multifunctionality (Hector and Bagchi, 2007). Plant taxonomic diversity was an important predictor of the multifunctionality related to soil processes in the semiarid grassland restoration (Guo et al., 2021).However, taxonomic diversity metrics emphasized the number of species per se and neglected the underlying information of species identity (Díaz et al., 2004). Recent work suggested that functional diversity indices, which were calculated based on plant traits related to plant ecological strategies, were better at predicting ecosystem multifunctionality, especially those related to nutrient cycling, biological productivity and soil fertility, than species richness in drylands (Gross et al., 2017). The analysis of data from 123 drylands found that both phylogenetic and functional diversity can serve as robust predictors of ecosystem multifunctionality (Le Bagousse-Pinguet et al., 2019). Given that in most cases it is impractical to measure all plant traits related to every ecosystem function, some researchers have advocated using phylogenetic diversity to substitute for functional diversity in predicting ecosystem functions (Srivastava et al., 2012). However, an experiment along an aridity gradient across the grasslands of Inner Mongolia found that phylogenetic diversity was not as effective as functional diversity in predicting EMF related to the cycling and storage of soil carbon, nitrogen, and phosphorus (Yan et al., 2020). Moreover, phylogenetic diversity might be applicable only when the relevant plant traits show a phylogenetic signal (Srivastava et al., 2012;Narwani et al., 2015); hence, it is not a sound proxy for functional diversity when the traits are not phylogenetically conserved (Steudel et al., 2016). Thus, it remains ambiguous which facets of biodiversity are most apt to predict EMF.
As succession progresses, the different facets of diversity undergo differential patterns of change (Lohbeck et al., 2012;Vellend et al., 2013;Yuan et al., 2016). Both increasing and decreasing trends of taxonomic diversity during grassland succession have been reported (Purschke et al., 2013;Lazzaro et al., 2020). The phylogenetic alpha diversity significantly increased from the early-mid to late-successional stages but with no significant changes detected in the early and early-mid successional stages in an arable-to-grassland chronosequence more than 270-years-long (Purschke et al., 2013). How functional diversity indices shift during succession has also been investigated and observed (Lohbeck et al., 2012;Zhang et al., 2015). For example, the community weighted mean (CWM) of specific leaf area (SLA) decreased, whereas the CWM of leaf dry matter content (LDMC) increased during succession in both grassland and forest ecosystems, indicating a shift of resourceuse strategies from acquisition to conservation during ecological succession (Garnier et al., 2004;Buzzard et al., 2016;Kelemen et al., 2017). Meanwhile, the individual ecosystem functions and multifunctionality were found augmented during grassland restoration in the Loess Plateau (Guo et al., 2021). The relationship of taxonomic diversity and multifunctionality during succession in temperate forest were recently studied (Lucas-Borja and Delgado-Baquerizo, 2019). However, the effect of different facets of biodiversity on EMF along the course of succession in semi-arid grasslands has yet to be resolved.
The Loess Plateau in China was one of the most degraded regions worldwide (Zhao et al., 2013). It is characterized by semiarid region, where soil water content (SWC) is the main limiting factor of ecosystem processes. SWC was mainly controlled by precipitation in semi-arid regions. The constraint of precipitation influenced both the speed of succession and the accumulation of soil organic carbon (SOC) during succession in semi-arid area (Gabarrón-Galeote et al., 2015). Ecological restoration projects such as the "Grain for Green Program" (started in 1999) have enabled the natural restoration of abandoned cropland in this region (He et al., 2016;Liu et al., 2017). This long-term restoration background has generated a series of abandoned fields with different stand ages, providing us with an excellent opportunity to explore the successional dynamics of biodiversity and ecosystem multifunctionality during the grassland restoration process. In this study, we investigated the three facets of biodiversity and their correlation with EMF in naturally restored grasslands along a successional gradient. We proposed a conceptual model to rigorously test the relationships between different aspects of biodiversity and EMF ( Figure 1). Specifically, the objectives of our study were to (1) determine the dynamics of the three key facets of biodiversity (i.e., taxonomic, functional, and phylogenetic diversity) and EMF during succession; and (2) explore which aspects of diversity can better explain the changed EMF in semi-arid grassland restoration.

Study Sites
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).

Experimental Design and Sampling
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  Lat-Long, latitude and longitude; MAP, mean annual precipitation; Precip, precipitation during the growing season; SWC, soil water content; BD, bulk density. Soil property parameters are for soil depths 0-30 cm. a The precipitation during the growing season from May to August at four sites in 2017 (YG, and AS), 2018 (SM) and 2020 (SD).
Frontiers in Environmental Science | www.frontiersin.org January 2022 | Volume 9 | Article 791801 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.

Ecosystem Function and Plant Trait Measurements
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 3 ) 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.

Taxonomic Diversity
We calculated the Shannon-Wiener diversity index (H), Simpson's dominance index (D), and Pielou evenness index (E) to quantify the taxonomic diversity. The equations of these indices are as follows: where P i is the proportional abundance of species i, and S denotes the number of species in each quadrat. The taxonomic diversity indices were calculated using the "diversity" function of the "vegan" package in R software (Oksanen et al., 2015).

Functional Diversity
The most widely used functional diversity indices based on single-trait are CWM and functional variance (FDvar). The CWM indices calculate the mean values of plant traits in the community, using species relative abundance as weighting factor (Díaz et al., 2007). It represents the trait values of the dominant species, and can be used to test the mass ratio hypothesis (Grime, 1998). The mass ratio hypothesis proposed that the contribution of a species to the ecosystem function is proportional to its biomass, and the functioning of ecosystem mainly depends on the values of the functional traits of the dominant species (Grime, 1998). In contrast, FDvar allows for a test of the niche complementarity hypothesis that the diverse distribution of functional traits can improve resource use in a complementary way (Mason et al., 2003;Conti and Díaz, 2013). We calculated CWM and FDvar, in each quadrat, based on plant trait values and the relative abundance of each dominant species.
where S is the total number of species in the quadrat, w i is the relative abundance of the ith species and the biomass for each species as the relative abundance, and x i is the trait value of the ith species (Garnier et al., 2004).
In these equations, five is a scaling factor used to define the index over a range of 0-1, V is the weighted variance of character x, w i is the relative abundance of the ith species and x i is the trait value of the ith species (Mason et al., 2003).

Phylogenetic Diversity
First, we standardized the species' names and taxonomic nomenclature according to the Plant List (http://www. theplantlist.org/) using the "plantlist" package in R v3.4.3. We constructed a phylogenetic tree using the phylogeny proposed by Zanne et al. (2014), which was built based on seven gene regions and had branch lengths in Phylomatic Version3 (http:// phylodiversity.net/phylomatic/) (Webb and Donoghue, 2005). Overall, the constructed phylogenetic tree contained 62 species (Supplementary Figure S2). We calculated Faith's phylogenetic diversity (PD) using the "picante" package in R. The PD index, defined as the cumulative branch length of a phylogenetic tree connecting all species in the community (Faith, 1992), represents the phylogenetic difference among species and has been applied in many ecological studies (Steudel et al., 2016;Xie et al., 2018). In addition, we also evaluated the phylogenetic signals of the seven functional traits, to evaluate the linkage between functional and phylogenetic diversity. The phylogenetic signal was calculated as Frontiers in Environmental Science | www.frontiersin.org January 2022 | Volume 9 | Article 791801 the "K statistic" and "Pagel's λ", both using the "picante" package (Yan et al., 2020).

Ecosystem Multifunctionality
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 .
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.

Statistical Analyses
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 (χ 2 ) test, root mean square error of approximation (RMSEA), and comparative fit index (CFI). A good model fit was indicated by χ 2 values associated with a p-value > 0.05 (suggesting that observed and expected covariance matrices are statistically different), χ 2 / 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.

Changes in the Diversity Indices During Succession
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).

Changes in EMF During Succession
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).

Structural Equation Model
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 (χ 2 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).

Shifts in Biodiversity During Succession
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 . We found that the species evenness strengthened during Frontiers in Environmental Science | www.frontiersin.org January 2022 | Volume 9 | Article 791801 6 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   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).

Improvement in Ecosystem Multifunctionality During Succession
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 Frontiers in Environmental Science | www.frontiersin.org January 2022 | Volume 9 | Article 791801 might have been constrained by less precipitation at SM (a drier site).

Impact of Functional and Phylogenetic Diversity on Ecosystem Multifunctionality
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 fastgrowing 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).

Effects of Soil Water Content
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 . 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.

CONCLUSION
Our study uncovered dynamic changes in three aspects of biodiversity (i.e., taxonomic, functional, and phylogenetic diversity) and EMF and their relationships during succession in the restoration of semi-arid grassland on the Loess Plateau. The taxonomic diversity and dominance index showed no significant changes with stand age. But the increasing CWM_LDMC and decreasing CWM_SLA over time reflected a shift in the resourceuse strategy from acquisition to conservation along the successional gradient. Notably, EMF was enhanced with stand age during grassland restoration. Stand age, SWC, and diversity indices together explained 56.0% of the variation in EMF. The FDvar_LN and CWMs of Ht and LDMC were the key functional indicators to predict changes in EMF across successional stages, providing strong evidence in support of the mass ratio hypothesis to explain the relationships between functional diversity and ecosystem multifunctionality during secondary succession. Functional and phylogenetic diversity can predict multifunctionality much better than taxonomic diversity. The quantification of multifunctionality in abandoned fields provides a more comprehensive assessment of ecological restoration in semi-arid grassland. We suggest that multidimensional biodiversity indices should be jointly considered to better gauge and understand the relationships between biodiversity and ecosystem functions in semi-arid regions.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.