Diversity-interaction model
To disentangle plant guilds effects on SOC, we built a Diversity-Interaction Model (DIM) (Kirwan et al. 2009) in two steps. Preliminary normality tests (distribution histogram and normal Q-Q plot, and the Shapiro Wilk W test (Shapiro et al. 1968)) showed that the variable SOC was not normaly distributed (W = 0.948; p-value < 0.001). We log-transformed SOC to prevent a breach of the normality assumption by the residuals of the models (W = 0.99; p-value = 0.18). From now on, it should be understood that the mention of SOC in the context of modelling refers always to log(SOC).
As a first step, to control the wide climatic and topographical variability of the Pyrenees (Garcia-Pausas et al. 2007, 2017) we fitted a linear model (Legendre & Legendre 1998) for SOC considering only environmental drivers other than plant guilds. We included both main effects and pair-wise interactions. We selected several candidate terms by a semi-automatic procedure by a genetic algorithm included in the R package glmulti (Calcagno & Mazancourt 2010). We selected the best model including only main effects of the predictors, using SOC as the response variable. Then we modelled the residuals of that model with the same predictors to look for significant candidate level-two interaction terms. Once we obtained a set of candidate model terms (with both main effects and interactions), those were included into a single model on which we performed backward-forward selection. We used several methods to compare and determine the final model, including the corrected Akaike information criterion (AICc; (Symonds & Moussalli 2011)), the adjusted determination coefficient (Radj2) and model comparison techniques with the “anova()” function in R, using Chi-square tests to test whether the reduction in the residual sum of squares was statistically significant.
At the second step, to complete the DIM, we added to our null model the guild proportion variables and the guild pair-wise interactions. Because all guild proportions sum 100%, we built a model without intercept at this step (Kirwan et al. 2009). We discarded curved relationships between the model terms and SOC by comparing this final DIM model with the corresponding generalized interaction model (Connolly et al.2013) using a likelihood test (Chi-squared = 1.27 x 10-3; p-value = 0.97; Hothorn et al., 2019).
In order to show the effects of plant guilds we calculated the predicted effects of the model for all possible combinations of grasses, forbs and legumes with the emmeans package (Lenth et al. 2019), fixing the rest of the variables at their mean value. Sedges were fixed at 0% because of their low effects on SOC (Table 1), their low mean proportion (2%) and to obtain clearer plots. We built a ternary plot with the ggtern package (Hamilton & Ferry 2018) to show predicted SOC variation across plant guild proportions. We also plotted SOC variation across a legume proportion gradient under: 1) a balanced scenario where grasses and forbs accounted equal proportions; 2) a forb-dominance scenario where forb:grass proportions were 80:20; 3) a grass-dominance scenario were forb:grass proportion was 20:80; and 4) a co-dominance scenario where grasses and forbs accounted equal proportions (50:50), of the remaining proportion unaccounted by legumes.