Materials and Methods
Data preparation. To investigate the integrative effect of grazing intensity, duration, livestock type and climatic factor on multiple biodiversity and ecosystem functioning, we searched peer-reviewed publications during 1900–2019 using Google Scholar, Web of Science, and China Knowledge Resource Integrated Database with the search terms: (grazing) AND (diversity or richness) AND (ecosystem functioning) OR (productivity/stability/soil properties /soil nutrients/carbon exchange).
We examined the identified publications according to the following criteria: (a) field-manipulated experiments included a control (no grazing treatment), which had similar conditions to the grazing plots, such as microclimatic factors, vegetation and soil types. (b) Grazing intensity was clearly quantified in terms of vegetation utilization rate or stocking rate. (c) Experimental duration, livestock type and climatic characteristics (i.e., mean annual temperature and precipitation) were clearly indicated. (d) The means, standard deviation, and the number of replicates were explicitly provided. We obtained the table-form data directly and extracted graphical data using Getdata software (GetData Pty Ltd, Kogarah NSW, Australia) from the original publications. In total, 104 publications that investigated the effects of grazing intensity on biodiversity and individual ecosystem functions were used to establish a global dataset containing 494 independent observations (Fig. S1).
In this dataset, mean annual precipitation (MAP) ranged from 76 to 1834 mm and mean annual temperature (MAT) ranged from −3°C to 26.6°C. Aridity index was calculated as mean annual evaporation (provided by the original literatures or local weather records of the experimental site) /MAP to quantify dry conditions. Livestock type was categorized as small herbivores (i.e., sheep and goats in 62 publications), and large herbivores in 35 publications, including cattle, yak, horse (1 publication), bison (1 publication) and wildebeest (1 publication), while both cattle and sheep were used in 7 publications. Grazing duration was represented by experimental duration and ranged between 1 and 75 years, which was calculated from the establishment of the control (grazing exclusion treatment), because there is usually an untraceable long-term grazing history on natural grasslands. In addition, three grazing intensities were classified (light, moderate and heavy grazing) according to the description in the original publications. We also calculated the percentage change in aboveground biomass (AGB) from the 73 publications containing AGB in the dataset to verify the level of grazing intensity (Tang et al.2019). Our results showed that the percentage decline in AGB gradually increased with rising grazing intensity, indicating that the grazing intensity estimates in the original publications were reliable (Fig. S2).
Our dataset collected 5 biodiversity components (i.e., plant species diversity, diversity of plant functional group, plant functional diversity, insect diversity and soil microbial diversity). Among them, plant species diversity included 5 metrics: species richness, Shannon-Wiener, Margalef’s, Pielou’s and Simpson index. Diversity of plant functional groups was indicated by Shannon-Wiener, Margalef’s and Pielou’s index. Plant functional diversity was indicated by Rao’s quadratic entropy’s, evenness and dispersion index. Insect diversity (including herbivores and predators) was indicated by richness and Shannon-Wiener index. Soil microbial diversity (represented by bacteria diversity) was indicated by richness, Shannon-Wiener and Pielou’s index. In addition, 12 individual ecosystem functions were collected in our dataset, including aboveground biomass (AGB), belowground biomass (BGB), temporal stability of vegetation productivity (TS), soil organic carbon (SOC), soil total nitrogen (STN), soil available nitrogen (SAN), soil total phosphorus (STP), soil available phosphorus (SAP), soil moisture (SM), net ecosystem productivity (NEP), ecosystem respiration (ER) and gross ecosystem productivity (GEP). For each publication, AGB and BGB were measured at the same period of peak vegetation biomass, and soil properties were surveyed at 0-30 cm depth.
Data analysis. We employed the natural log-transformed response ratio (lnRR) to evaluate grazing effects on all biodiversity metrics and individual ecosystem functions, which was calculated as
ln(RR)=ln(Xgrazing / Xcontrol)
where Xgrazing and Xcontrol represent the observed values of the selected variables in grazed and control plots, respectively(Hedges et al.1999; Tang et al. 2019).
We weighted each observation using its sample size as in previous studies (Ma & Chen 2016; Tian et al. 2018; Chen & Chen 2019):
Weighting = Ngrazing × Ncontrol/(Ngrazing + Ncontrol)
Where Ngrazing and Ncontrol indicate the number of replications for the response variables in grazed and control plots, respectively.
We employed a linear mixed effect model to analyze the overall effects of grazing intensity on multiple biodiversity groups and individual ecosystem functions and their 95% confidence intervals (Fig. S3). We examined the links between plant species diversity and other biodiversity groups (Fig. S4), the relationships between aboveground biomass of plant community with other individual functions (Fig S5), as well as the relationship between plant species diversity and individual functions (Fig. S6) under grazing using the linear regression analysis, because most of the studies in our database have species diversity and aboveground biomass of plant community that can match other indicators. Then, we found that these individual biodiversity and ecosystem functions had the consistent negative responses to increasing grazing intensity (Fig.S3), and there were no possible trade-offs among different biodiversity groups, individual ecosystem functions and the BEF relationship (Fig. S4-S6). Moreover, although EMF had more significant correlations with plant diversity (including species diversity, functional groups diversity and functional diversity) than insect and soil microbial diversity (Fig. S9), the slope of the BEF relationship significantly increased with increasing number of biodiversity groups or individual functions, and so is the R value of correlations (Fig. S10). These results provide a key motivation for further exploration of multi-diversity and EMF in response to grazing intensity.
We calculated multi-diversity (including 5 biodiversity groups) and EMF (including 12 ecosystem functions) using the averaging approach(Manning et al. 2018; Wang et al. 2019), with each variable naturally log transformed (lnRR) and weighted according to the number of indicators provided before averaging in each study. Then, we analyzed the effects of grazing intensity and its interactions with grazing duration, livestock type and aridity index on multi-diversity and EMF, with all these factors taken as fixed effects, and each “study” set as a random effect to account for possible autocorrelation among observations in each study (Fig. 2). As a supplementary method, we took the classified indicators of different biodiversity groups and individual ecosystem functions as the random effects, and nested them into “study” using linear mixed effect model to avoid the possible autocorrelation among biodiversity groups or individual functions, and then we detected the similar results as the averaging approach (Fig. S7, S8). For the analysis, all variables were examined for normality and homogeneity. The treatment effects were considered significant at α = 0.05, or if the 95% CIs does not cover zero. Linear regression was used to examine the relationships between lnRRs of individual biodiversity or multi-diversity and EMF under different grazing intensities, and the interaction term “grazing intensity×lnRR (multi-diversity)” specifically tests whether the effects of grazing intensity on EMF are dependent on multi-diversity. These analyses were conducted using the linear mixed effect model of lme4 package (Bates et al. 2015), and the correlations between multiple biodiversity and ecosystem functioning also were checked using model II regression of lmodel2 package.
Structural equation model. We employed a structural equation model (SEM) to evaluate the simultaneous effects of grazing intensity, duration, livestock type, aridity index and their interactions on the lnRRs of EMF directly and indirectly via lnRRs of multi-diversity (conceptual model), and then we selected the final model by excluding the insignificant factors (i.e., aridity index as well as the interactive effects of grazing intensity and other factors) based on goodness-of-fit statistics and lowest AIC value. The SEM was implemented using the “piecewiseSEM 1.2.1” package to account for the random effects of “Study” (Lefcheck, 2016). All statistical analyses were performed in r 3.5.2 (R Core Team, 2018).