Data analysis
Univariate analyses: grass cover and species richness
The relationship between herbivore abundance and species richness, habitat, bedrock, and their mutual two-way interactions (predictors) on grass cover and species richness (response variables) were explored by the linear mixed-effect models – LMMs (Table 1). Data on abundance and species richness of herbivores were merged for both seasons, as separate models for dry, wet and both seasons merged gave similar results (not shown). Triplet identity, which reflected the spatial clustering of the plots, was set as a random variable. For grass cover and species richness, we used LMMs with normal distribution; the response variable grass species richness was square-root transformed to improve normality and homogeneity of variance. Herbivore abundance was log-transformed in LMM to reduce the skewness of data. To test possible non-linear effects of herbivore abundance and herbivore species richness, we added their quadratic terms to the models; however, as they were not significant, they were not included in the final models. The significance of the terms in LMM models was obtained via anova function. The marginal coefficient of determination R 2, which represents the variance explained by the fixed effects, was computed using r.squaredGLMM function from the “MuMIn“ package (Bartoń 2022). Marginal R 2 for individual terms (predictors) was computed by comparing models with and without the particular term, obtained by model backward simplification. Here we present only R 2 for significant terms. Regression lines in figures were plotted based on simplified LMM models that contained only the predictor and response variables shown on the given graph. Regression estimates (marginal means and their confidence intervals) presented in graphs (Fig. 5, 7 and 8) were computed byggpredict function from the “ggeffect” package (Lüdecke 2018). The differences among individual habitats in the simplified models with only the predictor and response variables were further tested by the post-hoc Tukey HSD pairwise comparison of estimated marginal means (Lenth 2021). In addition to the R-base packages, we further used the packages “nlme” (Pinheiro et al. 2021) for fitting the linear mixed-effect models, and package “emmeans” was used for subsequent multiple comparisons among significant terms (Lenth 2021). Graphs were plotted using the “tidyverse” (Wickham et al. 2019), “grid” and “ggpubr” (Kassambara 2020) packages. All computations were done in the programme R 4.2.0 (R Core Team 2022).
Multivariate analyses: grass community composition and mean grass traits
In the first group of analyses, grass species covers were response variables, and (i) the total number of records and the species richness of herbivores , or (ii) the number of records of individual herbivores were predictors (Table 1). To filter out the effect of habitat, bedrock, and spatial arrangement given by the mutual position of permanent plots, we set them as covariates in all analyses. Spatial effects (autocorrelation of plots) were identified by the PCNM analysis (ter Braak and Šmilauer 2012), where GPS coordinates represented the plot position, from which a matrix of spatial PCoA vectors was calculated; here we used the scores of the first three PCoA vectors as spatial covariates. Because plots were grouped in triplets in the field, we used a hierarchical split-plot design, where triplets were set as whole plots, while plots within a triplet represented split plots. Both triplets and plots within triplets were permuted freely, and the significance was tested by Monte-Carlo permutation tests. In the second group of analyses, the response variables were the mean values of traits of the grass community, so-called community-weighted means, calculated as values of species traits weighted by their relative abundance within each plot; the mean values for all traits were computed for each plot (Šmilauer and Lepš 2014). Covariates and the permutation scheme remained the same as in the first group of analyses. All multivariate analyses were performed in Canoco 5 software (Šmilauer and Lepš 2014).