Statistical analyses
All the analyses were performed in R v.3.6.3 (R core team, 2020). We
used linear mixed-effects models from the R package “nlme” (Pinheiroet al. 2017) to account for the nested structure of the data. We
looked at the treatment effects on each plant diversity metric and
stability using the R syntax: lme (y ~ nut*fen, random=
~1|site). Alpha stability, spatial asynchrony,
and gamma stability were log-transformed to improve normality and
homogeneity of variance. To look at whether the effects of herbivore
exclusion and its interaction with nutrient addition on stability
increase as grazing intensity increases, we rerun the above models but
adding the herbivore index as a covariate in the models.
We built a structural equation model (SEM) using the function “psem”
from the R package piecewiseSEM (Lefcheck 2016) to evaluate the direct
and indirect effects of nutrient addition, herbivore exclusion, and
their interaction on alpha and gamma stability. An initial model (Fig.
S4) was built based on prior knowledge (Grman et al. 2010; Wilcoxet al. 2017; Zhang et al. 2019; Hodapp et al. 2018;
Gilbert et al. 2020). Rationales for each link in the initial SEM
are summarized in Table S2. To fit the SEM, we used the function “lme”
with site as the random effect when the component models tested only the
treatment effects, and we used function “lm” when the component models
tested the relative contribution of both treatments and plant diversity
metrics to stability to take into account the diversity gradient among
sites. Alpha stability, spatial
asynchrony, and gamma stability were log-transformed to improve
normality and homogeneity of variance. Several studies suggest that
abiotic variables such as rainfall, temperature, and soil parameters
impact grassland stability (Zelikova et al. 2014; García-Palacioset al. 2018; Zhang et al. 2018; Gilbert et al.2020). We, therefore, performed another SEM including temporal
variability in water balance and spatial variability in soil chemistry
following Gilbert et al. (2020) to test the robustness of our
results. See online supplementary text for more details.
We further analyzed which component (abundance gradients or balanced
variation) of temporal community dissimilarity was more related to alpha
stability and which component of spatial community dissimilarity was
more related to spatial asynchrony. We fitted “lm” models where alpha
stability was the response variable, temporal community dissimilarity
(or each of its components), and its interaction with treatments as the
explanatory variables. We fitted similar models for spatial asynchrony.
We define that an effect is marginally significant when
0.05<p<0.1, while significant when p≤0.05.