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−2), 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 ×
104 hm2, of which the grassland area
is 117.57 × 104 hm2. The usable
grassland area is 108.53 × 104 hm2,
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{PS}=S\),
\(\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; Pi is
the relative importance value of species i , calculated by the
formula (relative abundance + relative height + relative coverage) / 3;
and Nij 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: Nij = mij /
Mi , in which Mi = \(\sum_{j}^{r}m_{\text{ij}}\), where mij 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 Bi , 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−1, TN), organic carbon
(g·kg−1, OC), total phosphorus
(g·kg−1, TP), alkali hydrolyzed nitrogen
(mg·kg−1, AHN), available phosphorous
(mg·kg−1, AP), ammonium nitrogen
(mg·kg−1, AMN), nitrate nitrogen
(mg·kg−1, 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−2, AGB), plant belowground
biomass (g·m−2, 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 2) 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).