The study site is located in the town of Dawu, Maqin County of Qinghai Province, China, the core zone of the Three Rivers Headwater Area (Figure 1). The geographical position of the study site is 32^°31′–35°37′N in latitude, 96°54′–101°51′E in longitude, and the average elevation is over 4,000 m. The annual average temperature is 0.54°C (from 1997 to 2017 statistics) and the annual average precipitation is 518.2 mm (The coefficient of variation of the inter-divisional dynamics is 24%). The average temperature in the coldest month (January) is −11.1°C, and the average temperature in the hottest month (July) is 10.7°C. The soil type is loam (Gao et al., 2019). The primary vegetation type in the study area is an alpine meadow with Kobresia pgymaea as the dominant species, and Kobresia humilis, Stipa aliena, Festuca ovina, Poa pratensis, Elymus nutans, and Leontopodium nanum, as companion species. Harsh environments make alpine meadow habitats sensitive and fragile to climate change and human disturbances. Around one-quarter of the alpine meadow in the study sites have been degraded as “Black Beach” with severe soil erosion due to lowered coverage of pioneer plants such as Aconitum flavum, Ligularia virgaurea, Pedicularis kansuensis, and Potentilla anserina (Shang et al., 2018).

The severely degraded alpine meadow, “Black Beach” has been largely rebuilt with revegetated grasslands. Cultivated grassland (also called artificial grassland) is a new grassland community formed by planting one or more plants on the completely degraded grassland by using comprehensive agricultural technology. Over past decades, large areas of the “Black Beach” degraded grasslands have been restored with cultivated grassland under the governmental programs “Returning Grazing Land into Grassland” and “Grazing Ban” (Dong et al., 2013). With the support of a series of these restoration projects, numerous demonstration sites of cultivated grasslands with different restoration times have been established in the study area (Gao et al., 2019).

According to relevant information and as reported in the literature, in this study, we did the samplings based on the chronosequence/space-for-time approach in the cultivated grasslands, which were established in 2000 (A1), 2002 (A3), 2004 (A4), 2006 (A6), 2008 (A8), 2010 (A10), 2012 (A12), 2014 (A14), 2015 (A16), and 2017 (A18) and severely degraded alpine meadow (Black Beach, A0), and non-degraded healthy alpine meadow (ND). All the sampling plots in this study were located in the same demonstration study sample area accompanied by similar topography and geographical location. The planting area of each cultivated grassland was 50 hm². The cultivated grasslands by sowing E. nutans were all degraded grasslands on “Black Beach” and were built through deep turning, harrowing, fertilizing [45 kg/hm² of diammonium phosphate (NH₃HPO₄)], sowing (sowing amount was 45 kg/hm² local seed), covering soil, crushing, fencing, and other agronomic measures.

In mid-August 2017, five quadrats with an area of 1 m × 1 m were randomly selected from each plot for plant sampling and species identification. According to the classification method on the functional groups of alpine grassland plants (Wang et al., 2004), and in combination with the morphological and physiological characteristics of alpine grassland plants, we divided the species into four functional groups: Grass, Cyperaceae, Leguminosae, and the remaining species were classified as Forb (Yu et al., 2017; Zhao et al., 2017). The number of species, density, height, coverage, and frequency of individual species were estimated. In each plot, one subplot in the size of 0.5 m × 0.5 m was selected for sampling aboveground plant biomass, which was brought back to the laboratory and dried in an oven at 65°C to constant weight for measurement. After plant sampling, soil core samples were taken to a depth of 20 cm in each subplot.

Based on commonly used soil properties analysis test theory in the laboratory, we conducted experiments on the measurement of the following soil properties for each soil sample: total organic carbon (TC) and nitrogen content (TN) were tested by the elemental analyzer (Euro EA 3000, Italy), available phosphorus (AP) and available potassium (AK) were tested by an Inductive Coupled Plasma Emission Spectrometer (Germany). We also quantified soil moisture.

We used the variables related to diversity, stability, and the plant interspecific association (network analysis) mentioned above for data analysis. In detail, for characterization of species diversity and dominance structure, Pielou’s evenness index (E), Shannon–Wiener diversity index (H), and Simpson’s dominance (λ) index were calculated. Associated coefficient (AC) was used to characterize the degree of association among species (Fang et al., 2012). The range of AC is from -1 to 1. The closer the AC value is to 1, the stronger the positive association between species is. The closer the AC value is to -1, the stronger the negative association between species is. When the AC value is 0, the species are completely independent. So, we classified the inter-species relationships into three trends, positive inter-species association, negative inter-species association, and random inter-species association, respectively. The calculation of the above indices was conducted in R using the “vegan” package.

Stability (PS) was calculated as the inverse coefficient of variation, which was calculated via the ratio of the average biomass to the standard deviation (Strecker et al., 2016; Zhou et al., 2019). The importance value (IV) of individual species was calculated via the formula of IV = (C′+H′+P′)/3, where IV is the important value of each species, C′ is relative coverage, H′ is relative height, and P′ is relative frequency.

The network analysis method was used to identify the species connection in the plant community. In the network, each species was considered as a node. We used the SPAA package in R to calculate the inter-species relationship matrix, which was treated as the basic data of sub-group analysis. The modular analysis was performed to build the block model and simplify through the CONCER algorithm (convergent correlation) (Xie et al., 2019). Before drawing two complete sub-networks, we obtained the simplest network diagram that maintained the integrity of the network using UCINET 6.5 software, which is a social network analysis program containing a graphical illustration tool called Net-Draw, a program used widely to obtain the simplest network diagram of the plant community (Johnson, 1987; Thomas et al., 2007). The density is an indicator of the closeness of membership and betweenness centrality represented the ability of a member to act as an intermediary, that is, to what extent a node is located in the middle of other nodes in the network (Qin and Li, 2009). In the network of vegetation communities, the density reflects the closeness of the inter-species relationship in a vegetation community. The networks of interacting species are usually described in terms of hardiness and stability (Kleyer et al., 2018). Network density (D) was calculated as the ratio between the total actual lines presented in the network graphics and the potential lines (Wolfe, 1995). Betweenness Centrality (C_B) in the network was calculated via the formula (Wolfe, 1995):

Where g represents the number of nodes in the network diagram, g_ik represents the shortcut of member i to member k, C_B(n_i) is the intermediary of node n_i, and C_B(n*) is the maximum node intermediation. Based on this, we calculated the network density (SND) and betweenness centrality (SNB) of the plant interspecific linkage network.

The data were processed and analyzed using Microsoft Excel 2016, SPSS 22.0, R version 3.5.1, and R Studio version 1.1.456. Network analysis was performed using UCINET 6.5 software. Network analysis graphs were drawn with Net Draw.

We used 28 indicators including density (D), network indexes containing the complexity of plant interspecific association (network density: SND) and betweenness of plant interspecific linkage network (SNB), diversity indexes, which include the evenness indexes of the total community (E-T), Grass (E-G), Cyperaceae (E-C), Legume (E-L), and Forb (E-F); the diversity of whole community (H-T), Grass (H-G), Cyperaceae (H-C), Legume (H-L), and Forb (H-F); the proportional dominance of Grass (G_λ), Cyperaceae (C_λ), Legume (L_λ), and Fords(F_λ); the proportion of inter-specific associations of plant community containing positive inter-species associations (Positive); negative inter-species associations (Negative); random inter-species associations (Random); and stability indices which include the stability of the whole community (T-PS), Grass (G-PS), Cyperaceae (C-PS), Legume (L-PS), and Forb (F-PS) and soil properties containing the carbon (C) and nitrogen (N), available phosphorus (AP), available potassium (AK), and soil moisture of the soil samples.

To test the above-mentioned indicators as the influence of independent variables on the overall stability and to avoid redundancy, we used correlation analysis to filter the variables. For the same type of indicators with a correlation coefficient greater than 0.70, one of them was retained (Zhou et al., 2019). Finally, a structural equation model (SEM) was adopted to test the hypothetical model, which was established and tested using AMOS 4.0.

In each plot, the species with a 100% probability of occurrence are the grass E. nutans, the legume Astragalus membranaceus var. membranaceus, and the Forb Potentilla multifida. The dominant species was E. nutans (an indigenous grass planted) along restoration years except for the 8th year. The sub-dominant species changed greatly with the restoration time and the range of average species number of each plot (Table 1). Compared with “Black Beach,” the species richness of cultivated grassland decreased by 9.2% in the first recovery year and increased from the second restoration year until the 8th restoration year. Over 10 years of restoration, plant species richness in the cultivated grassland tended to continuously increase. Through all restoration years, the number of species in the plant community of artificial grassland increased. Over 16 years, there was a high similarity of species composition between the cultivated grasslands at different restoration years and the non-degraded native grassland, and low similarity of species composition between the cultivated grasslands at different restoration years and the “Black Beach.”

In the first 3 years after the cultivated grassland was built to replace the extremely degraded grassland, the proportion of comprehensive dominance of other functional groups except Cyperaceae increased. In the next 7 years, the proportion of dominance in Cyperaceae increased, while the proportion of other functional groups showed a fluctuating decrease trend. With the increase of recovery time, the dominance of Forb was close to 10%, while the dominance proportion of the other three functional groups was close to 30%, respectively (Table 1).

We classified the inter-species relationships into three categories, positive inter-species association, negative inter-species association, and random inter-species association (Figure 2A). The proportion of species with random associations was higher than those of species with positive and negative associations in the plant community of all types of grasslands except “Black Beach.” The proportion of species with positive associations was slightly greater than those with negative associations in all types of grasslands. Throughout the restoration time, the proportions of species with both positive and negative associations generally tended to decrease and those of species with random associations tended to increase, although there were great fluctuations of inter-species relations after 4 and 10 restoration years.

The SNB of extremely degraded grasslands and cultivated grasslands planted for 1 year was greater than 4. The SNB values of other cultivated grasslands were all about 2–3, while the value of non-degraded grasslands was the minimum of 1.268. The complexity of plant interspecific association (species network density) was positive for non-degraded alpine meadow and cultivated grassland communities were positive in the 6th, 12th, and 16th years, while in the “Black Beach” degraded grassland it was negative. Plant interspecific association indicated that a few species highly influenced the entire network (Table 2).

In the first 3 years of grassland restoration, the plant stability in the revegetated areas was significantly lower than that of the “Black Beach” degraded grassland (Figure 2B). The stability of Grass was increased after 6 years of restoration. The change in the productivity of the Forb was quite low. The similarity coefficients between the stability of Grass, Cyperaceae, Leguminosae, Forb, and the whole community were -0.158, 0.531, 0.659, and 0.368 respectively. The trend of stability of the whole community was similar to that of legumes, which was reflected by the greatest correlation among all plant groups.

The structural equation model was used to test the relationship and mediation between variables to study the influence mechanisms of plant community composition and diversity on grassland stability during the restoration process of alpine grassland. Since there was strong collinearity in the evenness and diversity of plant communities in different functional groups, eight factors including the evenness index of Grass (E-G), Cyperaceae (E-C), Leguminosae (E-L), Forb (E-F), Grass (H-G), Cyperaceae (H-C), Legumes (H-L), and Forb (H-F) were extracted and expressed as four principal components. We defined them as D-L, D-C, D-G, D-F, respectively, representing the diversity information of these four plant functional groups (Legume, Cyperaceae, Grass, Forb) (Table 3). Similarly, using the same analytical method, the total soil carbon (TC) and total soil nitrogen (TN), available phosphorus (AP), available potassium (AK), and soil moisture were extracted as two principal components. We defined them as soil nutrient and moisture, respectively (Table 4).

The results of the structural equation model indicated that compared with other functional groups, grass diversity had a more significant effect on the stability of the whole community during the restoration process of artificial grasslands (Figure 3). The restoration process of artificial grassland can directly affect or indirectly affect the species composition characteristics of the community (i.e., the network structure based on the interspecific relationship) by changing the soil nutrient conditions, and thus significantly affect the overall stability of the community.