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