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