The experiment was carried on at the Institute of Tibetan Plateau at Southwest Minzu University (32°48'N, 102°33'E, 3,500 m above mean sea level), located about 4 km north of Hongyuan County, Sichuan, China. The mean annual temperature and precipitation is 1.1°C and 700 mm, respectively (Liu et al., 2020). The coldest month is January (~9.7°C) and the warmest month is July (11.1°C). The majority of rainfall is concentrated in warm season (June to September). The soil type is alpine meadow soil (Chinese Soil Taxonomy Research Group, 1995). The type of vegetation is classified as alpine meadow and comprises sedges (most frequently Kobresia cristata, Scirpus pumilus), grasses (most frequently Elymus nutans, Poa pratensis, and Agrostis species), and forbs (dominated by Saussurea species and Anemone species). The study site has been used for trial plots under different land use patterns caused by different grazing practices since 2007.

There are various vegetation structure and composition in the experimental areas that have different grazing disturbance histories. The common approach referred to as space-for-time substitution (Pickett, 1989) was used to infer successional patterns of alpine meadow. The advantage of this approach is that plant community variation can be investigated well without long-term observations of a single plot, which need take decades to finish (Li et al., 2015). The belt transect method was used to measure vegetation in the experimental area in August, 2018. We randomly set 10 belt transects with a length of 200 m each. Three quadrats (0.5 × 0.5 m) were randomly selected in every 10 m of each sample strip. In total, there were ~600 quadrats (samples).

To quantify the quantity and value of ecosystem functions, we chose 12 plant and soil variables including AGB, BGB, SR, plant density (PD), plant height (PH), ST, SM, SBD, STC, STN, STP, and STS in the current study. Based on the selected indicators, the plant growth index (PGI), soil carbon accumulation index (SCI), soil nutrient cycling index (SNI), and ecosystem multifunctional index (EMF) was calculated via an average method as formula 7:

where N is the number of measured functions and equal to 12 in this study, x_i represents the measured value of function i, f(xi) is the standardization of x_i via Z-score standardization.

Clustering detection (Ward method) was performed to reveal the trend and direction of plant community succession of alpine meadow (Supplementary Figure S1). The data for the vegetation and soil properties of the 30 quadrats associated with each stage of succession were used for data statistical analyses.

All data analyses were performed in R version 4.1.2 (R Core Team, 2021). Data distributions and normality was examined through Shapiro-Wilk goodness-of-fit test. First, general linear model from the “ASREML” package (Gilmour et al., 2015) to test the effects of succession stage on the SR, EV, AGB, BGB, PD, PH, SBD, SM, ST, STC, STN, STP, PGI, SCI, SNI, EC_SP, EMF, and the biomass proportion of toxic plants. Thereafter, Tukey's HSD test in the “AGRICOLAE” package (de Mendiburu, 2014) was conducted to assess the differences among means. Second, correlation analysis was performed using the “GGCOR” package (Huang, 2019) to acquire the relationships among SR, EV, PGI, SCI, SNI, and EMF, along with their relationships to EC_SP (Mantel test). Third, simple regression analysis was used to compute the relationships between EC_SP and EMF at different successional stages, along with their relationships to the biomass proportion of toxic plants. Fourth, gradient regression analysis in “GBM” package was used to analysis the relative importance of biotic and abiotic factors to EC_SP and EMF. In the nonlinear model, Friedman's H-statistic was calculated to evaluate the relative strength of the interaction effect (Li et al., 2021). Finally, heatmap analysis from the “PHEATMAP” to realize correlations between SR, EC, ecosystem function indexes (PGI, SCI, SNI, and EMF) and plant, soil properties. All figures were drawn using Origin 2021b and the “GGPLOT2” package.

Biomass ratios of plant functional groups was used to quantify the stages of succession (Figure 1). The succession series is as follows: Stage 1 (sedges, 46%; desirable forbs, 28%; grasses, 15%), Stage 2 (desirable forbs, 36%; grasses, 25%; sedges, 24%), Stage 3 (desirable forbs, 38%; undesirable toxic plants, 29%; sedges, 23%), and Stage 4 (undesirable toxic plants, 44%; desirable forbs, 37%; grasses, 11%). From the tall Kobresia community with erect growth to the short undesirable toxic forbs community (e.g., Gentianaceae, Ranunculaceae) with prostrate growth; the vertical, perennial desirable forbs community represented the transition between these two stages (Figure 1).

The SR, EV, AGB, PD, PH, and STN of Stage 4 was significantly (P < 0.05) lower than that of the other three stages (Figures 2, 3). In contrast, SBD and ST in Stage 4 was significantly (P < 0.05) higher than that of the other three stages (Figure 3). The gamma diversity of each stage (from Stage 1, 2, 3, and 4) was 61, 55, 47, and 49 species, respectively (Figure 2C). There was no significant difference (P > 0.05) in SM among four stages (Figure 3).

The SCI and SNI of Stage 3 were both significantly (P < 0.05) higher than that for the other three stages (Figures 4B,C). The EMF and extent of EC_SP were both at their lowest levels in Stage 4 (Figures 4D,E). Compared with the EMF in Stage 4, this index increased by 102.6, 89.8, and 207.6% in Stage 1, 2, and 3, respectively (Figure 4D). Compared with the extent of EC_SP in Stage 4, there was a 20.0, 16.8, and 21.2% increase in coupling in Stage 1, 2, and 3, respectively (Figure 4E). Biomass proportion of toxic plants in Stage 3 and Stage 4 was significantly (P < 0.05) higher than that in Stage 1 and Stage 2 (Figure 4F).

In general, there was a significant positive relationship (R² = 0.57, P < 0.05) between the extent of EC_SP and EMF (Figure 5B). Relationships between EC_SP and EMF were at variance in different stages (Figure 5A). Coefficients of correlation between EC_SP and EMF in Stage 1 and Stage 2 were higher than those in Stage 3 and Stage 4 (Figure 5A). EC_SP and EMF were both significantly (P < 0.05) and positively related to SCI and SNI in Stage 1 and 2 and to the soil nutrient index in Stage 4 (Figure 5A). In Stage 3 and Stage 4, EC_SP and EMF both had significant relationships with PGI and SR (Figure 5A).

Quantification of sequence of community succession is central for raveling out the function and structure of alpine meadow (Bagchi et al., 2013; Avolio et al., 2015). The best way to characterize and identify the process of succession is by analyzing the changes in vegetation composition over the course of a long-term, single-site experiment (Kahmen and Poschlod, 2004; Abbas et al., 2019). However, long-term records of trajectories of community change do not exist, either empirically or experimentally, probably because, at least for plants, the process of succession from pioneer to climax species may take 25–150 years to complete (Porensky et al., 2017). Generally, using a space series to replace a time course is an accepted method for analyzing community succession (Duan et al., 2008; Hu et al., 2016). Plant functional group (PFG) pathway could eliminate the disadvantages of individual plant species approach because of each PFG has its specific and vital role in the grassland ecosystem by altering the stability and functions (Bermejo et al., 2012; Sun and Wang, 2016). Accordingly, we recommend using the biomass proportion of plant functional groups as an effective and integrated index for evaluating the succession trajectory of alpine meadows. Our results demonstrated that the community successional series of alpine meadow were divided into four stages with a trend over time from tall erect-growth sedges community (e.g., Kobresia) to short prostrate-growth undesirable toxic forbs community (e.g., Gentianaceae, Ranunculacese) (Figure 2). These findings were also consistent with the actual situation determined by experienced observation of the field work (Zhang Z. et al., 2021) and can be viewed as a general degenerate community succession series for alpine meadows for restoration practices. The succession series can contribute to restoration practices in three aspects as following: (i) further clarification of community structure development, (ii) more cognition of species interactions and of nutrient dynamics, and (iii) better understanding of the species transformation between successional stages and of how these stages fit together into successional trajectories (Walker and del Moral, 2009; Zhang et al., 2015).

During the process of plant succession induced by grazing, the structure and function of an alpine meadow are affected by changes in the species composition and richness, soil nutrient availability, and a number of abiotic factors including environmental filtering and dispersal mechanisms (Niu et al., 2016). In this study, the changes in across the successional stages can be categorized into the three following types: (i) The decay type is characterized by a decrease in the dominance of a given species or dominance in conjunction with an increase in grazing intensity, as exemplified by Kobresia graminifolia, Gramineae species (Poa pratensis, Stipa aliena, and Festuca sinensis), and Fritillaria cirrhosa; (ii) The growth type is characterized by an increase in the dominance of a given species in conjunction with increased grazing intensity. This was noted for Anemone rivularis, Gentiana macrophylla, and Anaphalis lacteal and (iii) The fluctuating type, which included Carex thibetica, Scirpus pumilus, Elymus nutans, and Potentilla fragarioides, undergoes changes that are not consistent across the stages (Figure 2).

In a system such as this alpine meadow, the grazing intensity, and selective foraging of the yaks should be the major factors driving change in the species composition (Dong et al., 2015). Rational grazing can promote the regeneration and development of erect bunch grasses (e.g., Scirpus pumilus, Elymus nutans, and Poa pratensis) by the fragmentation of tussocks, deposition of manure, removal of standing dead plant tissues, and the stimulation of soil water and nutrient transformation and utilization (Sun et al., 2017; Wang et al., 2018). Also, toxic plants and therophyte are easily germinate and become dominant in the open space among the tussocks resulted from livestock trampling or at excrement patches (Zhang et al., 2020), which is consistent with their dominance under heavy grazing practices (i.e., Stage 4 in this study). During early succession stages, there is no selective advantage for livestock to specialize on specific species, for most of them are edible and desirable. Hence, more plenty of rapid growth plant species with low defense chemicals would be dominant owing to the low consumption of resources to herbivory (Wang et al., 2010; Milchunas and Vandever, 2013). In contrast, resource availability is relatively at low level in the late succession. Given these conditions, slow-growing, well-defended toxic weeds should dominate due to the great cost of resources to herbivory (Bosc and Pauw, 2020).

Plant regressive succession weakens the nutrient cycle and energy flow of ecosystem, ultimately leading to discordance between two dominant components (i.e., soil and plant), as well as in maladjustment of ecosystem multifunctionality in alpine meadow (Wang et al., 2020a). Special standards for grassland restoration emphasize the need for comprehensive evaluations (Gann et al., 2019), yet, plant–soil feedback studies in grassland degradation succession remains controversial and how successions stages shape ecosystem multifunction is still unclear. In consequence, using comprehensive indicators such as the extent of ecosystem coupling, ecosystem multifunctionality, and their relationships should provide a broader understanding and sight of the succession process and mechanism induced by grazing.

The difference in EC_SP, EMF and their relationships may reflect the differences in comprehensive changes of plant succession (Figure 5; Supplementary Figure S2). In this study, the values of EC_SP and EMF in Stage 4 were the lowest compared to other three stages (Figure 4). This most likely occurred because EC_SP and EMF were driven by different indicators at different stages (Cortois et al., 2016). In the Stage 1 and 2, soil factors (i.e., STC, STN, and SM) played dominant roles in influencing EC_SP and EMF. Meanwhile, both plant (PD and PH) and soil factors (STC, STN, SM, and SBD) determined the performance of EC_SP and EMF. in the Stage 4 (Supplementary Figure S2). The site of Stage 3 was associated with the highest values for the SCI and SNI, which may help to explain the highest level of EC_SP and EMF (Figure 5A). Taken together, high ecosystem coupling may show a higher sensitivity response to STC, STN, SM, and STP (Figure 6A). High ecosystem multifunctionality may show a higher sensitivity response to STC, SM, PD, and AGB (Figure 6B). Previous studies had used observational data to indicated that ecosystem functions differ in their sensitivity to grazing pressure (Ren et al., 2018) and restoration practices (Wang et al., 2020b). In this case, the lower ecosystem coupling and multifunctionality in Stage 4 induced by grazing may be caused by several mechanisms. First, heavy foraging activity of yak encouraged unpalatable broad-leaved plants (e.g., Euphorbia esula and Gentiana macrophylla) to invade. Second, overgrazing by yak leads to alpine meadow degeneration which result in rodent infestation and further grassland degeneration. Furthermore, the deeper and nutrient-poor soil was transferred to the ground level and mixed with topsoil through burrowing activity of rodent, which could contribute to the decrease of SOC and STN in the topsoil (Yu et al., 2017). Also, the latest research found that microbes also mediate the shift in ecosystem multifunctionality (Wang et al., 2021). To further illustrate the mechanism, how soil microbe affects ecosystem coupling and ecosystem multifunctionality in each succession stage needs to be further explored. In addition, our results also suggested that coefficients of correlation between ecosystem coupling and the ecosystem multifunctionality in the earlier stages were higher than those in the later stages (Figure 5B).

Toxic plants can exhibit self-protective mechanisms in degraded grasslands and promote positive feedback between plants and soil (Zhang et al., 2020). For instance, toxic plants with well-developed root systems can stabilize sand and extract nutrients from the coarser soil (Estapé et al., 2013; Li et al., 2014). Our study also emphasized the critical role of toxic plants during the late stage of succession (Stage 4). We found that both ecosystem coupling and ecosystem multifunctionality showed a significant increase with an increase in the biomass proportion of toxic plants in Stage 4 (Figure 7). In the case of this alpine meadow, we suggest that the expansion of toxic plants is a consequence of their greater suitability form a successional perspective. We thus conclude that the heavily degraded region of this alpine meadow needs toxic plant species, and these species are vital for ecological restoration, as well as for maintaining ecosystem function and coupling.