Conceptual framework: (a) the biotic community is composed of three trophic levels—plants, rodents, and microbes. (b) Multitrophic diversity is predicted to enhance ecosystem mutifunctionality. (c) Biotic and abiotic variables have direct and indirect positive and negative effects on ecosystem multifunctionality.

2 MATERIALS AND METHODS

2.1 Experimental site

The experimental site is located in Maqin County, in the southern part of Qinghai Province in the hinterland of the Qinghai–Tibet Plateau, which is situated at the source of the Yellow River and in the core area of the Sanjiangyuan (Table S1). The altitude is 4100 m above sea level, which means the site lies within a typical continental plateau monsoon climate, with strong solar radiation (annual total radiation of 6194 MJ·m⁻²), 2493.6 h of sunshine, and no absolute frostless period (Chang et al., 2023). The mean annual temperature is approximately 0.4°C, with low temperatures and large daily temperature differences. The annual precipitation is 531.6 mm, mostly concentrated in June and September. In the last decade, the highest average temperature was 9.9°C and the lowest was −6.1°C, with rain and heat occurring at the same time, which is favorable for pasture growth (Xiang, 2023). The total area of Maqin County is 134.60 × 10⁴ hm², of which the grassland area is 117.57 × 10⁴ hm². The usable grassland area is 108.53 × 10⁴ hm², which accounts for 92.3% of the total grassland area, and the grassland type is mainly alpine meadow grass, within which the main dominant plants are Carex alatauensis , C. parvula , C. capillifolia , C. tibetikobresia , and C. hughii (Dou et al., 2023).

2.2 Experimental design

In 2023, our experiment were conducted in five distinct vegetation types on an alpine grassland (Table S1). Under each of these specific vegetation types, we selected five spatial replicates of sampling plots (15 m × 15 m), each separated by a distance of at least 2 km. In each plot, we randomly selected five quadrats (0.5 m × 0.5 m) for the vegetation survey. The aboveground biomass was collected by pruning all plants within the quadrats and putting them in envelopes. Three points were evenly selected on the diagonal of each treatment plot as soil sampling points, and 0–10 cm layers of soil were obtained from the quadrat with a soil drill (diameter: 5 cm). The samples from the same layer were mixed with the three sampling points of each treatment, and then the soil roots were sieved through a 0.5 mm mesh, before packing the soil into soil sample bags and taking them to the laboratory for calculation of the soil nutrient indexes. In addition, soil samples of 0–5 cm were obtained from the quadrats with a soil drill (diameter: 1 cm) near the soil sampling points, mixed evenly, packed into Ziplock bags, labeled, and stored at −80°C in a freezer for microbial sequencing analysis.

In the rodent survey, circular quadrats with a radius of 30 m and similar site conditions were established in different types of grassland over rodent holes. The holes were then plugged and landfilled, and the total number of holes and effective rodent holes were investigated. Because the quadrats were relatively large, to reduce the systematic error of test sampling, the circular quadrats were evenly divided into three equal parts along the radius (Figure S1). After gathering the statistics, the rodents were captured in the sample plots using the rope noose method, and the species were counted before then being released (Appendix S1).

2.3 Biotic variables

The plant diversity indices selected in this study included Margalef (PS ), Shannon–Wiener (PH’ ), Simpson (PD ), and Niche breadth (B ), to estimate the diversity of community structure, the calculation formulae of which are as follows (Wen et al., 2020): rodent diversity indices selected in this study included the number of rodent species (RS ), number of rodent holes (NRH ), number of effective holes (NEH ), and proportion of effective holes (PEH ). The microbial diversity indices selected in this study included the Chao1 index of bacteria and fungi (BC and FC ), the Shannon–Wiener index of bacteria and fungi (BH’ and FH’), and the S impson index of bacteria and fungi (BD and FD , Appendix S1). Lastly, the multitrophic diversity and species diversity (e.g., plant diversity) were calculated by averaging the standardized diversity indices measured across trophic levels (Soliveres et al., 2016; Luo et al., 2022).

\(\text{PH}^{\prime}=-\sum_{i=1}^{S}P_{i}\ln P_{i}\),

\(\text{PD}=1-\sum_{i=1}^{S}P_{i}^{2}\),

\(B_{i}=-\sum_{j}^{r}{N_{\text{ij}}\text{log}N_{\text{ij}}}\),

where S is the number of species; P_i is the relative importance value of species i , calculated by the formula (relative abundance + relative height + relative coverage) / 3; and N_ij is the frequency of resource utilization of species i in sample j , normalized by the total frequency of resource utilization of species i across all samples: N_ij = m_ij / M_i , in which M_i = \(\sum_{j}^{r}m_{\text{ij}}\), where m_ij is the dominance degree of species i on resource j , which is equivalent to the importance value of species i in the sample, and r represents the total number of samples. The larger the value of B_i , the wider the niche, and the greater the total amount of resources utilized by the species, indicating a stronger competitive ability.

2.4 Abiotic variables

Abiotic variables included altitude and soil pH, which may affect ecosystem functions or multifunctionality directly or indirectly (Peters et al., 2019).

2.5 Ecosystem functions

Fifteen ecosystem functions were measured. There were moisture content (%, MC), total nitrogen (g·kg⁻¹, TN), organic carbon (g·kg⁻¹, OC), total phosphorus (g·kg⁻¹, TP), alkali hydrolyzed nitrogen (mg·kg⁻¹, AHN), available phosphorous (mg·kg⁻¹, AP), ammonium nitrogen (mg·kg⁻¹, AMN), nitrate nitrogen (mg·kg⁻¹, NN), the ratio of organic carbon to total nitrogen (C: N), the ratio of organic carbon to total phosphorus (C: P), the ratio of total nitrogen to total phosphorus (N: P), the ratio of organic carbon to total nitrogen to total phosphorus (C: N: P), plant aboveground biomass (g·m⁻², AGB), plant belowground biomass (g·m⁻², BGB), and the ratio of plant aboveground biomass to belowground biomass (A: B). Details of the data collection for these different ecosystem functions are provided in Appendix S1.

2.6 Ecosystem multifunctionality

As mentioned above, a total of 15 important ecosystem functions were selected to characterize the ecosystem multifunctionality. Specifically, these functions are good indicators of water conservation (e.g., moisture content), soil fertility (e.g., total nitrogen), nutrition cycling and transformation (e.g., the ratio of organic carbon to total nitrogen), and community productivity (e.g., plant aboveground biomass). Since there is no standardized method for studying ecosystem multifunctionality, we used an intuitive averaging method for evaluation (Hooper and Vitousek, 1998; Maestre et al., 2012; Byrnes et al., 2014). We first normalized all ecosystem functions to a comparable range of 0–1, and then the ecosystem functions (EF) were min–max-transformed: EF = [rawEF − min (rawEF)] / [max(rawEF) − min(rawEF)] (Argens et al., 2023). The EMF and ecosystem functions (e.g., water conservation) were then calculated by averaging the normalized ecosystem functions.

2.7 Statistical analysis

In this study, the effects of different grassland types on biotic (e.g., plant, rodent, and microbial diversity indices) and abiotic variables (e.g., ecosystem functions) were assessed using linear mixed-effects models with the ‘lme4’ package. Post hoc tests were conducted employing the ‘Tukey’ method in the ’emmeans’ package (Dusza et al., 2022). Pearson correlation analysis was performed on the relationships of species diversity, keystone species, with ecosystem functions. Random Forest modeling, a machine-learning algorithm, was applied to select the most important species in different grassland types and the most important variables affecting ecosystem multifunctionality with the “rfPermute” package (Archer, 2022). The importance of genera abundance was estimated by calculating the percentage increases in the mean squared error (MSE) of variables. Higher values of MSE indicate more significant variables. We also used variation partitioning analysis to quantify the relative importance of three groups of factors as predictors of ecosystem functions and the multifunctionality via the “vegan” package in R (Oksanen et al., 2020). Additionally, we compared the effects of species diversity and multitrophic diversity on ecosystem functions and multifunctionality and evaluated the explanatory power (R ²) of single and multiple trophic levels from general linear models. All predictors and response variables were standardized to range from 0 to 1 before analysis, including (a) multitrophic diversity, (b) pH, and (c) altitude. Finally, ”vegan” and ”piecewise SEM” in the R package were used to quantify the effects of biotic and abiotic variables on ecosystem multifunctionality. To further improve the fit of the models, we created alternative models by progressively removing nonsignificant paths.

3 RESULTS

3.1 Species diversity and composition

Results showed that species diversity indices had significant differences among grassland types. Specifically, the plant diversity indices of GL5 were the highest, and the fungal diversity indices of GL4 were the highest (Figures S2 and S3). The rodent species of GL5 were the highest, with Ochotona curzoniae being the absolute dominant group of grassland rodents in the study area (Table S3). Also, there were contrasting linkages among species diversity indices. The number of rodent species was positively correlated with the number of plant species and bacterial Chao1, whereas fungal diversity indices were negatively correlated with plant and rodent diversity indices (Figure S4). Dominant species were composed of different biotic communities among grassland types. Grassland types were classified by the dominant plants as C. moorcroftii , C. alatauensis , Poa annua , Ligularia virgaurea, and C. alatauensis withDasiphora fruticosa (Figure 1a). The keystone species in the plant community had different niche breadths among sample sites (Figure 1b). The dominant species of microbial and rodent communities changed with the plant community (Figures S5 and S6).