The Combination of Plant Diversity and Soil Microbial Diversity Directly and Actively Drives the Multifunctionality of Grassland Ecosystems in the Middle Part of the Northern Slopes of the Tian Shan under Grazing Disturbance

Abstract

It is well known that biodiversity and ecosystem multifunctionality (EMF) guarantee the well-being of human society. Most studies have focused on the relationship between biodiversity and ecosystem function, and less is known about the individual and combined effects of above- and below-ground biodiversity on ecosystem multifunctionality under grazing disturbance. The aim of our study was to investigate the relationship between plant and soil microbial (bacterial and fungal) diversity and ecosystem multifunctionality under grazing disturbance by using multiple methods to assess ecosystem multifunctionality. We conducted experiments in desert grasslands on the northern slopes of the Tian Shan Mountains and compared the relationship between ecosystem multifunctionality and biodiversity assessed by different methods under light grazing and heavy grazing. Our results showed that at the heavy grazing level, ecosystem multifunctionality calculated by the mean method and plant diversity, soil fungal diversity, soil bacterial diversity and soil fertility calculated by the single function method showed a significant decrease (p < 0.05), but grass productivity was significantly increased (p < 0.05). Among them, ecosystem multifunctionality, soil carbon storage function and soil fertility all showed significant positive correlations with plant diversity and soil microbial diversity (p < 0.05). We calculated that ecosystem multifunctionality also essentially showed positive correlation with plant diversity and soil microbial diversity using the multi-threshold method, and the effect curve was approximately a single-peaked curve, first increasing and then decreasing. Finally, we used plant diversity, soil fungal diversity and soil bacterial diversity under grazing disturbance as biotic factors and soil pH as an abiotic factor to construct structural equation models, and we found that grazing can have direct effects on ecosystem multifunctionality and indirect effects on ecosystem multifunctionality through above- and below-ground biodiversity. Our study emphasizes the importance of the combination of above- and below-ground biodiversity in maintaining the multifunctionality of desert grassland ecosystems on the northern slopes of the Tian Shan Mountains. A moderate reduction in grazing intensity can better conserve biodiversity and improve ecosystem multifunctionality, and it is a feasible strategy to maintain sustainable management of desert grasslands.

1. Introduction

Grasslands are important terrestrial ecosystems across the Earth’s surface and are one of the most significant parts of China’s terrestrial ecosystems [1,2]. Grassland ecosystems provide a variety of services for humanity, including material, environmental and cultural functions. They are the support systems that support and sustain life on Earth and are the basis for people’s survival and development [3]. In recent years, the loss of biodiversity has accelerated due to climate change and anthropogenic disturbances, greatly affecting ecosystem functional services [4,5]. However, ecosystem functions are essential for human survival and ensure the well-being and sustainable development of human societies [6]. Scholars have begun to find that focusing on a single function alone can easily lead to overlooking the trade-offs between ecosystem functions [7]. As a result, scholars have gradually begun to consider multiple ecosystem functions together [5] and proposed the concept of ecosystem multifunctionality (hereafter EMF). EMF considers how biodiversity affects multiple ecosystem functions simultaneously, whereby policy makers and managers can adopt measures to maintain ecosystem functions at the most appropriate level to maximize resource use and ecosystem sustainability.

EMF refers to the ability of an ecosystem to provide multiple services and functions [8,9], including net primary material production, carbon storage and sinks, climate regulation, soil and water conservation and wind and sand control, among other ecosystem functions. Biodiversity is the synthesis of ecological complexes formed by organisms and their environments and the various ecological processes associated with them, including animals, plants, microorganisms and the genes they possess and the complex ecosystems they form with their habitats [10]. A number of studies have shown that loss of biodiversity has multiple effects on ecosystem functions [11,12], such as decreased productivity, decreased soil nutrients, imbalance in nutrient cycling functions, decreased soil carbon accumulation and decreased pollen transfer capacity [13,14,15]. Currently, exploring the relationship between biodiversity and ecosystem multifunctionality (hereafter BEMF) is one of the research hotspots in ecology [16]. In fact, the earliest research on BEMF can be found in the early 19th century when George Sinclair explored the relationship between species diversity and biomass using grassland experiments in England. The results showed that biomass increased with species diversity, and his results were later cited by Darwin [17]. However, this is only a study of the relationship between a single ecosystem function and biodiversity, and biodiversity impacts multiple ecosystem functions simultaneously. For example, an increase in plant diversity can lead to an increase in productivity and system stability [18]. Since 2007, when Hector and Bagchi [8] first quantified the effects of plant species diversity on multiple ecosystem processes simultaneously, there has been a growing body of research on BEMF. However, relatively few studies have been conducted on BEMF under disturbance [19]. Disturbance is a phenomenon that exists everywhere in nature and affects the structure and function of ecosystems [20]. In BEMF studies, disturbance is an important driver that can directly and indirectly affect biodiversity and ecosystem multifunctionality, and BEMF relationships are highly context dependent [21]. Natural disturbances such as temperature changes increase the homogeneity of plant communities as well as the abundance of bacterial communities, thus increasing EMF [22]. In addition, alterations to the frequency of precipitation change the relative abundance of cyanobacterial communities, thus decreasing EMF [23]. Anthropogenic disturbances such as long-term fertilization can increase soil microbial community diversity and thus increase EMF [24]. Mowing, enclosure and grazing are common uses of grassland ecosystems, and a study found that grazing significantly reduced plant diversity, and thus EMF, compared with mowing and enclosure [25]. Thus, although biodiversity has a regulatory effect on EMF under disturbance conditions, the current regulatory mechanism of biodiversity on EMF is not perfect due to the complexity of both ecosystems and external disturbances. Therefore, a more detailed study on the regulation mechanism of biodiversity on EMF under disturbance conditions is needed, which will be an important reference for the conservation and management of ecosystems.

In summary, there are many studies on the relationship between disturbance and BEMF [26,27], but most of them mainly focus on the regulatory pathways of plant diversity on ecosystem multifunctionality [28] and do not involve subsurface biodiversity. In addition, studies on BEMF have mainly focused on alpine meadow ecosystems and forest ecosystems [29,30], and relatively few studies have been conducted on desert grasslands in the Tian Shan Mountains of Xinjiang, China. However, the desert grasslands in the middle part of the northern slope of the Tian Shan Mountains are important spring and autumn grasslands in Xinjiang. Desert grasslands are relatively widely distributed and an important type of grassland in Xinjiang, and they play an important role in maintaining the stability of grassland ecosystems in Xinjiang. Finally, there are various methods for evaluating ecosystem multifunctionality, and different methods have their own advantages and disadvantages, whereas most of the methods for evaluating BEMF relationships focus on the mean method and single function method. Relatively few studies have used multiple methods to evaluate ecosystem multifunctionality in an integrated manner. Therefore, in order to bridge the above knowledge gap, this paper synthesizes the shortcomings of previous studies and adopts an experiment with controlled grazing intensity. Above-ground and below-ground biodiversity were considered together under grazing disturbance, and the relationship between biodiversity and ecosystem multifunctionality was explored by using different methods for evaluating ecosystem multifunctionality (to avoid methodological limitations). The aim of this study was to test the hypothesis that the combination of above- and below-ground biodiversity under grazing disturbance could better maintain ecosystem function and multifunctionality.

2. Materials and Methods

2.1. Overview of the Study Area

The area of this study is located in the middle of the northern slope of the Tian Shan Mountains in the Asheli Kazakh National Township, Changji City, Xinjiang Uygur Autonomous Region, People’s Republic of China, with an altitude of 2200–2400 m [31] and latitude and longitude of 43°59′–43°60′ N, 86°67′–86°71′ E. The climate of the study area is dry, with more precipitation only in the summer at 150–200 mm per year. The soil type of the mountain is gray-brown desert soil; the grass type is Ili sericea desert grass; the most common plant species is Seriphidium transiliense, and the associated plant species are Tulipa iliensis, Ceratocephalus orthoceras, Malcolmia africana, Ceratocarpus arenarius, Salsola collina and Chenopodium album.

2.2. Experimental Design

The experiment was established in 2015 to study the effects of different management regimes on grassland ecosystems. The grazing method in this study area is free grazing, and the grazing time is from April to May each year. The grazing livestock are a mixture of sheep, cattle and horses, with sheep as the main species and cattle and horses as supplementary species. Among them, the sheep breed is mainly Xinjiang fine wool sheep.

We selected six herders’ grazing sample plots under the guidance of experts as the replication group for the experiment, and each grazing sample plot needed to be as similar as possible in terms of topography, slope and grazing history. The area of each grazing sample plot was roughly 0.33 km². Regarding the classification of grazing intensity in the study area, on the one hand, we classified the grazing sample sites according to the composition of the grassland plant community and the dominant types of plant species [32]. On the other hand, the sites were classified according to their distance from the herders’ settlements, and the sites 500 m away from the herders’ settlements were classified as heavily grazed (HG). Since the grazing method of herders in this study area is settled free grazing, the sample plots are closer to the herders’ residences, and the livestock are more likely to gnaw the grassland frequently. Combined with the changes in grassland plant community composition and dominant plant species, the plant communities in the heavily grazed sample sites were seriously degraded; the dominant plant species was Salsola collina, and the relative coverage of the grassland was less than 70%. The sample site 2000 m away from the herders’ residence was classified as lightly grazed (LG). This sample site is far away from the herders’ settlement, and livestock are restricted by the herders and do not appear frequently in this sample site. In addition, the plant community structure of the lightly grazed sample site in previous studies was a two-layer canopy structure with a relative cover greater than 95% [33,34], whereas the study area in this paper is a desert grassland with few plant species and low biomass. The dominant plant species in the upper layer of the lightly grazed (LG) sample site was Seriphidium transiliense, and the lower layer was Festuca ovina, whose plant individuals were relatively small, and the relative cover of grasses in this sample site was greater than 95%. In addition, the grass type of each grazing sample site in this study was desert grassland, and the differences in the dominant species of plant communities in the two treatments were caused by the differences in grazing intensity. The experimental design is shown in Figure 1.

2.3. Sample Collection and Index Determination

We set up a parallel sampling strip to create grassland sampling points in each of the two grazing intensity treatments in the above experimental area and set up 71 m × 1 m sample squares (spaced about 50 m apart) in each parallel strip and conducted vegetation surveys in August 2020. First, we made sure to avoid livestock manure areas when selecting sample squares to avoid experimental errors. Second, we needed to record the plant species composition in each sample square, measure the natural height of plants, estimate the cover of the sample plant community by visual inspection and record the number of each plant by counting. Finally, we cut all plants in the sample squares flush with the ground, weighed them and placed them in marked envelopes and brought them back to the laboratory.

Meanwhile, based on the grassland vegetation survey, we used the soil auger method (3 replicates per sample square) to collect soil from the 0–10 cm soil layer. Then, we mixed soil samples from 7 samples of the same strip to make 1 soil sample and sieved the soil samples using a 2 mm sieve to remove plant roots and gravel from the soil. Finally, we divided the obtained soil samples into 2 portions and placed 1 soil sample in a sterile centrifuge tube in an ice box, collected and stored at −20 °C in the laboratory for characterizing the soil microbial community. The other one was placed in a sealed bag and labeled, taken back to the laboratory after the sampling process was completed, and then dried naturally and used to determine the soil’s physicochemical properties.

2.4. Determination of Soil Physical and Chemical Property Indicators

Regarding the methods for measuring soil physical and chemical property indicators, this study referred to Bao’s “Soil Agrochemical Analysis” [35]. We used a potentiometric method to determine the soil pH with a water–soil ratio of 5:1. Additionally, we determined the soil organic carbon content by a potassium dichromate concentrated sulfuric acid external heating method; total nitrogen content by the Kjeldahl method; total phosphorus by the sodium hydroxide alkali fusion–molybdenum antimony anti-colorimetric method; total potassium by the sodium hydroxide alkali fusion–flame photometric method; fast-acting nitrogen by the sodium hydroxide alkali diffusion method; fast-acting phosphorus by the sodium bicarbonate leaching–molybdenum antimony anti-colorimetric method; and fast-acting potassium by the ammonium acetate leaching–flame photometric method.

2.5. Soil Bacterial DNA Extraction and Sequencing

We extracted the total soil microbial DNA according to the instructions of the MP-Fast DNATMSpin kit for soil (Omega Bio-tek, Norcross, GA, USA), and after the genomic DNA extraction, the extracted genomic DNA was detected using 1% agarose gel electrophoresis. Immediately after PCR amplification of the corresponding regions of the bacterial 16 s rRNA gene and fungal ITS rRNA gene using their respective primers, the reaction conditions were 95 °C pre-denaturation for 3 min, 27 cycles (95 °C denaturation for 30 s, 55 °C annealing for 30 s, 72 °C extension for 45 s) and finally 72 °C extension for 10 min. PCR products were mixed and detected by 2% agarose gel electrophoresis, and the PCR products were detected and quantified by the QuantiFluorTM-ST blue fluorescence quantification system (Promega). Finally, the purified amplified fragments were constructed into a library of PE250 according to the Illumina MiSeq platform (Illumina, San Diego, CA, USA) guidelines and sequenced using Illumina’s MiSeq PE250 platform. The library construction and sequencing tasks for our study were performed by Shanghai Meiji Biomedical Technology Co. (Shanghai, China).

2.6. Quantification and Evaluation of Ecosystem Multifunctionality

We referred to previous scholars’ studies on ecosystem multifunctionality [36,37,38], while considering the regulating, supplying and supporting service functions of grassland ecosystems. Finally, we selected 6 categories of functions (20 variable indicators), including soil fungal diversity (hereafter SFD), soil bacterial diversity (hereafter SBD), plant diversity (hereafter PD), soil carbon storage (hereafter SCS), soil fertility (hereafter SF) and grassland productivity (hereafter GP), to establish the ecosystem multifunctionality evaluation system in our paper (

Table 1. Grassland ecosystem multifunctionality evaluation system.

Type	Ecosystem Functioning	Ecosystem Function Indicators
Supply function	Grassland productivity (GP)	Above-ground biomass
Regulating function	Soil fungal diversity (SFD)	Fungal Shannon index; Fungal Simpson index; Fungal ACE index; Fungal Chaol index
	Soil bacterial diversity (SBD)	Bacterial Shannon index; Bacterial Simpson index; Bacterial ACE index; Bacterial Chaol index
	Plant diversity (PD)	Plant Shannon index; Plant Simpson index; Plant Margalef index; Plant Pielou index
Support function	Soil carbon storage (SCS)	Soil organic carbon;
	Soil fertility (SF)	Soil total nitrogen; Soil available nitrogen; Soil total phosphorus; Soil available phosphorus; Soil total potassium; Soil available potassium

).

2.7. Statistical Analyses

(1)

We calculated the plant diversity index based on the results of the plant community survey:

$$\mathrm{Shannon}\text{–}\mathrm{Wiener}\mathrm{diversity}\mathrm{index}:\mathrm{H}=-\sum {\mathit{P}}_{\mathit{i}}{\ln \mathit{P}}_{\mathit{i}}{(\mathit{P}}_{\mathit{i}}{=\mathit{N}}_{\mathit{i}}/\mathit{N})$$

$$\mathrm{Simpson}\mathrm{ecological}\mathrm{advantage}\mathrm{index}\text{:} \mathit{D}=1-\sum {\mathit{P}}_{\mathit{i}}^2$$

$$\mathrm{Pielou}\mathrm{uniformity}\mathrm{index}:\mathit{E}=\mathit{H}/\ln \mathit{S}$$

$$\mathrm{Margalef}\mathrm{richness}\mathrm{index}:\mathrm{Ma}=(\mathit{S}-1)/\ln \mathit{N}$$

In the formula, P_i is the ratio of the number of individuals of the ith species to the total number of individuals of all species in the community, N is the total number of individuals of all species in the sample, and S is the number of species in the sample. In addition, we generated above-ground biodiversity indices (hereafter “plant biodiversity”) for each sample using the averaging method. This was achieved by averaging the standardized values of the Shannon, Simpson, Margale and Pielou indices (mean = 0, SD = 1).

(2)

Processing of soil microbial data

We first concatenated the data using Flash software and then used Trimmomatic software for quality control to obtain a valid sequence. Secondly, the received sequence was based on the similarity of 97% as the standard, and the corresponding α diversity index was obtained through the Mothur software. Finally, the data were compared to those in the Silva (Release128) database, the species of different sequences was classified and marked, and then the community composition of various samples was counted to analyze the diversity of species in different samples. In addition, we evaluated the Shannon, Simpson, Chao1 and ACE indices of soil fungi and soil bacteria using the averaging method to generate a single index of soil fungal or soil bacterial diversity. This was achieved by averaging the standardized values (mean = 0, SD = 1) of soil bacterial and soil fungal diversity, as has been performed elsewhere, to generate a single value of subsurface biodiversity [15,39].

(3)

Calculation of ecosystem multifunctionality index

Currently, there is no universally accepted method to measure EMF. Different methods have their own advantages and disadvantages, and the measured results of each are different. Therefore, in order to avoid the errors caused by experimental methods, this study combines multiple methods to measure ecosystem multifunctionality. We used the single function method to calculate the individual functions of the ecosystem, but this method cannot reflect the overall functions of the ecosystem. However, recent findings suggest that focusing only on single functions and interspecific relationships underestimates the importance of biodiversity to ecosystem function, and that species-rich communities are more likely to maintain higher levels of ecosystem multifunctionality. Thus, we then used the mean value method to quantify ecosystem multifunctionality and to quantify ecosystem function as a whole [40]. The disadvantages of the averaging method are that the weights between functions are less easy to measure, and the default in the calculation is that the decrease in one function can be compensated by the increase in another function, and it is impossible to determine whether the EMF is averaged from several functions with similar highs and lows or from several functions with unequal highs and lows [41,42,43]. To compensate for the shortcomings of the mean method, we finally used the multi-threshold method to assess ecosystem multifunctionality. The multi-threshold method does not take into account the existence of trade-offs and interactions between functions, and it is also good for obtaining the number of functions that reach the threshold. It has a wide range of applicability [44], and it is good for analyzing the relationship between biodiversity and ecosystem multifunctionality under grazing disturbance over the entire threshold range (0–100%). The specific methodology is as follows.

We first standardized the 20 indicators in Table 1, using the min–max standardization method [45] to put the indicators in the same order of magnitude. The top 5% of observations for each variable indicator were averaged as the maximum value of the indicator, and the bottom 5% of observations were averaged as the minimum value of the indicator.

$$\mathrm{Min-max}\mathrm{standardized}\mathrm{calculation}:{\mathit{f}}_{\mathit{ij}}=\left({\mathit{x}}_{\mathit{ij}}-{\mathit{min}}_{\mathit{ij}}\right)/\left({\mathit{max}}_{\mathit{ij}}-{\mathit{min}}_{\mathit{ij}}\right)$$

where f_ij is the normalized value of the type j ecosystem function variable of plot i, and x_ij is the actual measurement value of the type j ecosystem function variable of plot i. min_ij is the minimum value of the j ecosystem function variable across all plots of the same factor, and max_ij is the maximum value of the j ecosystem function variable across all plots of the same factor.

$$\mathrm{The}\mathrm{single}\mathrm{function}\mathrm{method}\mathit{F}:{\mathit{F}}_{\mathit{ij}}={{\displaystyle \sum }}_{\mathit{j}}^{\mathit{n}}{\mathit{f}}_{\mathit{ij}}/\mathit{n}$$

$$\mathrm{The}\mathrm{average}\mathrm{method}\mathrm{EMF}:{\mathrm{EMF}}_{\mathit{i}}=1/\mathit{N}{{\displaystyle \sum }}_1^{\mathit{N}}{\mathit{f}}_{\mathit{ij}}$$

$$\mathrm{The}\mathrm{multi-threshold}\mathrm{method}\mathrm{EMF}:{\mathrm{EMF}}_{\mathit{i}}={{\displaystyle \sum }}_1^{\mathit{N}}\left({\mathit{f}}_{\mathit{ij}}{\text{>}\mathit{t}}_{\mathit{i}}\right)$$

where F_ij is the functional index of the jth function of sample site i; n is the number of indicators of ecosystem variables contained in this function; EMF_i is the ecosystem multifunctionality index of sample site i, calculated based on the standardized average of indicators of all variables in the sample site; N is the number of all ecosystem functions contained in sample site i; and t_i is the threshold value.

(4)

Data analysis

We used Excel 2019 software (Microsoft Office) to organize and analyze data, and we used the independent-samples t test with SPSS software (version 26.0) to analyze significant differences in ecosystem function and multifunctional indices of grassland ecosystems under different grazing disturbances. The R multifunc package was used to calculate multi-threshold and averaging methods, and a general linear model was developed using R software (version 4.1.1) (to analyze the regression relationship between ecosystem function and multifunctionality and biodiversity. Finally, we used structural equation modeling to analyze how grazing disturbances can indirectly change EMF by altering biodiversity.

3. Results

3.1. Effect of Different Grazing Intensity on Plant Diversity and Soil Microbial Index

We performed a significant difference analysis for plant diversity and soil microbial diversity (

Table 2. Differences in plant diversity and soil microbial indices under different grazing treatments; different letters indicate significant differences between grazing treatments (p < 0.05).

Treatment	LG	HG
Margalef index	0.87 ± 0.33 a	0.5 ± 0.16 b
Simpson index	0.49 ± 0.18 a	0.36 ± 0.11 b
Shannon index	0.97 ± 0.40 a	0.63 ± 0.17 b
Pielou index	0.67 ± 0.15 a	0.49 ± 0.12 b
Bacterial Chao1 index	2375.67 ± 162.80 a	2256.71 ± 98.99 b
Bacterial Ace index	2391.24 ± 164.43 a	2261.94 ± 98.72 b
Bacterial Shannon index	7.67 ± 1.62 a	7.73 ± 1.59 a
Bacterial Simpson index	0.5 ± 0.51 a	0.5 ± 0.52 a
Fungal Chao1 index	814 ± 140.32 a	657.98 ± 186.70 b
Fungal Ace index	801.76 ± 126.17 a	656.14 ± 174.98 b
Fungal Shannon index	5.55 ± 0.95 a	4.86 ± 0.61 b
Fungal Simpson index	0.64 ± 0.45 a	0.32 ± 0.43 a

). All plant diversity indices showed significantly higher values (p < 0.05) in the light grazing treatment than in the heavy grazing treatment. The Chao1 index and ACE index of soil microorganisms under both grazing treatments showed significantly higher values (p < 0.05) for light grazing than heavy grazing. The Simpson index of soil microorganisms was not significantly different under the two grazing treatments. The Shannon index of soil fungi also showed significantly higher values (p < 0.05) for light grazing than heavy grazing, but the Shannon index of soil bacteria did not differ significantly under the two grazing treatments.

3.2. Effect of Different Grazing Intensity on Soil Chemical Composition

We analyzed the differences in soil physicochemical properties under different grazing treatments for significance, and the results are shown in

Table 3. Differences in soil chemical composition under different grazing treatments; different letters indicate significant differences between grazing treatments (p < 0.05).

Treatment	LG	HG
Soil organic carbon (g/kg)	11.52 ± 3.09 a	9.69 ± 2.98 b
Total nitrogen (g/kg)	1.48 ± 0.42 a	0.73 ± 0.14 b
Total phosphorus (g/kg)	0.86 ± 0.03 a	0.78 ± 0.02 b
Total potassium (g/kg)	23.76 ± 0.66 a	21.45 ± 1.06 b
Available nitrogen (g/kg)	201.32 ± 67.03 a	79.69 ± 15.38 b
Available Phosphorus (g/kg)	17.39 ± 2.53 a	12.25 ± 1.54 b
Available potassium (g/kg)	653.35 ± 142.64 a	453.11 ± 57.69 b
pH	8.86 ± 0.52 a	8.95 ± 0.33 a

. Except for soil pH, the remaining soil chemical components showed significantly higher values in the light grazing treatment than in the heavy grazing treatment (p < 0.05), and there was no significant difference in soil pH between the two grazing treatments.

3.3. Effects of Different Grazing Intensities on Single Ecosystem Functions and Multifunctionality

We calculated the ecosystem individual function index and EMF index using the single function method and the mean method, respectively (Figure 2a). The grassland ecosystem multifunctional index (EMF) ranged from 0.216724 to 0.569721, with a mean value of 0.381154. All functional indices except GP exhibited higher values for LG than HG. Among them, PD, FD and BD exhibited significantly higher values for LG than HG (p < 0.05), in addition to EMF and SF exhibiting significantly higher values for LG than HG (p < 0.01). In contrast, GP showed significantly lower values for LG than HG (p < 0.05), and the remaining SCS showed no significant difference between grazing methods. Finally, we used the mean values of multifunctional index and single function index under different grazing to plot Figure 2b, and we found that in terms of ecosystem function contribution, the largest contribution to EMF at LG level was BD, and the largest contribution to EMF at the HG level was SCS.

3.4. Relationship between Biodiversity and Ecosystem Function and Multifunctionality under Different Grazing Intensity

We used the single function method to calculate the single function of the ecosystem and the average method to calculate the EMF to study the relationship between biodiversity and ecosystem function. We took soil pH as an abiotic factor and combined this with plant diversity, soil fungal diversity and soil bacterial diversity to explore the relationships among the four factors and ecosystem function in an arid desert grassland. The final models featured different relationships between predictors and the response for the most significant results.

As shown in Figure 3, plant diversity, soil fungal diversity and soil bacterial diversity were all correlated with ecosystem function and multifunctionality. Among them, SCS and SF showed significant positive correlations (p < 0.05) with plant, soil fungal and soil bacterial diversity in terms of single ecosystem functions (Figure 3a,b). GP showed significant negative correlations (p < 0.05) with plant diversity (Figure 3(c1)), whereas no significant correlations were found with soil fungal diversity and soil bacterial diversity. EMF showed significant positive correlations (p < 0.05) with plant, soil fungal and soil bacterial diversity (Figure 3(d1–d3)). However, no significant relationship was found between soil pH and ecosystem function and multifunctionality. Overall, the correlation between ecosystem function and biodiversity was greater than that with abiotic factors, indicating the important role of biodiversity in maintaining EMF.

The results for EMF calculated using the multi-threshold method are shown in Figure 4. The effect of biodiversity in relation to ecosystem multifunctionality is essentially positive, and the effect curve is, approximately, a single-peaked curve which increases and then decreases. The maximum effect of plant diversity on multifunctionality (Tmde) is 51% (Figure 4a), and the maximum effect of plant diversity (Rmde) is 0.42 (meaning that for each additional species, the number of functions will increase by 0.42). The threshold value (Tmde) for the maximum effect of fungal diversity on multifunctionality was 47%, and the maximum effect of fungal diversity (Rmde) was 0.52 (Figure 4b). The threshold value (Tmde) for the maximum effect of bacterial diversity on multifunctionality was 42%, when bacterial diversity had a maximum effect (Rmde) of 0.3806 (Figure 4c). The threshold value for the maximum effect of soil pH on multifunctionality (Tmde) was 32%, and the maximum effect of soil pH (Rmde) was 0.79, meaning that for each additional species, the number of functions would increase by 0.79 (Figure 4d). It is worth noting that compared with biodiversity, when the threshold value exceeds 90%, the effect value of pH in relation to EMF is less than zero, and the number of ecosystem functions shows a decreasing trend.

3.5. Effect of Biotic and Abiotic Factors on Ecosystem Multifunctionality under Different Grazing Intensity

Finally, to gain insight into the mechanism of biodiversity regulation of EMF, we modeled structural equations under grazing disturbance conditions with soil pH as an abiotic factor and combined this with biotic factors such as plant diversity, soil fungal diversity and soil bacterial diversity. Overall, based on the chi-square test p > 0.05, RMSEA = 0.00, NFI = 0.999, and the rest of the fit indicators met the model fit criteria, indicating a good fit for this model. Structural equation modeling indicated that grazing could directly affect EMF and explain 22.6% of the variation in EMF. In addition, this explains that grazing can indirectly affect EMF through plant diversity, soil fungal diversity and soil bacterial diversity. However, the abiotic factor, soil pH, had no direct effect on EMF, and its indirect effect on EMF was mediated through plant diversity (β = −0.55, standardized coefficient) and soil bacterial diversity (β = −0.46, standardized coefficient). Among different aspects of biodiversity, plant diversity had the greatest effect on EMF (p < 0.001, β = 0.52, standardized coefficient), followed by fungal diversity, and bacterial diversity had the least effect.

4. Discussion

4.1. Effect of Biotic and Abiotic Factors on Ecosystem Multifunctionality under Different Grazing Intensities

The impact of grazing on grassland ecosystems is mainly through livestock feeding, trampling and excretion, which further affects ecosystem functions, including causing a decrease in nutrient cycling [46], productivity [47] and soil carbon storage [48]. In this study, both ecosystem single function and multifunctionality were calculated using the single function method and the averaging method. Except for the GP function, all the ecosystem single functions showed higher values of LG than HG. PD can reflect the degree of disturbance to the community through competition between species [49], and its decrease is mainly due to the grazing process of livestock foraging and trampling on the plant community structure, resulting in the reduction in dominant species in the plant community [50]. At the same time, the excretions of livestock will increase with the change in grazing disturbance, which in turn will indirectly affect the plant community. SFD and SBD can be important indicators for evaluating the quality of the soil environment [51,52] and play a crucial role in the cyclic transformation of soil nutrients. The grazing process of livestock causes changes in the soil microenvironment [53,54], leading to an impact on the structure of the soil microbial community. The decrease in soil fertility function is due, on the one hand, mainly to the fact that grazing reduces the apoplastic material of plants, which leads to a decrease in the nutrients returned to the soil. On the other hand, the trampling effect of livestock changes the soil’s properties [55,56] and reduces the decomposition of soil microorganisms, thus causing a decrease in soil fertility. The EMF index in this study was averaged from the combination of ecosystem single functions, and the decrease in ecosystem single functions under heavy grazing naturally causes a decrease in the EMF index as well. In addition, in our paper, the GP function showed a significantly lower value for LG than HG. On the one hand, this may be due to the fact that grazing stimulates plant regeneration, and fewer edible forage species are present under heavy grazing; in addition, toxic grasses that are not eaten by livestock, such as Salsola collina, which is large and can grow rapidly, can become the dominant species in grassland [57]. On the other hand, the study area is a desert grassland, and the grassland ecosystem is not well stabilized, with fewer plant species and lower productivity compared with other grassland types. In addition, regarding the contribution of each ecosystem function to EMF, the highest contribution to EMF was made by SBD and SCS under light and heavy grazing, respectively, indicating that the contribution of each ecological function to the ecosystem multifunctionality index was different under different grazing intensities; with the increase in grazing intensity, the contribution to ecosystem function was affected to different degrees, and then the contribution of single ecosystem functions changed. As the grazing intensity increased, it affected the ecosystem functions to different degrees, which in turn showed that the contribution of a single function of the ecosystem changed.

4.2. Relationship between Biodiversity and Ecosystem Multifunctionality under Grazing Disturbance

In our paper, we used various methods (single function method, mean method and multi-threshold method) to calculate EMF, and we explored the relationship between biodiversity (plants, microorganisms) and EMF in desert grassland ecosystems on the northern slopes of the Xinjiang Tian Shan Mountains. Our results were obtained using the single function method to calculate ecosystem single functions and using a general linear model showing that both SF and SCS were significantly and positively correlated with plant diversity, bacterial diversity and fungal diversity under different grazing conditions, similar to the results of previous studies [58,59,60]. SF and SCS represent the nutrient properties of the soil, and plants and soil microorganisms are closely linked to the soil. Plants and soil microorganisms depend on soil fertility and soil carbon accumulation for growth and development, and further carbon sequestration by plants and decomposition by soil microorganisms generate feedback on soil fertility and soil carbon accumulation [61]. In addition, GP was significantly and negatively correlated with plant diversity, a finding that differs from those of most scholars [4,62]. However, some scholars have also shown that GP is negatively or not correlated with plant diversity [63,64,65,66]. In general, plant diversity determines the productivity of grasslands [67,68], and an increase in plant diversity can be directly expressed as an increase in community biomass and apoplankton biomass; high community biomass is conducive to maintaining ecosystem stability [69,70,71] and resisting invasive alien organisms. The reasons for this situation may be, on the one hand, that the trampling, feeding and excretions of livestock can affect the growth and development of plants, as well as cause plant death [72,73], and can also indirectly affect the interspecific competition of plants [74]. On the other hand, regarding the excrement of livestock as the main method of nutrient input to grazing land, the accumulation of excrement will directly change the soil nutrient and moisture conditions [75,76] and play a role in improving the plant growth environment. The interaction of climate, animals, plants and microorganisms under the influence of long-term grazing causes changes in the grassland ecosystem environment [77,78,79]. Because of the selectivity of livestock feeding, grazing reduces plant populations that are more palatable in the plant community but has less effect on plant populations that are not preferred by livestock [57]. In addition, because the dominant species are more widely distributed in grasslands, they are subject to greater trampling by livestock relative to other plant populations [80,81,82]. Furthermore, because of the competitive relationship of plant populations [83,84], a vicious circle can be created for the growth of the dominant species in the plant community. In the long run, grazing will change the spatial distribution pattern of grassland plant populations. The above reasons lead to a significant negative correlation between plant diversity and grassland productivity under grazing disturbance. Based on this, the management of degraded grasslands needs to take into account the threshold of grazing intensity, as well as the presence of livestock not consuming forage and toxic grasses, so that the grasslands can be used and protected more efficiently.

To reasonably assess the importance of biodiversity to ecosystem function, we calculated a significant positive correlation (p < 0.05) between EMF and biodiversity using both the mean method and the multiple threshold method. Our results suggest that maintaining the multifunctionality of desert grassland ecosystems on the northern slopes of the Tian Shan requires higher biodiversity, as shown in previous studies [85,86,87]. On the one hand, a single aspect of biodiversity cannot drive all ecosystem functions, and different species may drive different ecosystem functions [88,89]. In the Tien Shan desert grassland ecosystem, communities composed of plants, soil bacteria and soil fungi interact with each other and can drive different ecosystem functions simultaneously. On the other hand, species communities live together in an ecosystem with relatively limited resources, and there may be ecological niche differences between communities, when ecosystems with more species can make better use of resources [68] and therefore can increase EMF. In the application of the multi-threshold method, Byrnes [44] pointed out that biodiversity can effectively drive multifunctionality when Tmin is low, Tmax is high and both Tmde and Rmde are large. In our study, Tmin and Tmax for plant diversity were 9% and 83%, and Tmde and Rmde were 51% and 0.4158, respectively. Tmin and Tmax for fungal diversity were 10% and 60%, and Tmde and Rmde were 47% and 0.5220, respectively. Tmin and Tmax for bacterial diversity were 21% and 97%, and Tmde and Rmde were 42% and 0.3806, respectively. The Tmin and Tmax for bacterial diversity were 21% and 97%, and Tmde and Rmde were 42% and 0.3806, respectively. In addition, the possible percentages of maximum diversity effect achieved for all three types of species diversity, Pmde, were above 65%, indicating that plant diversity, soil fungal diversity and soil bacterial diversity can drive ecosystem multifunctionality well. This is not quite consistent with the results of van der Plas [90] in European forests. The main reason may be that desert grassland ecosystems, the subject of this paper, are less stable than forest ecosystems and are vulnerable to environmental factors, and this result also leads to the relative importance of species diversity in desert grasslands to EMF compared with European forests. In addition, due to the coordination of trade-offs among ecosystem functions, a species supporting a high level of some functions will inevitably reduce the level of others [36], a result that is also reflected in our paper. At a high threshold level, high biodiversity does not guarantee that all ecosystem functions perform at the highest level at the same time, a result that can also be verified in the studies by Zavaleta [40] and Byrnes [44]. Based on such results, in the subsequent restoration and management of degraded grassland ecosystems, attention should be paid to the conservation of biodiversity so that biodiversity and EMF are at a high level of coordination. The impact of biodiversity loss on the overall ecosystem function should be reasonably avoided to ensure the simultaneous performance of multiple ecosystem functions and to ensure the sustainable use of the ecosystem to a greater extent.

The results of our structural equation modeling suggest that grazing can affect EMF directly and also indirectly through plant diversity, soil fungal diversity and soil bacterial diversity. In addition, we also found that combining above- and below-ground biodiversity can better predict EMF. This result can also be corroborated with the above-mentioned results suggesting that biodiversity can predict EMF and is also consistent with the results of previous scholars [91,92,93]. On the one hand, as grazing intensity increases, grazing can have an impact on plant, soil fungal and soil bacterial community characteristics [94,95,96], resulting in impacts on both plant and soil microbe-driven functions. For example, grassland productivity will decrease due to grazing by livestock [97], and soil nutrients and nutrient cycling driven by soil microorganisms will also both be reduced by grazing [94,98]. On the other hand, plants and soil microorganisms interact with each other under grazing disturbance; the root systems of plant communities provide large amounts of nutritious organic material for the growth and development of soil microorganisms [99], and microorganisms provide nutrients to plant communities for growth and development by consuming organic material [100]. In addition, there is a plant–soil feedback (PSF) relationship between plant communities and soil microbial communities, in which soil microbes can strongly influence plant interspecific relationships, community dynamics and ecosystem function [101,102,103], resulting in changes in the characteristics of plant populations. Our findings provide strong evidence that above- and below-ground biodiversity under grazing disturbance can effectively communicate above- and below-ground ecological functions, and combining evaluations of above- and below-ground biodiversity can better sustain changes in EMF.

4.3. Weaknesses of the Study

We only used soil pH as a representative abiotic factor to explore the effect on EMF, but the direct effect of soil pH on EMF was not significant in the study (Figure 5) and only influenced EMF through plant diversity and fungal diversity. Soil pH was not significantly correlated with either ecosystem function or multifunctionality. The relationship between EMF and soil pH calculated by the multi-threshold method was not significantly representative, and soil pH did not appear to drive EMF. Therefore, more environmental factors should be used to represent abiotic factors in subsequent experimental studies in order for the drivers of EMF to be better explored. In addition, no blank control was set for the measurement of grazing disturbance in our study, which needs to be combined with ungrazed sample plots in the follow-up experiments so that the management system of degraded grassland restoration can be better explored.

5. Conclusions

Our study provides new insights into the relationship between EMF and biodiversity in desert grasslands under grazing disturbance, highlighting the relative importance of plant diversity and soil microbes in maintaining EMF. We found that grazing can reduce EMF directly and indirectly through above- and below-ground biodiversity, and that a combination of high above- and below-ground biodiversity can better maintain grassland ecosystem function and multifunctionality on the northern slopes of Tian Shan. Therefore, in the subsequent evaluation of grassland ecosystem function and multifunctionality, we need to consider both above-ground and below-ground biodiversity, and we should adjust grazing duration and optimize grazing intensity. However, in this study, only soil bacterial diversity and soil fungal diversity were used, which could not fully represent belowground biodiversity, and biodiversity has multidimensional and multiscale attributes; more biodiversity attributes should be considered in the follow-up study in combination with multifunctionality. In addition, more environmental factors should be used to characterize abiotic factors so as to better predict EMF and to also better explore the effects of grazing on grassland EMF.