Study area
The study was conducted in the middle part of the Northern Slope of Tianshan Mountain, Changji City, Ashli Kazakh Ethnic Township (43°59´~ 43°60´N, 86°67´~ 86°71´E), at an altitude of 2200–2400 m (See Fig. 1 for location of the study site) (Li et al. 2021). The region has a dry, typical continental climate with more precipitation but only in summer, an annual rainfall of 450–500 mm, annual temperature of 6.5°C, and frost-free period of 100 d. The soil type is mountain chernozem soil, and the grassland type is mountain meadow-grassland. The common species mainly include Stipa capillata, Carex stenocarpa, Festuca Ovina, Alchemilla Tianschanica, Potentilla Bifurca, and Heteropappus altaicus. The grazing mode in the study area is mainly free grazing, and the main grazing animals are mixed sheep, cattle, and horses, supplemented by cattle and horses. The sheep breeds are mainly Xinjiang fine wool sheep or Kazakh sheep suitable for local nomads.
Experimental design
In this paper, six grazing sites were selected as replicates that were as similar as possible in terms of grassland topography, soil conditions, grassland area, and livestock species and numbers. Then, the grazing intensity was specifically classified according to the dominant plant community species in the grazed grasslands, and each grazing site was divided into three grazing experimental groups: lightly grazed grassland, heavily grazed grassland, and control group (ungrazed grassland) according to previous studies (Li et al. 1997; Yi et al. 2012; Wen et al. 2013). Among them, the dominant species of the grazed grassland plant community was Achnatherum Inebrians, and this grazing site was designated as a heavily grazed (HG) area. The dominant species of the main plant community of the grazing site were Carex stenocarpa and Stipa capillataw, designating this grazing site as a light grazing (LG) area. Finally, the fenced grassland (ungrazed grassland) will be used as an anthropogenic control sample (CK) in this paper in 2018; the dominant species of plant communities in this grassland were Stipa capillata, Festuca ovina, and Poa pratensis. The type of grassland in the three grazing experimental groups mentioned above was mountain meadow grassland, and the differences in the dominant species of plant communities in each experimental group were all caused by the differences in grazing intensity (Fig. 2).
Sampling and measurement
Parallel sampling strips were set up as grassland sampling points in each of the three grazing experimental groups in the above experimental area, seven 1×1 m sample squares (spaced about 50 m) were set up in each parallel sampling strip, and then followed by vegetation surveys in August 2020 and 2021. The sample plots were selected to avoid areas with livestock manure as much as possible to avoid experimental errors. After the sample plots were set up, the plant species composition of the grassland sample plots was recorded, and the natural height of the grassland plants, grassland cover (visual method) and the number of each species were measured. Finally, all plants in the grassland sample plots were cut flush to the ground, weighed fresh and placed in named envelope bags, then brought back to the laboratory to be dried at 65°C in order to determine the dry weight of above-ground biomass.
Meanwhile, based on the grassland vegetation survey, soil from 0–10 cm soil layer was collected by soil auger method (three replicates per sample square) in this paper. Then, the soil samples from seven sample squares were evenly mixed into one soil sample, and the soil was removed from the adulterated plant roots and gravels with a 2-mm sieve. Finally, the total soil sample was divided into two parts, one of which was placed in sterile centrifuge tubes, which in turn were then placed into an ice box. These boxes were immediately transported back to the laboratory and stored at -20°C to characterize the soil microbial community. The other soil samples were placed in sealed bags and labeled, and returned to the laboratory to assessing the physicochemical properties of the soil after natural air drying.
Sample collection and index determination
Determination of soil physical properties
According to the soil agro-chemical analysis undertaken by Bao (2000), the physical and chemical properties of soil such as organic carbon, soil total nitrogen, soil total phosphorus, soil total potassium, soil available nitrogen, soil available phosphorus, soil available potassium, and soil pH were typically measured (Table 1).
Table 1
Methods of determination of soil physical and chemical indexes
Indicators
|
Assay method
|
Organic carbon
|
Potassium dichromate-concentrated sulfuric acid external heating method
|
Total nitrogen
|
Kjeldahl method
|
Total phosphorus
|
Sodium hydroxide alkali fusion-molybdenum antimony anti-colorimetry
|
Total potassium
|
Sodium hydroxide alkali melt-flame spectrophotometry
|
Available nitrogen
|
Sodium hydroxide alkaline diffusion method
|
Available phosphorus
|
Sodium bicarbonate leaching-molybdenum antimony anti-colorimetric method
|
Available potassium
|
Ammonium acetate extraction - flame photometry
|
pH value
|
Potentiometric method (water-soil ratio 5:1)
|
Extraction and sequencing of soil microbial DNA
In this study, total DNA was extracted from soil microorganisms using the MP-FAST DNATMSpin Kit for Soil (Omega Bio-Tek, Norcross, GA, USA). The extracted genomic DNA was detected with 1% agarose gel electrophoresis. The corresponding regions of the bacterial 16S rRNA gene and fungal ITS rRNA gene were amplified with PCR using respective primers. The reaction conditions were as follows: pre-denaturation at 95°C for 3 min, 27 cycles (denaturation at 95°C for 30 s, annealing at 55°C for 30 s, extension at 72°C for 45 s), and a final extension at 72°C for 10 min. The PCR products were pooled and detected using 2% agarose gel electrophoresis. The PCR products were quantified using a QuantifluorTM-ST Blue Fluorescence Quantification System (Promega Company). The PE250 library was constructed from the purified amplified fragments according to an Illumina MiSeq platform (Illumina, San Diego, USA). Sequencing was undertaken with an Illumina Misqe PE250 platform. Shanghai Meiji Biomedical Technology Co., Ltd completed the library construction and sequencing tasks; Mothur software was used to calculate the α diversity index of the sample.
Quantification of ecosystem multifunctionality
Based on a previous system of evaluating EMF (Jing et al. 2015; Xiong et al. 2016; Wang et al. 2019), five types of functions (13 variable indicators) related to grassland ecosystem were selected in this paper, taking into account the regulatory, supply and support service functions of grassland ecosystems. These included water regulation (MR), soil fertility (SF), nutrient transformation and cycling (NC), soil carbon storage (SCS), and grassland productivity (GP). Finally, this study also established an evaluation system of mountain meadow-grassland EMF on the north slope of Tianshan Mountain (Table 2).
Table 2
Evaluation system of grassland ecosystem multifunctionality
Type
|
Ecosystem functioning
|
Ecosystem function indicators
|
Tuning services
|
Moisture regulation (MR)
|
Soil bulk density; Soil moisture content
|
Supply services
|
Grassland productivity (GP)
|
Aboveground biomass; Underground biomass
|
Support services
|
Nutrient conversion and cycling (NC)
|
Soil carbon-nitrogen ratio; Soil pH
|
Soil fertility (SF)
|
Soil total nitrogen; Soil available nitrogen; Soil total potassium; Soil available potassium; Soil total phosphorus; Soil available phosphorus
|
Soil carbon storage (SCS)
|
Soil organic carbon
|
Soil nutrients (SS)
|
Soil conductivity
|
Statistical analysis and data processing
Calculation of diversity index of the plant community
Shannon-Wiener Diversity Index: H′ = -∑PilnPi(Pi = Ni/N) (1)
In the formula, Pi 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.
Soil microbial data processing
In this paper, the data were spliced using FLASH software and quality controlled using Trimmomatic software to obtain valid sequences. The obtained sequences were clustered according to a 97% similarity criterion. Then the corresponding diversity index of α could be obtained by Mothur software. Comparisons were made with the Silva (Release128) database, species classification of different sequences were made, and the sequences annotated. Finally, the community composition of each sample was counted, and the species diversity in each sample was analyzed.
Calculation of the Ecosystem Multifunctionality Index
In this study, the averaging method was used to evaluate the electromagnetic fields as follows:
All indicators (Table 2) were standardized to the same order of magnitude using the min-max standardization method. The maximum value of an indicator is the average of the top 5% observed value of each indicator, and the minimum value of an indicator is the bottom 5% observed value of each indicator:
Min-max normalization calculation:\({f}_{ij}=\left({x}_{ij}-{min}_{ij}\right)/\left({max}_{ij}-{min}_{ij}\right)\) (2)
f ij is the normalized value of the type j ecosystem function variable of sample i, xij is the measurement value of the type j ecosystem function variable of sample i, minij is the minimum value of the j ecosystem function variable across all plots of the same factor, and maxij is the maximum value of the j ecosystem function variable across all plots of the same factor.
Ecosystem single-function index: Calculates the single-function index that reflects the different functions of the ecosystem (Table 2 for details).
Ecosystem Single Function Index F:\({F}_{ij}=\sum _{j}^{n}{f}_{ij}/n\) (3)
F ij is the functional index of the j function of the sample i, and n is the number of ecosystem variable indicators contained in the function.
index of EMF:\({EMF}_{i}=\sum _{1}^{N}{f}_{ij}/N\) (4)
EMF i is the ecosystem multifunctionality index of sample i, calculated based on the standardized mean of all variable indicators in the sample, N is the number of ecosystem functions contained in sample i.
Analysis of trade-offs and coordination among ecosystem functions
This study used Pearson correlation analysis to study trade-offs and coordination between different functions. If there was a significant positive correlation between the two functions (r > 0, p < 0.05), it indicated that the two functions are synergistic. If there was a significant negative correlation between the two functions (r < 0, p < 0.05), it indicated that the two functions were a trade-off relationship. If they were not correlated, no association was present.
The data were classified and analyzed by using Excel 2019 software. To analyze the differences in the multifunctional index of lightly grazed, heavily grazed and ungrazed grassland ecosystems, this study used the Duncan's test in SPSS 26.0 software to test for significant differences between the single functional and EMF index. Analysis of Pearson correlations between ecological functions were undertaken in R 4.1.1 software. The relationship between the diversity of plant and soil microbial and ecosystem function was analyzed using Pearson correlation analysis. Finally, this study used structural equation modeling (SEM) to analyze the mediating effect of biodiversity on EMF under grazing conditions. SEM fitting tests were performed on plant diversity, bacterial diversity, and fungal diversity under different grazing intensities to determine the direct and indirect effects of grazing intensity on EMF.