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