Identification of thresholds in the response of diversity and
function to land use
Before testing for thresholds, we checked whether any variables showed
non-linear responses to land use intensity and its three components. To
detect non-linear responses, we compared three models: linear
regressions, quadratic regressions, and generalized additive models
(GAM) (Hastie & Tibshirani 1990). Quadratic models allow for non-linear
effects, with one turning point in the relationship between the
variables, while GAMs allow for complex non-linear responses, with
several turning points. Both quadratic models and GAMs are often used to
identify ecological thresholds (Berdugo et al. 2020). The best
fitting model for each response was selected using the Akaike
Information Criterion. Specifically, we considered that quadratic models
and GAMs had better fits than the linear model when they had an AIC
values at least 2 units lower than the linear model (ΔAIC =
AIClin – AICqua,GAM > 2)
(Burnham & Anderson 2003). Then, for variables that showed a nonlinear
response, we tested whether they also had a threshold (i.e. an
abrupt shift) in their response to LUI, and its three components, using
segmented regression (Muggeo 2003). Segmented regressions were only
applied to those variables that showed non-linear responses because they
force break points in the explanatory variable, potentially leading to
overfitting and spurious thresholds, if the response is linear. In
addition, if GAMs suggested the existence of several turning points, we
fitted segmented regressions with more than one break point. When the
segmented regression model was significant and indicated a better fit
than the linear model (based on ΔAIC > 2), we considered
that the variable showed a threshold response. To control for the effect
of region on diversities and functions, we used the residuals of each
variable from a model including region as the explanatory variable. All
variables were standardized by subtracting the mean and dividing by the
standard deviation before running the models. We used the gam(Hastie 2019) and segmented (Muggeo 2008) packages in R.
As land use intensity varied over time in each plot we always tested the
response of each diversity and function to the land use intensity the
previous year (e.g. the response of plant diversity in 2013 to
grazing in 2012). Similarly, when variables were measured in several
years, the mean value for the time period was used for both explanatory
and response variables (e.g. the response of average plant
diversity during 2011-2013 to average grazing during 2010-2012). We
decided to use the land use value from the previous year because the
study area is continuously managed and the different taxa and ecosystem
functions may be measured before much of the land use in a given year
has occurred. However, although LUI changes over time, the changes are
not sufficient to dramatically alter the land use levels of the plots:
both LUI and its components showed strong positive correlations between
years (mean r ± standard error; LUI r = 0.7 ± 0.01; grazing r = 0.65 ±
0.01; mowing r = 0.74 ± 0.01; fertilization r = 0.71 ± 0.01) and all
years were strongly related (r > 0.74) to the mean values
over the total time period (Supplementary Figure S4).