Statistical analysis
Species & phylogenetic diversity
We analyzed the effect of grazing intensity on SR (overall and for the shrubs, herbs, and graminoids), H’, Evar, phylogenetic dispersion and abundance-weighted phylogenetic dispersion using linear mixed effects models with site as a random effect (following Begley-Miller et al. 2014) and grazing, year, a grazing:year interaction, and wet/dry as fixed effects. First, we tested the grazing:year interaction using a likelihood ratio test. If the interaction was not significant (p > 0.05), then we removed it from the model and tested the remaining fixed effects. If a significant effect of grazing or year was found (p < 0.05), then a post-hoc analysis was performed with Tukey pair-wise comparisons. To control for the possibility that changes in phylogenetic dispersion may result from the transition of communities from dwarf shrubs to graminoids, we also tested the effect of grazing on phylogenetic dispersion and abundance-weighted phylogenetic dispersion with the same procedure described above, while controlling for the proportion of species that are graminoids (for phylogenetic dispersion) and the proportion of the relative abundance that are graminoids (for abundance-weighted phylogenetic dispersion). All analyses were performed in R version 3.6.0 (R Core Team 2019).
Community structure
To compare the community structure of plots we used the Bray-Curtis dissimilarity index, calculated using vegdist from the R package “vegan” (Oksanen et al. 2016), as a measurement of the distance between plant communities based on our relative species abundance data. To partition variance within the distance matrix, we used a non-parametric permutational multivariate analysis of variance (PERMANOVA), as implemented in the vegan function adonis . Significance values and pseudo F-statistics were obtained from permutations (n = 1000) restricted within each site due to our nested experimental design. Given that this technique allowed us to perform a multivariate analysis, we include grazing, year, a grazing:year interaction, and wet/dry as covariates. When significant values (p < 0.05) were obtained, we performed a post-hoc analysis with Bonferroni corrections to correct for multiple comparisons in the PERMANOVA.
To visualize and corroborate the results of the PERMANOVA, we used a non-metric multi-dimensional scaling (NMDS) from the functionmetaMDS in vegan. NMDS is an ordination technique that represents highly dimensional data by maximizing the correlation of ranked-distances between the original highly dimensional data and a two dimensional representation (Faith et al. 1987, Minchin 1987). A stress score is calculated as a measure of how accurately the two dimensional ordination represents the distances in the original data; stress scores < 0.2 are generally considered acceptable (Clarke 1993). Communities grouped closely together in the ordination space are interpreted as being more similar than those placed farther away.