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