Data analyses
In the first step of our analysis, we investigated the relationships between the intensive cover of agriculture, pasture, urbanization, afforestation, environmental (sediment heterogeneity and water quality deterioration), stream morphology (depth), biogeographical (latitude of stream sites) and climatic variables (MAT and MAP) and (i) taxonomic richness, (ii) functional diversity, and (iii) diversity of three trait categories (recruitment and life-history, resource and habitat use, and body size). To do so, we employed linear mixed-effects models (LMM) in the R ‘nlme’ package (Pinheiro et al., 2013). The two seasonal sampling periods were nested within each site that was nested within each study area (Amazonia and Uruguay) as a random effect. In the models of functional diversity, we included taxonomic richness as a predictor to account for confounding effects of the local species pool (Mayfield et al., 2010; see the Table S2). Using a stepwise regression procedure all the models were reduced to select the best model (i.e., the model with better predictors for each biodiversity attribute of each assemblage) with the lowest Akaike’s Information Criterion corrected for small sample size (AICc). With the best selected models, we then performed a model-averaging procedure based on AICc selection (∆AICc < 2) to determine parameter coefficients for the best final subset of predictors for each response variable (Tables S2 and S3). This procedure was performed using the dredge function in the MuMIn package (Bartoń, 2014). Visual inspection of residuals using graphical diagnostics (QQ plots and residual plots) revealed that the assumptions of normality and homoscedasticity were met. We assessed the multicollinearity of each predictor variable within models by calculating the variance inflation factor (VIF), and we removed all variables with VIF > 3 (here, latitude). All other variables had VIF < 3, indicating they were not highly correlated (Figure S4 and Table S3). A priori , we standardized all predictors (z-scored: centered to mean and divided by the SD) to interpret slope estimates on a comparable scale (Schielzeth, 2010). To infer the relative importance of each predictor on different biodiversity components, we computed the percentage of variance each explained according to Le Provost et al. (2020). To do so, we compared the absolute value of their standardized coefficients with the sum of all standardized regression coefficients of all predictors considered in the best models. This approach is analogous to a variance partitioning analysis because all predictors were previously transformed to z-scores (Le Bagousse-Pinguet et al., 2019).
We also applied structural equation models (SEM) to disentangle the direct and biodiversity-mediated indirect pathways by which intensive land-use cover (ILUC) influenced standing biomass. We fitted separate SEMs to taxonomic richness and functional diversity. Environmental (water quality, heterogeneity), stream morphology (depth), and climate (temperature and precipitation) covariates were also included in SEMs, and all paths specified were theoretically supported (Figure S1). For simplicity, we only use the ILUC, as it is strongly related to agriculture, pasture, and urbanization (Figure S2). To evaluate whether the ILUC effects differed between Amazonia and Uruguay, we constructed SEMs using multi-level analysis considering the two areas separately (results presented in supplementary material, Figure S5, Table S10). This approach allowed us to implement a model-wide interaction and to test each path interaction within the model. To reduce the number of predictors, we performed a model selection based on Akaike’s Information Criteria corrected for small sample size (AICc). In this selection, the full model (including all predictors) was compared with the reduced model using AICc (AICfullmodel – AICreducedmodel). We considered ΔAICc > 2 units as distinguishing between the full and the reduced models. Notably, the full and reduced final models differed in ΔAICc > 10 (Table S8). The lack of direct or indirect effect on standing biomass was used as a criterion to remove predictors. SEMs were constructed with the same random factor design as previous LMMs using linear mixed effects models (Pinheiro et al., 2018) with the ‘psem function’ from the ‘piecewiseSEM’ package (Lefcheck, 2016). We present the standardized coefficient for each path and estimated the indirect effects by coefficient multiplication. Path significance was obtained by maximum likelihood, and model fit was evaluated using Shipley’s test of d-separation through Fisher’s C-statistic (p > 0.05 indicates no missing path). All analyses were conducted in R version 4.1.1. (R Core Team, 2021).