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.