This research was conducted at the town of Xihai, Haiyan Country (36°44′-37°39′ N, 100°23′-101°20′ E), Qinghai Province, China, situated in the Qinghai Lake basin. With an average elevation of 3,100 m, this region has a typical plateau continental climate. The mean annual temperature and precipitation were 1.4°C and 330–370 mm, respectively. Specifically, the non-growing season in this region is from October to April and is long and cold, with the average temperature of the coldest month being −24.8°C. The growing season is from May to September, with an average temperature of 12.5°C in the hottest months. Moreover, the precipitation also mainly concentrated in the growing season. The soil is a clay loam and the plant community is dominated by Kobresia humilis (C. A. Mey. ex Trautv.) Sergiev, Leymus secalinus (Georgi)Tzvel., Elymus nutans Griseb. Nachr. Ges. Wiss. Gott., Carex aridula V. Krecz. and Potentilla acaulis L. Sp. Pl.

We set up a grazing experiment with local livestock yak and Tibetan sheep started in 2014 with a completely randomized block design. It had been proved that moderate intensity grazing can maintain higher aboveground biomass and plant diversity (Gao and Carmel, 2020; Liu et al., 2023). Based on moderate intensity grazing, which was represented as the feed intake of livestock for above-ground biomass is about 50–55%, five grazing treatments including single grazing and mixed grazing were set up, which were: single grazing for yaks (YG), single grazing for Tibetan sheep (SG), mixed grazing for yaks and Tibetan sheep at the ratio of 1:2 (MG12), 1:4 (MG14), 1:6 (MG16), respectively, and a control treatment (CK, grazing exclusion) in three blocks, there were 18 plots in total. By adjusting the area of each grazing plot according to the livestock numbers, the grazing intensity of each plot was consistent. A detailed description of the experimental site and design can be found in our previous studies (Yang et al., 2019).

In August 2019, four sites were randomly selected in each plot for plant community investigation and soil sampling. Aboveground biomass were harvested from a 0.5 m × 0.5 m quadrat in each sampling site. After plant community investigation, soil sampling was conducted in the same quadrat (Yang et al., 2019). Three random soil cores (0–10 cm, 10–20 cm and 20–30 cm) were taken in each quadrats. The samples from each quadrats were collected and mixed to form a composite sample by different depth. After careful removal of the surface litter, earthworms or bird droppings, each sample was mixed and divided into two parts. One part was air-dried for analysis of the soil physicochemical properties, and the other was stored in a freezer at −80°C and then sieved through a 2-mm sieve before analysis of the microbial communities.

The soil total carbon (TC), soil total nitrogen (TN), soil total phosphorus (TP), ammonium nitrogen (AN) and pH were measured following method of Huang (Huang et al., 2019).

Microbial community genomic DNA was extracted using the E.Z.N.A.® soil DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to manufacturer’s instructions. The DNA extract was checked on 1% agarose gel, and DNA concentration and purity were determined with NanoDrop 2000 UV–vis spectrophotometer (Thermo Scientific, Wilmington, United States).

The copies of nifH gene were quantified using ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The quantification was conducted with universal primer sets for nifH, PolF: TGC GAY CCS AAR GCB GAC TC; PolR: ATS GCC ATC ATY TCR CCG GA (Poly et al., 2001; Che et al., 2018). The 10 μL reaction systems contained: 4.4 μL of SYBR Green Mix, 0.3 μL of forward primer (20 μmol L⁻¹), 0.3 μL of reverse primer (20 μmol L⁻¹) and 5 μL of template DNA. The standard curve was constructed using plasmids inserted with the nifH gene fragment. The PCR runs started with an initial denaturation at 95°C for 10 min, followed by 40 cycles of 10 s at 95°C, 34 s at 60°C, 15 s at 95°C, 60 s at 60°C and 1 cycle of 30 s at 95°C, 15 s at 60°C. The specificities of PCR products were checked by melting curve analysis.

The PCR amplifications of nifH gene was conducted with barcoded universal primer set (F: 5’-AAAGGYGGWATCGGYAARTCCACCAC-3′; R: 5’-TTGTTSGCSGCRTACATSGCCATCAT-3′). Purified amplicons were pooled in equimolar and paired-end sequenced on an Illumina MiSeq PE300 platform (Illumina, San Diego, United States) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw nifH gene sequencing reads were demultiplexed, quality-filtered by fastp version 0.20.0 and merged by FLASH version 1.2.7 (Magoč and Salzberg, 2011). Operational taxonomic units (OTUs) with 97% similarity cutoff were clustered using UPARSE version 7.1, and chimeric sequences were identified and removed (Edgar, 2013). The taxonomy of each OTU representative sequence was analyzed by RDP Classifier version 2.2 against the 16S rRNA database using confidence threshold of 0.7 (Wang et al., 2007).

In total, 898,086 nifH gene sequences were obtained for the 54 soil samples, and the total sequence number was reduced to 859,258 after the quality filtering. The qualified sequence number of each sample ranged from 9,808 to 24,478, and thus the sequence number of each sample was rarefied to 9,808 for further analysis. With an identify cutoff of 97%, we obtained 2,421 OTUs. More details of the taxonomic assignment for nifH OTUs was conducted similarly as described by (Che et al., 2018). The data were analyzed on the online platform of Majorbio Cloud Platform.¹

One-way ANOVA was used to analyze the effects of grazing on soil properties, abundance of nifH genes and OTU numbers, followed by Duncan’s new multiple range. Data were tested for normality and homogeneity of variance before analysis using the Shapiro–Wilk test and Levene’s Test, respectively. Two-way ANOVA was performed testing the main and interactive effects of grazing and soil depth on soil properties. Mantel tests were performed to assess correlations between diazotrophic taxa and experimental factors. These analyses were performed with the IBM Statistical Package, SPSS version 25.0 (IBM, Armonk, NY, United States), and the histogram plots were plotted by Sigma plot software. Linear regression analyses and structural equation model (SEM) were performed to test the relationships between the OTU numbers and environmental factors using the ‘ggplot2’ and ‘piecewiseSEM’packages in R version 4.1.0.

Soil properties are shown in Table 1. A two-way ANOVA showed that the grazing treatments significantly affected all soil property parameters (p<0.01), and soil depth effect (0–10 cm, 10–20 cm and 20–30 cm) also influenced these parameters (p<0.01). Grazing treatments and soil depth had an interaction effects on soil ammonium nitrogen content (AN, p<0.01) (Supplementary Table S1). Compared with the no grazing CK, the YG and MG14 treatments significantly reduced the soil pH at 0–10 cm and 10–20 cm, respectively (p<0.05). Soil total carbon (TC) decreased significantly with the increase of soil layer in MG16 treatment (p<0.05). Soil total nitrogen (TN) under grazing treatment did not change significantly in soil layers 0–10 (p>0.05), but with the increase of soil layer, soil total nitrogen decreased significantly (p<0.05). In soil layers 10–20 and 20–30, soil total nitrogen under SG treatment was the highest, and soil total nitrogen under MG16 treatment was the lowest (Table 1). The grazing treatments significantly increased the soil AN (p<0.05), especially in the SG and YG treatments, and the AN are higher in the grazing alone treatment than that in the mixed grazing treatment. Moreover, compared to the CK, the grazing treatments also significantly increased the soil total phosphorus (TP, p<0.05), especially in the SG and YG treatment in soil depth of 0–10 cm (Table 1). In addition, the N/P ratio was lowest in SG treatment, while highest in MG12 treatment (p<0.05).

Compared with the CK, both grazing treatments significant decreased the nifH gene copy number in all soil depth (p<0.05) (Figure 1). The copy number of nifH gene in the CK treatment decreased significantly with the increase of soil layer (p<0.05), whereas no significant change was observed among the other grazing treatment (p>0.05). In the 0–10 cm from the soil surface, the average copy numbers of the six treatments ranged from 0.73 × 10⁵ to 9.19 × 10⁵ gene copies per g dry soil (Figure 1).

The CK and MG12 treatments had similar OTU numbers, and the YG and MG16 treatments had the lower the number of OTUs (p<0.05) (Figure 2). There were no significant differences among CK, SG, MG12 and MG14 (p>0.05) (Figure 2). OTUs were taxonomically classified into different order across all of the treatments. Figure 3 shows the relative abundances of the most abundant order in soils of the six grazing treatments. The order Rhodospirillales was the most abundant, containing 58–81% of the total nifH gene sequences in all of the soil samples. Different grazing treatments changed the relative abundances of these main order (Figure 3). Compared with CK, SG, YG and MG12 had similar dominant communities, while the relative abundances of Rhodospirillales decreased and Nostocales increased in the MG14 and MG16 treatments.

The Mantel test revealed that the diazotrophic community structures in soil surface were closely correlated with multiple soil variables, such as pH, TC, AN and TP in the 0–10 cm of soil (p<0.05) (Table 2). then the nifH gene copies only had a negative correlation with TP in 10–20 cm and 20–30 cm of soil layers (p<0.05). whereas had no significant correlation with aboveground biomass and root biomass (p>0.05) (Table 2). However, specifically, Rhizobiales and Bacillales are closely related to leguminous biomass, and they are related to soil TC and TP, respectively. In addition, Nostocales is related to forb biomass, and Opitutales is related to aboveground biomass (Figure 4). Based on the linear regression analysis, there are negative correlations between OTUs number, AN and TP in the 0–10 cm of soil (p<0.05) (Figure 5).

The effects of soil properties and plant functional group on the structures of diazotrophic communities were further analyzed using the redundancy analysis, based on the selected soil variables, biomass of plant functional group and OTU composition. These soil variables explained 47.5% of the variation, and the first two axes explained 30.43 and 17.07% of the total variation. According to the vectors, the diazotrophic communities of both CK and grazing treatments were associated with higher TP and AN values (Figure 6). In addition, the diazotrophic communities were also associated with leguminous biomass and forbs biomass. These biomass of plant functional group explained 58.53% of the variation, and the first two axes explained 41.9 and 16.63% of the total variation (Figure 7).

SEM were used to identify the key drivers on diazotrophic communities. Path analysis indicated that grazing had an impact on diazotrophs through the direct effects on leguminous biomass, AN and TP content. Grazing reduced the biomass of legumes (Supplementary Figure S1), but increased AN and TP content, and legumes were positively correlated with nitrogen-fixing bacteria, so nitrogen-fixing bacteria also decreased under grazing treatment (Figure 8).

Many microbially driven processes in soils can be impacted by land management practices and changes in the plant community (Lindsay et al., 2010; Meyer et al., 2013). Previous studies had described the diversity, composition and structure of the diazotroph community under different systems (Song et al., 2019; Xiao et al., 2020a; Barros et al., 2021). However, little information is available about the changes of diazotrophic community in Qinghai-Tibetan Plateau grassland soils with different grazing patterns. Thus, this study evaluated the influence of grazing patterns on diazotrophic communities in Qinghai-Tibetan Plateau grassland under moderate intensity grazing. The decreased of nifH gene copies under all grazing treatments were consistent with previous reports (Yan et al., 2021). In the present study, compared with CK, both singly and mixed grazing treatment reduced the nifH gene copies. This is consistent with Lindsay’s findings, they found that there was a trend for higher abundance of nifH in sites where grazing was excluded (Lindsay et al., 2010). The significant negative relationship between nifH gene abundance and soil physicochemistry (TC, AN and TP) under grazing treatment, suggested that frequent livestock grazing could inhibit dinitrogenase reductase and lead to a reduction in the biological capacity for nitrogen fixation. A large proportion of the nitrogen in dung and urine from livestock can provide high nitrogen inputs in various forms, therefore can reduce the ecosystem dependence on free-living nitrogen fixing organisms (Singh et al., 2011).

With regard to the effects of grazing patterns on the diazotrophic community composition, a previous study showed that intensive grazing may increased the relative abundances of Proteobacteria, Bacteroidetes and Firmicutes (Zhang et al., 2023), whereas another study found that grazing exclusion increased the relative abundance of these bacteria (Wang Z. et al., 2021). In our study, YG and MG16 significant reduced the OTUs number of diazotrophs, and the dominant diazotrophic communities belong to Proteobacteria, Firmicutes, Rhizobiales. The high proportions of Proteobacteria in CK and MG12 treatment indicate the mix grazing patterns with ratios of yak to sheep as 1:2 is closed to the natural grassland which without grazing in regulating soil microbial community. Meanwhile, the dominant bacteria had significant correlation with legume biomass suggested the abundance of nitrogen-fixing bacteria in soil is influenced by the foraging strategies of different livestock. Among the Proteobacteria and Firmicutes phylum, taxa such as Rhizobiales and Bacillales species are well known for forming root nodules with legumes and conducting symbiotic nitrogen fixation, and varying significantly in abundance among the types of land use (Yan et al., 2021).

Both of the nifH gene copies and OTUs had negative correlation with soil TP. Although diazotrophic communities were sensitive to phosphorus, a study in karst ecosystems reported that soil P availability played an important role in regulating N₂ fixation by increasing diazotroph diversity (Xiao et al., 2020b), it was inconsistent with our results. Our preliminary work found that the response of soil TP showed a strong dependence on livestock species (Yang et al., 2019), the responses of TP to grazing were significantly positively related to the proportion of sheep in the mixed livestock group. This suggests that changes in N₂-fixing communities may also be related to the proportion of sheep in mix grazing. In addition, grazing can also increase soil AN content due to the feces and urine of small herbivores have a relative high N content (Liu et al., 2023). Previous studies indicated that high available nitrogen content had a negative effect on soil diazotrophs as the nitrogenase enzyme is sensitive to ammonia (Che et al., 2018; Han et al., 2019). Moreover, the increased content of TP and AN enriched soil available nutrients and reduced the dependence of plants on nitrogen-fixing bacteria. Consequently, this could explain why nifH gene copies decreased in grazing treatments. Livestock grazing in grassy woodlands could potentially alter the nitrogen cycle (Lindsay et al., 2010), and grazing has been found to enhance the activity of soil nitrifying and denitrifying bacteria (Patra et al., 2005). However, herbivory can impose strong limits on diazotroph plant host abundance (Ritchie et al., 1998). Meanwhile, pasture management may affect the bacterial community through change the botanical composition and pasture productivity (Wakelin et al., 2009), and the grass roots can harbor abundant endophytic N₂-fixing microbes (Ritchie and Raina, 2016).

The type of land management and land use intensity has been identified as a major driver for microbial performance in soil (Meyer et al., 2013). In this study, grazing treatment had a significantly effects on soil AN. The soil AN and nifH were negatively correlated in different grazing patterns, suggesting that aboveground herbivory reduces the capacity for belowground nitrogen fixation (Lindsay et al., 2010). The composition and activities of diazotrophic communities are essential for the functioning of the soil nutrient cycles (Hsu and Buckley, 2009). Due to plants and soils can form a complex mutual feedback relationship, the root biomass and aboveground biomass of legume and herb also play an important role in regulating the diazotrophic communities. Plant communities and soil variables can potentially produce an evident interaction in their influence on the soil organic carbon and nitrogen contents (Yuan and Jiang, 2021). As the most important and common land use of grasslands, livestock grazing and especially over-grazing can alter the composition and productivity of the plant community and can accelerate the loss of soil fertility and the depletion of other resources (Wang B. et al., 2021). In the present study, the herbivore assemblage did not change the diversity of the N₂-fixing communities, but significantly reduced their abundance. Moreover, MG12 and CK had similar OTUs of diazotrophs, which revealed that the composition of N₂-fixing microorganisms could be maximally maintained when yak and Tibetan sheep were 1:2. Based on SEM, the effects of grazing on diazotrophs were achieved by changing leguminous biomass, soil AN and soil TP, respectively.