Grazing intensity
Following Borer et al. (2020) and Anderson et al. (2018),
we quantified grazing intensity from vertebrate herbivores at each site
using a herbivore index. Specifically, herbivore species
(>2 kg) that consume grassland biomass were documented at
each site by site PIs, and each species was assigned an importance value
from 1 (present, but low impact and frequency) to 5 (high impact and
frequency). An index value was calculated for each site as the sum of
herbivore importance values for all herbivores.
Plant diversity metrics and stability across scales
Following Hautier et al. (2020), we treated each 1 m2subplot as a “community” and the replicated subplots within a
treatment across blocks within a site as the “larger scale” sensu
Whittaker (1972) (see an illustration in Fig. S1). Plant diversity
metrics used in this study included alpha diversity, beta diversity,
Pielou’s evenness, and community dissimilarity metrics. Alpha diversity
is the average number of species (i.e. species richness) recorded in the
three subplots in each treatment at each site. Beta diversity is
calculated as gamma diversity/alpha diversity (i.e. multiplicative beta
diversity), where gamma diversity is the total number of species
recorded in three subplots in each treatment at each site. Pielou’s
evenness was calculated as H/ln (S), where H is Shannon’s diversity
index, and S is alpha diversity.
We calculated community dissimilarity (temporal and spatial community
dissimilarity) using Bray–Curtis dissimilarity metrics based on cover
data. Note that some researchers also refer to temporal and spatial
community dissimilarity as temporal and spatial beta diversity (e.g.
Chalcraft et al. 2008; Dornelas 2014). Community dissimilarity
can arise from two concurrent processes, namely abundance gradients and
balanced variation in abundance (Baselga 2017). Abundance gradients
arise from a simultaneous increase or decrease in the cover of all
species, reflecting changes in the total cover. Balanced variation
arises from replacement among species (i.e. decreases in the cover of
some species are compensated for by increases in other species),
reflecting changes in the relative cover. Temporal and spatial community
dissimilarity may impact gamma stability via alpha and spatial
asynchrony respectively, and their impact may depend on the driving
processes (see Fig. S2 for more details). Therefore, we also look at
which process is driving temporal and spatial community dissimilarity,
and their impact on alpha and gamma stability. Temporal community
dissimilarity of each treatment was calculated as dissimilarity of a
community through the 5-year experiments and averaged over 3 blocks.
Similarly, spatial community dissimilarity of each treatment was
calculated as dissimilarity of 3 blocks in each treatment each year and
averaged over the experimental years. Temporal/spatial community
dissimilarity and the partitioning of it into abundance gradients and
balanced variation were done using the function “beta.multi.abund”
from the R package betapart with the index.family of “Bray” (Baselga
& Orme 2012).
Stability at a given spatial scale was calculated as temporal
invariability: \(\frac{\mathbf{\mu}}{\mathbf{\sigma}}\) , where μ and σ
are the mean and standard deviation of aboveground biomass over the
experimental years. Alpha stability was the stability of aboveground
biomass averaged over three subplots in each treatment at each site;
gamma stability was the stability of total aboveground biomass in three
subplots in each treatment at each site (Wang et al. 2019; Hautier et
al. 2020). Biomass was not detrended because no clear trends were shown
over the 5-year experiment in most sites (Fig. S3). Also, previous
studies using NutNet data show that treatment effects (i.e. nutrient
addition) on stability are quantitatively the same with or without
detrending (Hautier et al. 2020). Spatial asynchrony was
calculated as 1/φ, φ =\(\frac{\sum_{i,j}w_{\text{ij}}}{\left(\sum_{i}\sqrt{w_{\text{ii}}}\right)^{2}}\), where wij is the temporal covariance of aboveground
biomass between local communities i and j, and wii is
the temporal variance of aboveground biomass of local community i. These
variables were calculated using the function “var.partition” (Wanget al. 2019).