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).