The grazing experiment was conducted in Nagqu city, Tibetan Autonomous Region of China, which is situated at the central QTP with latitude 31.441°N, longitude 92.017°E, and altitude 4,500 m a.s.l. The region has an arid and cold climate. The annual average temperature is 0.03°C and the annual average rainfall is 474.7 mm. In the growing season (May–September), the average rainfall is 414.7 mm and the average temperature is 7.8°C. The soil agro-type is alpine meadow soil with high proportions of clay. The main vegetation is an alpine meadow, with the plant community dominated by Kobresia pygmaea and accompanied by Kobresia humilis, Potentilla humilis, Potentilla saundersiana, Potentilla bifurca, Astragalus membranaceus, Leontopodium leontopodioides, Stipa capillata, and so on.

During the growing season of 2014–2017, the grazing experiment was carried out in a 30,000 m² paddock of the alpine meadow, which was enclosed with wire mesh and divided into 12 50 × 50 m plots (Supplementary Figure S1). Nine plots were randomly placed on the paddock for rotational grazing (RG) and were further divided into groups A, B, and C, with each group having three replicates. The other three plots which were subjected to grazing before 2014 were used as three replicates of grazing exclusion (GE). At the beginning of July 2014, 2015, 2016, and 2017, two yaks with an average weight of 105±10 kg were grazed in each replicate plot of group A for 7 days (six yaks in three replicates). Then, these yaks grazed in groups B and C in turn in the same way. Three rotation cycles lasted 63 days each year except for 2015, in which only two rotation cycles lasted 42 days due to drought, which led to the earlier withering of plants. Three plots of the same size as the RG and GE plots were randomly selected as continuous grazing (CG) plots outside the fence in 2016 and 2017. The area outside the fence contained a 37.5-ha paddock of the alpine meadow that was permanently grazed by approximately 300 yaks. In the RG and CG plots, the grazing intensity was eight young yaks per hectare.

On September 5, 2017, plant leaves of all species were collected following the method of Cornelissen et al. (2003) and Pérez-Harguindeguy et al. (2013). The sampling method has been specifically described previously (Li et al., 2019, 2020). In brief, at least 30 plants of each species were selected in each grazing regime; three to five fully expanded leaves of each individual were sampled for further leaf trait measurements. In the laboratory, the dry plant leaves of each species were ground by a ball mill (NM200; Retsch, Haan, Germany) for chemical analysis. Total carbon and nitrogen contents were measured by an elemental analyzer (Perkin-Elmer, Boston; MA, USA), and total phosphorus content was measured by an ICP-AES analyser (Thermo-Jarrell Ash Corp; MA, USA). In addition, six soil samples with a diameter of 5 cm were randomly collected at a depth of 0 to 15 cm in each plot, three samples were used to analyze the soil nutrient content and microbial community, and three samples were used to estimate the plant belowground biomass. The NH₄-N and NO₃-N of the soil samples were determined by a flow injection AutoAnalyser (AACE, Germany), and available phosphorus was determined by inductively coupled plasma spectrometry (SPECTRO ARCOS EOP, Germany). Soil total carbon content was measured using an elemental analyser (2400 II CHNS/O Elemental Analyser; Perkin-Elmer, Boston, MA, USA). Soil bulk density was determined with a cutting ring by sampling three cores in each plot at a depth of 0 to 15 cm. The detailed estimation of the aboveground and belowground biomass of the plant community can be found in our previous study (Li et al., 2020).

The soil microbial community was analyzed using high-throughput sequencing. Soil microbial DNA was extracted from each soil sample using a Fast DNA SPIN Kit for Soil (MP Biochemicals, Solon, OH, USA). The V4-V5 of bacterial 16S rRNA gene was amplified by polymerase chain reaction (95°C for 2 min, followed by 25 cycles at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 45 s) using the universal primer 515F (5′-GTGCCAGCMGCCGCGG-3′) and 907R (5′-CCGTCAATTCMTTTRAGTTT-3′). ITS rRNA of fungi was amplified by polymerase chain reaction (95°C for 2 min, followed by 32 cycles at 95°C for 30 s, annealing at 61°C for 30 s, and extension at 72°C for 45 s) with primer ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2R (5′-GCTGCGTTCTTCATCGATGC-3′). The PCR products were gel-purified with an AxyPre DNA Gel Extraction kit (Axygen Biosciences, Union City, CA, USA) and quantified by Qubit® 2.0 Fluorometer (Invitrogen Corp., Carlsbad, CA, USA), then a mixture of amplicons was used for sequencing on the Illumina MiSeq platform. The processes of controlling quality and trimming sequencing reads were performed following the method of Hong et al. (2015). Operational Taxonomic Units (OTUs) were built at 97% sequence similarity cutoff using UPARSE (Edgar et al., 2011). OTUs were assigned to taxonomic lineages using the Ribosomal Database Project (RDP) classifier within the Silva database (http://www.arb-silva.de). The total sequence number in each sample was rarefied to the minimum sequence number across all the samples. The rarefaction curves and the number of OTU reads of samples are shown in Supplementary Figure S6 and Table S3, respectively. The functional profiles of soil bacteria and fungi were assigned by FAPROTAX 1.2.3 (Louca et al., 2016) and FUNGuild (Nguyen et al., 2016), respectively.

Two-tailed paired t-tests and analysis of variance (ANOVA) were performed to test the grazing regime effect on plant species coverage and soil total carbon content. Permutational multivariate analysis of variance (PERMANOVA) based on Bray–Curtis dissimilarity indices were performed to examine the difference in plant community composition among grazing regimes. Kruskal–Wallis tests were performed to examine grazing regime effects on soil microbial diversity indices and microbial phyla relative abundance. Microbial community alpha diversity was represented by richness, Shannon, Simpson, PD whole tree, and Chao 1 indices. Richness was presented as the observed OTU number. Three nonparametric tests (multiple response permutation procedure, MRPP; PERMANOVA; and analysis of similarity, ANOSIM) based on Bray–Curtis dissimilarity indices were performed by R software to examine the difference in soil microbial taxonomic composition among grazing regimes. A phylogenetic tree was annotated and visualized in iTOL software (Letunic and Bork, 2019).

A null modeling-based approach by Ning et al. (2019) was performed to infer soil microbial community assembly mechanism, that is, calculating the normalized stochastic ratio for estimating the relative importance of stochasticity and determinism in shaping community structure. Specifically, the null model algorithm was used to qualify the significance of the observed difference in the microbial community from random expectation. If the observed microbial community is statistically different from the null expectation, the community is regarded as largely shaped by deterministic processes. Otherwise, it is considered to be dominated by stochastic processes. Also known as normalized stochasticity ratio with 50% as the boundary point between more deterministic (<50%) and more stochastic (>50%) assembly.

To explore the relationships among bacteria, fungi, archaea, and plant species, a co-occurrence network based on the relative abundance of microbial classes and plant species dominance data was constructed using the MENA pipeline (http://ieg2.ou.edu/MENA), which performs random matrix theory to identify thresholds for constructing a highly confident ecological network (Zhou et al., 2011; Deng et al., 2012), and the default settings of MENA were adopted.

The plant community composition under RG was significantly different compared to that under GE and CG (Supplementary Table S1). Specifically, Kobresia pygmaea was the dominant species among the 13 main species; it was a palatable species for yaks and its coverage under RG was higher than that under GE and CG (Figure 1). Another priming species was Stipa capillata, which was also a palatable species. However, its coverage under GE was higher than that under RG and CG.

The soil bacterial diversity indices showed no differences among the grazing regimes (Supplementary Figure S2). The Shannon, Chao 1, richness, and PD whole tree indices of the soil fungal community were all higher under CG than under RG, while the Simpson index showed no difference among the grazing regimes (Figure 2). The soil's whole microbial community composition did not show any difference among the grazing regimes (Table 1). Stochasticity dominated soil microbial community assembly, RG, and CG significantly decreased the relative importance of stochastic processes (Figure 3).

The soil bacterial functional composition under RG and CG was similar (Figure 4A), while the soil fungal functional composition under GE and RG was similar (Figure 4B). Specifically, aerobic chemoheterotrophy and chemoheterotrophy were the main bacterial functional guilds, which include most OTU, and their relative abundances were higher under GE than those under RG and CG. Lichenized and Undefined Saprotroph were the main fungal functional guilds, while most OTU belonged to Arbuscular Mycorrhizal, with low relative abundance.

To better characterize the grazing regime effect, we identified the dominant taxa (OTUs with a relative abundance ≥ 0.5%). Although 4684 bacterial OTUs and 5241 fungal OTUs were retrieved from all the grazed plots, only 42 (0.8%) and 84 (1.6%) OTUs were identified as the dominant taxa, respectively. For all the grazed plots, these dominant bacterial OTUs were mainly affiliated with Thermoleophilia (28.6%), Actinobacteria (23.8%), and Acidobacteria (26.2%) (Figure 5A); the dominant fungal OTUs were mainly affiliated with Sordariomycetes (26.2%) and Leotiomycetes (26.2%) (Figure 5B). The majority of the dominant taxa did not change with the grazing regime. Specifically, a greater relative abundance of nine dominant fungal taxa was observed under the CG than that under the RG, including OTU 52, OTU 5178, OTU 6, OTU 3067, OTU 5287, OTU 4, OTU 82, OTU 64, and OTU 27, which all belong to Ascomycota. However, a greater relative abundance of six dominant bacterial taxa was observed under the CG than that under the GE, including OTU 46, OTU 4, OTU 19, OTU 2050, OTU 7, and OTU 23.

At the phylum level, soil bacteria were overwhelmingly dominated by Actinobacteria, Proteobacteria, and Acidobacteria; soil fungi were overwhelmingly dominated by Ascomycota, Zygomycota, and Basidiomycota (Supplementary Figure S3). The abundance of Acidobacteria and Armatimonadetes was higher under CG than that under GE, and the abundance of Elusimicrobia was higher under RG and CG than that under GE (Supplementary Figure S4).

The co-occurrence network showed 104 nodes, including 76 bacteria, 16 fungi, one archaea, and 11 plant species (Figure 6A). Node with a maximum degree and eigenvector centrality was Deltaproteobacteria, and node with maximum betweenness and stress centrality was Cytophagia. There were 288 links in the network, with 25 negative and 10 positive links between plants and bacteria; four negative links occurred between plants and fungi; one negative link occurred between plants (S. capillata) and archaea; 22 positive and 18 negative links occurred between bacteria and fungi; five positive links and one negative link occurred between bacteria and archaea; and one positive link occurred between fungi and archaea. In plant species, K. pygmaea is rarely linked with microbes, but it linked four plant specialists that were not directly linked to microbes. S. capillata had most links with microbes, 11 with bacteria, one with fungus, and one with archaea, while only two links of them were positive.

Within each group, six positive and two negative links were detected among 11 plant species, whereas 139 positive and 49 negative links were detected among 76 bacterial classes and two positive and two negative links were observed among 16 fungal classes. Overall, for bacteria, Chloroflexi had the most links of 79 among them 62 were positive; for fungi, Ascomycota had the most links of 23 among them 11 were positive.

The topological roles of the nodes identified in the network were shown as a Z-P plot in Figure 6B. The majority (95.2%) of the nodes were peripherals with most of their links inside their modules. Among these peripherals, 46.5% had no links outside their modules (i.e., Pi = 0). Approximately 4.8% of the nodes were generalists, including one module hub and four connectors. However, no network hubs (supergeneralists) were observed in the network. The module hub was Cytophagia, and the connectors were ChloroflexiSHA26, Actinomycetia, Chloroflexia, and Microbotryomycetes. These four of the five microbial generalists showed no difference in their abundances among the grazing regimes (Supplementary Figure S5).

In addition, the Mantel test showed that node connectivity was positively related to the correlation coefficient between plant aboveground biomass, belowground biomass, community coverage, and soil NO₃-N, while other factors were irrelevant (Supplementary Table S2). Thus, nodes with high connectivity tended to be more sensitive to these factors.

Our results found that the aboveground plant community and belowground soil microbial community in the grassland ecosystem showed different responses to grazing disturbance, and the fungi and bacteria in the soil microbial community also behave differently. The grazing regime significantly affected plant community composition and is mainly reflected in two species, Kobresia pygmaea and Stipa capillata, which are both palatable species for yaks. The former is the dominant species, both GE and CG decreased their coverage compared to RG. However, RG and CG significantly decreased the coverage of Stipa capillata compared to GE; Stipa capillata therefore can be regarded as a species sensitive to grazing. In contrast, the grazing regime did not alter the soil microbial taxonomic composition, implying a stable soil microbial community in the grazed grassland. For soil microbial diversity, bacterial diversity indices were not affected by grazing regime, while only the Simpson index of the fungal diversity indices had no difference among grazing regimes; these results indicated that the diversity of the bacterial community was more tolerant, while the diversity of fungi was more sensitive to grazing. Previous studies obtained similar results showing that fungal community composition was more sensitive than bacterial community composition to disturbance or management regimes (Lauber et al., 2008) and vegetation change (Dassen et al., 2017) and that there is accumulating evidence that fungi are the first consumers of belowground inputs of plant-derived C (De Deyn et al., 2011; Ballhausen and de Boer, 2016).

The microbial community assembly process is important in coupling community composition with the ecosystem function, which includes two processes: determinism and stochasticity (Zhang et al., 2019). Deterministic processes are based on the theory of ecological niche, stressing the role of environmental filtering imposed by environmental factors in the assembly of the microbial community; stochastic process occurs through ecological drift, including random birth and death and dispersal events (Langenheder and Szekely, 2011). Both deterministic and stochastic processes occur simultaneously during the community assembly and the relative importance changes along with environmental conditions (Tang et al., 2021). We found stochastic mechanisms strongly structured soil microbial community assembly under all grazing regimes. Compared to grazing exclusion, rotational grazing and continuous grazing reduced the relative importance of stochastic ratio, indicating that grazing acted as an environmental filter that increased the deterministic soil microbial community assembly.

The soil bacterial functional composition under RG and CG was similar, while the soil fungal functional composition under GE and RG was similar. From GE to RG, then to CG, the interference increased gradually, implying that the type of interference determines the functional composition of bacteria, while the intensity of interference determines the functional composition of fungi. Aerobic chemoheterotrophy and chemoheterotrophy were the main bacterial functional guilds with rich OTUs, and their relative abundance was higher under GE than that under RG and CG, probably resulting from more organic matter under GE without grazing yak's intakes.

Only a few soil microbial phyla showed differences in abundance among the grazing regimes. The abundances of Acidobacteria and Armatimonadetes were higher under CG than that under GE, implying a more acidic environment under CG (Bardhan et al., 2012). Acidobacteria is mainly considered oligotrophs and are known to utilize a variety of carbon resources, even recalcitrant carbon substrates, and contribute to litter decomposition and carbon resource conversion, benefiting plant growth (Macdonald et al., 2015; Lladó et al., 2017). Armatimonadetes is related to carbohydrate transportation and metabolism and can also utilize a wide range of carbohydrates (Lee et al., 2014). Continuous removal of plant tissues by yaks under CG reduced plant carbon entering the soil (Supplementary Figure S7), probably stimulating soil microbial carbohydrate transportation and metabolism processes, and leading to the greater abundance of Acidobacteria and Armatimonadetes under CG. Different from our results, Yang et al. (2019) observed an equal relative abundance of Acidobacteria between grazing exclusion, winter grazing, and annual grazing on the Qinghai-Tibet Plateau. Macdonald et al. (2015) found rabbit grazing had a negative impact on the Acidobacteria and thought the result related to rabbit grazing increased C-allocation belowground through root exudation. Two other studies also observed a reduced relative abundance of Acidobacteria under grazing (Xun et al., 2018; Zhang and Fu, 2021). Different results may be related to climate, herbivore type, grazing method, and intensity, thus it is difficult to draw a generalized conclusion about the grazing effects on Acidobacteria. Armatimonadetes received little attention in the studies of grassland grazing, probably because it is not a dominant taxon. At the OTU level, nine dominant fungal OTUs were higher under CG than that under RG, as CG likely benefited slow-growing plant species that produce low-quality litter (high C/N), thus favoring soil fungi; six dominant bacterial OTUs were higher under CG than that under GE, as livestock fed on aboveground plants and returned nearly half of it as feces, thus favoring bacterial growth (Xun et al., 2018). In addition, the grazing increased Acidobacteria and Armatimonadetes that are related to C metabolism, while grazing did not change other microbial phyla.

For bacteria, Chloroflexi is diverse in terms of morphology, nutrition, and metabolic pathways, mediates the biogeochemical cycle of elements such as C, N, and S (Zarzycki et al., 2009; Sorokin et al., 2012, 2014), and had the most links in the network of plant species and soil microbial classes, including 62 positive and 17 negative links. For fungi, Ascomycota, which has the highest species diversity and can break down plant residues, had the most links in the network of plant species and soil microbial classes, including 11 positive and 12 negative links. The module hub of the network was Cytophagia (phylum Bacteroidia), which is an obligate anaerobic gut microbe in animals (Eckburg et al., 2005; Winter and Bäumler, 2014), and it was the hub likely due to the accumulation of dung and urine by grazing (Barik and Murugan, 2015); the connectors were ChloroflexiSHA26, Actinomycetia, Chloroflexia, and Microbotryomycetes. The abundance of four of the five microbial generalists was not affected by the grazing regime, and the little changes may have contributed to the stability of soil microbial composition.

Across plants, bacteria, fungi, and archaea, most of the links between plants and microbes were negative, and they may compete for nutrients, such as available nitrogen, in the infertile environment (Inselsbacher et al., 2010). Plant species were more linked to bacteria than fungi. There were 22 positive and 18 negative links between bacteria and fungi, five positive and one negative links occurred between bacteria and archaea, and only one positive link was detected between fungi and archaea, indicating that bacteria, fungi, and archaea exhibited cooperative interspecies relationships more than competitive interactions. Within each network component, there were six positive and two negative links in the plant community; 73.9% of the bacterial links were positive. Two positive and two negative links were found in 16 fungi, indicating that fungi were more linked with other components than themselves and that intra-group members were more independent. Overall, soil microbes exhibited cooperative behaviors such as cross-feeding, syntrophic interactions, mutualistic interactions, and competition with plant species, similar to intraspecific cooperation and interspecific competition. Microbial communities often act in consortia to synergistically degrade complex plant-derived compounds, with some microbes utilizing metabolites or taking advantage of broken-down products of extracellular enzymes produced by other taxa (Lynd et al., 2002; Alessi et al., 2017).

In the links between plants and microbes, two remarkable plant species were K. pygmaea and S. capillata. There were few links between K. pygmaea and microbes in the network, but it served as a bridge between microbes and four plant specialists that were not directly related to microbes. S. capillata was linked with 11 bacteria, one fungus, and one archaea, and most of the links were negative. Coincidentally, K. pygmaea and S. capillata were also the two special species that had different coverage among grazing regimes. K. pygmaea was a grazing-tolerant species and the dominant species in the grazed grassland, while S. capillata was a sensitive species that was rarely distributed in the grazed grassland but sharply increased after grazing was prohibited. K. pygmaea and S. capillata were the two extremes in the plant community, the dominant species K. pygmaea acted as a bridge between the plant community and microbes, and the sensitive species S. capillata was more independent. It was found that different species play distinct roles in the community and were generally divided into three types. One type was only widely linked to microbes; the other type was plant specialists, which were not directly linked with microbes; and the third type was the bridge between microbes and plant specialists.

Among the environmental factors, plant community aboveground biomass, belowground biomass, coverage, and soil NO₃-N were closely related to the node connectivity of the network. Nodes with more links tended to be more sensitive to these factors but not to the other factors, such as plant community diversity, soil NH₄-N, available P, and bulk density. This result implied a close relationship between the soil microbial community and plant community production, in line with previous findings (Van Der Heijden et al., 2008; Yang et al., 2017; Xun et al., 2018; Chen et al., 2020). NO₃-N is much more mobile than NH₄-N in soil and can be easily lost by leaching into groundwater and surface runoff (Scherer-Lorenzen et al., 2003); thus, NO₃-N tends to be the object of competition between plants and microbes and maybe the major cause why many of the links between plants and microbes were negative.