Grazing Exclusion Changed the Complexity and Keystone Species of Alpine Meadows on the Qinghai-Tibetan Plateau

Grazing exclusion is an effective approach to restore degraded grasslands. However, the effects of grazing exclusion on keystone species and the complexity of plant community were poorly investigated. Here, we conducted a field survey among different grazing exclusion durations, i.e., Grazing, grazing exclusion below 5 years, grazing exclusion with 5 years, grazing exclusion with 7 years, and grazing exclusion over 7 years, in alpine meadows on the central Qinghai-Tibetan Plateau (QTP). The complexity and keystone species of alpine meadows were analyzed by a network analysis. The results showed the following: (1) The species richness did not change, but aboveground biomass and the coverage of the plant community tended to increase with the extension of the grazing exclusion duration. (2) The soil nutrients, i.e., total nitrogen, total organic carbon, available nitrogen, and available potassium, remained stable, while the soil bulk density decreased under grazing exclusion conditions. (3) There was a hump-shaped change of the complexity (i.e., average connectivity and average clustering coefficient) of the plant community along with the extension of the grazing exclusion duration. Moreover, the keystone species were different among the grazing exclusion treatments. Based on the complexity of the plant community and the changes of keystone species, the optimum duration of grazing exclusion for alpine meadows should be between 5 and 7 years. Our results suggest that besides the productivity, the change of the complexity and keystone species of plant community should be considered when grazing exclusion is adopted to restore the degraded alpine meadows.


INTRODUCTION
Grazing exclusion is an effective approach to restore degraded grassland ecosystems around the world, which could promote rapid recovery of the productivity of a grassland ecosystem (Al-Rowaily et al., 2015;Fedrigo et al., 2018;Golodets et al., 2010). The Qinghai-Tibetan Plateau (QTP) is called the roof of the world and plays crucial roles in the supply of water resources and animal by-products for the population of Asia. Land degradation is an urgent ecological issue on the QTP (Teng et al., 2018) and seriously affects the life quality of people in this region Li X.-W. et al., 2015). To restore the degraded alpine meadows, a grazing exclusion strategy was adopted since 2001 on the QTP (Dong et al., 2007). The effects of grazing exclusion on plant diversity and productivity were widely discussed and concluded that grazing exclusion is a benefit to the restoration of degraded meadow ecosystems on the QTP Wang et al., 2018). Studies reported that the positive effects of grazing exclusion on diversity and productivity of alpine meadow could be detected at a specific duration of grazing exclusion Cheng et al., 2016). From the perspective of plant recovery, some studies identified that the optimum duration of grazing exclusion for alpine meadows on the QTP is 6 years  or 4 years (Sun et al., 2020), and others suggested that grazing exclusion should be ceased after 6-10 years (Xiong et al., 2016). A 5-year experiment revealed that grazing exclusion could affect plant communities (i.e., phenology and reproduction) by changing soil physical properties (e.g., soil moisture) . However, how grazing exclusion causes change of plant community structure is not clearly understood yet. In particular, changes in the complexity and keystone species of the plant community, which typically reflect the change of community structure, are still unknown under grazing exclusion.
Since the first concept of keystone species was proposed (Paine, 1969), ecologists developed the concept and found that the keystone species may not be the dominant species but play crucial roles in the structure and functioning of biotic systems at either ecosystem-level or community level (Mills et al., 1993;Davic, 2003;Lu et al., 2017;Hale and Koprowski, 2018). There are two traditional ways to identify keystone species, i.e., experimental removal and comparative method (Power et al., 1996). The removal method is devastating and requires a precondition: there are limitless species in an ecosystem, which is impossible in a real one. The reliability of the results of the comparative method could be reduced by the heterogeneous conditions at different locations (Power et al., 1996).
The connection among components of the ecosystem, which could be used to describe the complexity of the ecosystem (Pimm, 1984), is a primary character of an ecosystem. This character provides a new perspective to identify the keystone species of an ecosystem. The network analysis, which is a powerful approach to evaluate the connection of species in a complex ecosystem (Strogatz, 2001), has been used to assess community structure and identify the key species according to the highly connected nodes in a given ecological network (Banerjee et al., 2018;Meyer et al., 2020). In comparison, though the complexity and keystone species are inherent characteristics of the plant community and could evaluate the efficiency of grazing exclusion in degraded grasslands, there are few studies that discussed these issues based on the network analysis.
In this study, we used the network analysis to identify the complexity (measured by the average connectivity and average clustering coefficient of the network) and the keystone species of plant communities along a grazing exclusion gradient in the central QTP. We hypothesized that the duration of grazing exclusion could strongly affect the complexity and keystone species of alpine meadows. Moreover, changes in the complexity and keystone species could indicate the optimum duration of grazing exclusion.

Field Investigation and Laboratory Measurement
This study was conducted in the central QTP, i.e., Nagqu County of the Tibet Autonomous Region (Figure 1), which distributes a mass of alpine meadows that are dominated by Kobresia pygmaea. During 1956-2013, the mean annual temperature and precipitation of this region are around -2 • C and 600 mm, respectively. Many alpine meadows are fenced to excluding livestock grazing for land restoration in this region. In 2014 summer (i.e., the growing season of plants), we investigated four groups of grazing exclusion pastures, i.e., grazing exclusion over 7 years (GEO7), grazing exclusion 7 years (GE7), grazing exclusion 5 years (GE5), and grazing exclusion below 5 years (GEB5). Meanwhile, the continuously grazed plots were investigated outside each fenced pasture (denoted by Grazing) (Figure 1). Three individual pastures were investigated in the group of GEO7, GE7, GE5, and GEB5, respectively. In each sampling pasture, around five quadrats (1 m × 1 m) were randomly arranged over 20 m away from the fence to record the plant species as well as the canopy cover (the percentage of plant cover in the 1 m 2 quadrat) and height (measured five times to get the average value) of each species (Figure 1). In total, 122 quadrats were investigated. Moreover, the herbage was clipped at ground level from the quadrat. The clipped herbage was oven-dried for 48 h at 70 • C and measured with electronic scales (0.01 g) to evaluate the aboveground biomass of the plant community in each sampling quadrat. In total, 122 quadrats were investigated ( Table 1).
At each clipped quadrat, a ∼40 cm 3 soil sample was collected using a cutting ring (with a diameter of 2.5 cm and a depth of 5 cm) to measure the bulk density (BD) according to the method used by Li et al. (2018). A total of 122 BD samples were collected. Meanwhile, 0-to 15-cm soil cores with a diameter of 5 cm were collected to test the contents of soil nutrients. A total of 122 soil nutrients samples were collected. The total nitrogen (TN), available nitrogen (AN), total organic carbon (TOC), available phosphorus (AP), and available potassium (AK) were determined according to the methods used by Wang et al. (2014).

Network Construction and Statistical Analysis
There were three steps to construct the plant community network in this study, i.e., calculation of the importance value (IV) of species, construction of the nodes and edges, and network analysis and visualization.
Firstly, the IV of each plant species was calculated according to its relative cover (RC) and relative height (RH) in a community . The IV of species i was calculated by the formula below: Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 1 | Location of the study area. GEO7, GE7, GE5, and GEB5 represent grazing exclusion over 7 years, 7 years, 5 years, and below 5 years, respectively.
Where RC i = C i / C i , RH i = H i / H i . C i and H i are the cover and height of species i, respectively. C i and H i are the total cover and total height of a community (i.e., a 1 m × 1 m quadrat), respectively. The summed IV equals 1 in a community. Therefore, high IV means high dominance of a plant species in the community. The IV of plant species was summed according to the category of different functional groups (i.e., Gramineae, Cyperaceae, Leguminosae, and Forbs).
Secondly, the nodes and edges that were used in the network analysis were constructed for Grazing, GEB5, GE5, GE7, and GEO7, respectively. The plant species were eliminated if their frequency was less than 5% in communities, and the remaining species were designed as the nodes of the network. Pearson's correlation coefficient between each two plant species (frequency over 5%) was calculated based on the IV of species. To simplify TABLE 1 | Details on the samples conducted in continuously grazed pastures (Grazing), grazing exclusion below 5 years (GEB5), grazing exclusion 5 years (GE5), grazing exclusion 7 years (GE7), and grazing exclusion over 7 years (GEO7).

Treatments
Number the network, the significantly correlated relations (we set the significance level as α = 0.2, because of the limitation of sampling size) was taken to be designed as the edges of the network. The absolute value of Pearson's correlation coefficient was used as the weight of edges. Based on the nodes and edges mentioned above, an undirected random network was constructed and visualized in Gephi 0.9.2. The node degree is the number of edges incident to the node ( Table 2), which is a primary index to imply the importance of the node in a network (Deng et al., 2012), while the hub index, which is calculated based on the in-degree metric (a kind of node degree) of the node, could measure the importance and authority of a node in a network (Kleinberg, 1999). Therefore, the node (i.e., species) that had the highest hub value was identified as the keystone species of a plant community in this study. The average connectivity (or average degree) of a network indicates the strength of species interactions in a community ( Table 2), and a higher average connectivity indicates a more complex community (Deng et al., 2012). The average clustering coefficient indicates the organization level of a network ( Table 2), and a higher average clustering coefficient means a higher organization level (Watts and Strogatz, 1998;Latapy, 2008). Thus, the average connectivity index and average clustering coefficient indices were used to analyze the complexity of plant communities along the grazing exclusion gradient in this study.
The differences of species richness, coverage, aboveground biomass, IV of functional groups, and IV of keystone species among treatments were detected by one-way ANOVA in

Indexes
Formula or explanation Average connectivity* ac = n i=1 k i n , where k i is degree of node i. n is the number of nodes.

Hub value §
The hub value of node i is calculated by an iterative algorithm based on the in-degree of i.
, where l i is the number of links between neighbors of node i. k i is the number of neighbors of node i.
where cc i is the clustering coefficient of node i. *Formula is cited from Deng et al. (2012). § More details please refer to Kleinberg (1999).
PASW Statistics 18.0 (IBM Inc., United States). We used regression analysis to represent the relationships between species richness and aboveground biomass and the complexity of plant community, and the relationship between the IV of keystone species and soil properties by PASW Statistics 18.0 (IBM Inc., United States). To explore the changes of plant community composition under grazing exclusion, a non-metric multidimensional scaling (NMDS) was conducted based on the IV of species to order sampling sites (Legendre and Gallagher, 2001). The NMDS was conducted by the vegan package (Oksanen et al., 2019) in R (v.4.0.3; R Core Team, 2020).

Overall Changes of Plant Communities and Soil Properties Under Grazing Exclusion
Seven families, i.e., Cyperaceae, Gramineae, Rosaceae, Compositae, Ranunculaceae, Gentianaceae, and Leguminosae, were repeatedly observed at different grazing exclusion pastures. Generally, these families could be classified into four functional groups: Gramineae, Cyperaceae, Leguminosae, and Forbs. The IV of Cyperaceae in GEO7 was significantly higher than that in GEB5, GE5, and GE7 (F = 9.97, p < 0.001), but did not significantly different from grazing treatment ( Table 3). The IV of Gramineae did not significantly change along the grazing exclusion duration (F = 0.30, p = 0.88). The highest IV of Leguminosae was detected in GE7, while the lowest one was detected in GE5 (Table 3). For Forbs, the highest IV was detected in GEB5 (Table 3). Overall, the IV of each functional group was ordered as Cyperaceae > Forbs > Gramineae ∼ Leguminosae in each treatment (Table 3).
Species richness did not change, but the coverage and aboveground biomass of the plant community trended to increase with the duration of grazing exclusion (Figure 2).
The NMDS analysis revealed that the structure of plant community in GE5 and GE7 was different from those in grazing Different lowercase letters among functional groups in each treatment (i.e., each row) means significant difference; different uppercase letters among treatments for each functional group (i.e., each column) means significant difference (α = 0.05). Grazing exclusion below 5 years, GEB5; Grazing exclusion 5 years, GE5; Grazing exclusion 7 years, GE7; Grazing exclusion over 7 years, GEO7. Frontiers in Ecology and Evolution | www.frontiersin.org conditions, while the plant community under either GEB5 or GEO7 were similar to those under grazing conditions (Figure 3).
The BD trended to decrease with the extension of the duration of grazing exclusion ( Table 4). The content of AP increased, while other soil nutrients (i.e., TN, AN, TOC, and AK) did not change among grazing exclusion durations (Table 4).

Changes in the Complexity of Plant Communities Along the Grazing Exclusion Gradient
The number of qualified connections, i.e., the significant correlation at α = 0.2, between species increased in GEB5, GE5, and GE7 when compared with the grazing condition. The number of qualified connections between species was similar to grazing condition when grazing exclusion over 7 years (GEO7) (Figure 4). Accordingly, the complexity of plant communities increased gradually when the duration of grazing exclusion was under 7 years. Then, the complexity of plant communities decreased when the duration of grazing exclusion was over 7 years (Figures 5A,D).
A quadratic model fitted the relationship between species richness and the average connectivity of the plant community. The species richness decreased, firstly, and then increased with the average connectivity of the plant community (Figures 5B,E). The aboveground biomass of plant community trended to increase with the average connectivity of plant communities ( Figure 5C).
The changing trend of the average clustering coefficient, which represented the organization level of the plant community, was similar to the changing trend of the average connectivity of the plant community among the grazing exclusion durations (Figure 5D). Species richness was higher in either low or high organization level than in a moderate organization level ( Figure 5E). A linear model fitted the relationship between aboveground biomass and the organization level of the plant community ( Figure 5F).

The Changes of Keystone Species With the Grazing Exclusion Gradient
The keystone species was Potentilla multifida under the treatments of Grazing and GEB5 (Figures 4A,B, 6). It was Kobresia pygmaea, Kobresia humilis, and Lancea tibetica in GE5, GE7, and GEO7, respectively (Figures 4C-E, 6). The IV of P. multifida and K. humilis was changed significantly, while the IV of K. pygmaea and L. tibetica did not change significantly with the increase of grazing exclusion duration ( Table 5).

The Relationship Between the Optimum Duration of Grazing Exclusion and the Succession of Plant Community
It is widely reported that grazing exclusion has little effects on species richness in grassland ecosystems (Mayer et al., 2009;Rong et al., 2014;Xiong et al., 2016). The species richness might decrease when the duration of grazing exclusion is much longer, e.g., over 8 years (Shi et al., 2013;Zou et al., 2016;Su et al., 2017). In this study, we found that the species richness did not significantly differ among each grazing exclusion treatment, which might be due to the limitation of the duration of grazing exclusion. Alternatively, this might be attributed to the loss of seed bank of some sexual reproductive species caused by irrational grazing management .
Grazing exclusion could strain competition among plant species and change the plant community composition (Amiaud et al., 2007;Harris et al., 2015;Odriozola et al., 2017;Loydi, 2019). We found that the structure of the plant community was different among GE5, GE7, and continuous grazing plots, but it was similar between GEO7 plots and continuous grazing plots. This phenomenon was consistent with experimental results obtained from a New Mexico short-grass prairie (Holechek et al., 2006). Stable species richness and similar complexity level mean a similar competition situation in plant communities, which might be the reason why the plant community was similar between long-term grazing exclusion and continuous grazing pastures in this study. Corresponding to the distribution of plant community detected by the NMDS, the complexity of plant community showed a hump-shaped change with the increase duration of grazing exclusion. The complexity caused by the increasing interactions among species could improve the stability of the plant community (Pimm, 1984). Therefore, the change of the complexity of plant community detected in this study indicated that the optimum duration of grazing exclusion was around 7 years on the central QTP.
It seems that soil nutrients had little help to explain the change of keystone species along a grazing exclusion gradient because it is generally reported that grazing exclusion does not change the contents of soil nutrients in alpine meadow ecosystems on the QTP (Lu et al., 2015;Zhang et al., 2015), which is just like what found in the present study. The change of soil physical properties under the grazing exclusion was generally observed on the QTP . In this study, we found that the BD decreased with the increase of grazing exclusion duration, which means soil moisture would vary among these pastures. Our previous study revealed that soil moisture is a limiting factor for the growth of meadow plants on the central QTP (Ganjurjav et al., 2016). Therefore, the change in keystone species might be explained by a soil moisture-related reason: the species that identified as keystone species might play crucial roles to maintain soil moisture. Under a grazing condition, the selective ingestion of livestock and wild herbivores removes a mass of plant coverage. This could increase the evaporation of soil moisture. The remaining unpalatable species become important covering to maintain soil moisture by reducing soil evaporation and surface soil temperature in alpine meadow ecosystems (Chen et al., 2016;Li et al., 2016). Consequently, the interactions between unpalatable species and other species would be inevitably enhanced. P. multifidi is unpalatable and wide-blade on alpine meadow ecosystems on the QTP. It was detected as the keystone species in Grazing and GEB5 pastures in this study. P. multifidi could play crucial roles in maintaining soil moisture under grazing and short-term grazing exclusion pastures where plant cover was low. Then, it could change the structure of plant communities on the QTP by stimulating the growth of plants whose height is low and the root is shallow Wang et al., 2019).
After a 5-year grazing exclusion, the coverage of plant community was over 70% in this study. The increase of coverage would decrease the soil evaporation. In this condition, K. pygmaea and K. humilis were detected as keystone species under GE5 and GE7, respectively. Cyperaceae species, e.g., K. pygmaea or K. humilis, are dominant species and sometimes are treated as key species for alpine meadows on the QTP (Kaiser et al., 2008;Miehe et al., 2011;Wang et al., 2017). This indicated that the alpine meadow was back to a healthy status after 5-7 years' grazing exclusion.
In agreement with other studies (e.g., Li et al., 2018), we found that long-term grazing exclusion could induce a decrease in soil BD. Consequently, surface soil would become lacunose and the strength of evaporation would increase. Species that grow close to the ground and wide-blade are more efficient to prevent surface water loss than other species under a condition of lacunose soil structure. L. tibetica, which grows very close to the ground and is wide-blade, was identified as a keystone species in the community of GEO7 in this study. Further studies in the future that focus on the role of P. multifida and L. tibetica in the maintaining of soil moisture would help to support the theory mentioned in this study. Moreover, K. humilis was a keystone species as well as a dominant species in GE7. In other treatment, the plant community was dominated by K. humilis whose IV was over 0.3, while the IV of keystone species (i.e., K. pygmaea, P. multifida, and L. tibetica) was low (major less than 0.1). These results suggested that dominance might not be a necessary condition to identify the keystone species.
Generally, the keystone species played important roles in constructing an ecological community (Paine, 1969;Meyer et al., 2020). The shift of keystone species figured that the appropriate duration of grazing exclusion was around 7 years, indicating that the succession of plant community inherently determines the optimum duration of grazing exclusion on the QTP.

A Comprehensive Strategy Should Be Adopted to Manage Alpine Meadows on the QTP
Overgrazing is blamed to be responsible for the degradation of alpine meadows on the QTP (Harris, 2010). Grazing exclusion, which is generally adopted to restore degraded grassland ecosystems around the world (Bi et al., 2018), is used to improve the coverage and productivity of alpine meadows on the QTP. However, both field investigations and remote-sensing results revealed that the relationship between the restoration level of alpine meadows and the duration of grazing exclusion is nonlinear on the QTP Wu et al., 2017). Studies showed that there is an optimum duration of grazing exclusion, FIGURE 4 | The plant community networks under Grazing (A), grazing exclusion below 5 years (GEB5, B), grazing exclusion with 5 years (GE5, C), grazing exclusion with 7 years (GE7, D), and grazing exclusion over 7 years (GEO7, E) in the central QTP. The network was constructed using high-frequency species (frequency over 5%). Plant species is the node, and the link between significant-correlated species is the edge in a network. Plant species were separated and colored according to modules, and the size of dots was determined by the hub index of species in each network. e.g., 6 years from the perspective of the plant biomass, on the QTP . In this study, we found that the optimum duration of grazing exclusion was around 7 years when it was considered from the perspective of the succession of the plant community. The similar conclusions on optimum duration, which were obtained from different perspectives, suggest that the duration should be considered when a grazing exclusion strategy was adopted on the QTP.
It is reported that temporary grazing exclusion, which was conducted in specific seasons (such as spring or fall), was beneficial to the rapid recovery of the plant productivity (Fedrigo et al., 2018). Evidence showed that reasonable livestock grazing is   beneficial to maintain the health of alpine meadow ecosystems (Medina-Roldán et al., 2012;Zhang et al., 2019). Moreover, historically, nomadism was a traditional and sustainable way to use alpine meadow systems on the QTP. Therefore, we suggest that the grazing exclusion strategy should combine with a rotational grazing system to maintain the health of meadow ecosystems on the QTP. Besides, complexity is the most typical nature of an ecosystem, which ensures the ability of an ecosystem against external disturbances (Commoner, 1971). It widely confirms that more species means higher complexity of an ecosystem (Chapin et al., 2000). In this study, there was an inverted parabolic relationship between the complexity and species richness, which means that the complexity of plant community might be determined not only by species richness but also by other factors, such as connections among species. This indicated that in addition to increasing the species number, building effective correlations among species might be another key approach to improve the quality of a degraded alpine meadow on the QTP.

CONCLUSION
Change in the complexity of plant community revealed that the optimum duration of grazing exclusion for alpine meadows was around 7 years on the central QTP. Shifting of the keystone species represented a succession process of plant community with the increase of the grazing exclusion duration. Therefore, the duration of grazing exclusion could strongly affect the complexity and keystone species of alpine meadows. Moreover, changes in the complexity and keystone species could indicate the optimum duration of grazing exclusion in the central QTP.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
HG, YZ, QG, and SD designed the experiment. YZ and HG conducted the field survey. All authors participated in the writing of the manuscript.