Statistical analyses
Prior to running models and analyses, we used variance inflation factor
(VIF) analysis to assess multicollinearity among predictor variables,
and sequentially removed those with the highest VIF until all had a VIF
< 3 (Zuur et al. 2010). Total number of flowing day, distance
to the nearest perennial site, and the length of upstream perennial
reaches were collinear with other variables and were thus removed from
analyses. We also used Pearson’s correlation coefficients to assess
correlations between FP and remaining connectivity variables (i.e. DS
and %UC). To characterize site environmental conditions we performed a
principal component analysis (PCA), including all environmental
variables (except FP, DS and %UC) averaged across season for each site.
We further used a permutational analyses of variance (PERMANOVA based on
Euclidean distance) to assess if environmental characteristics differed
among intermittent and perennial sites.
To address H1 and H2, we designed a set of linear mixed effect models
(Table 1) to estimate the effects of environmental, hydrological and
connectivity variables on the responses of the instream leaf stock
quantity (dry mass), leaf litter quality (%N, %C and C/N),
invertebrate community structure (taxa richness, and leaf-shredder
abundance) and decomposition rates driven by invertebrates (coarse-mesh
bags; CK) and microorganisms (fine-mesh bags: FK; Table 1 & S2). To
test H3 and H4, in each model, we also included the riparian response
(equivalent to instream response; e.g. riparian CK for instream CK) and
its interaction with FP and DS as continuous predictor allowing us to
assess if relationship between instream and riparian processes changed
with flow intermittence and along the river network (Table 1). To test
H5, we included either richness or leaf-shredder abundance as predictor
and their interaction with FP and DS in CK models only, allowing us to
assess if the link between communities and functioning changed with
hydrological and spatial variables (Table 1 & S2). Polynomial
regressions were in some case attempted, when data exploration revealed
clear non-linear relationships.
We used a model-averaging approach to select the most influential set of
predictors and to calculate parameter estimates and significance for
each response variables. This way, all potential combinations of
predictors were assessed but only parameters of models with a Δ Akaike’s
information criteria corrected for small-sample sizes (AICc) <
2 from the model with the lowest AICc were selected and their estimates
averaged (Anderson 2008). Pseudo Nagelkerke’s R2 were
also calculated to estimate model goodness-of-fit (Nagelkerke 1991) and
averaged when more than one model was selected. To account for the
non-independence of samples from individual sites, we included site as a
random factors in each model. All predictors was centred around their
mean and scaled (i.e. divided by their standard deviation). In a few
cases, response variables were transformed, and outliers removed (Table
2) to ensure model assumptions of heteroscedasticity and normal
distribution of residuals. To reduce the disproportionate effect of a
few large DS values, DS was log-transformed in all models.
To address H3 and H4 in greater detail, we used another set of LMM to
determine if leaf stock characteristics, community composition and
decomposition differed between habitats (riparian vs. instream) and if
differences depended on site location (headwater vs. mainstem) and flow
regimes (intermittent vs. perennial). In these models, instream and
riparian responses were treated together as response variable and the
categorical variables habitat, location and flow regime and their
interactions were used as predictor. We also calculated pair-wise
Euclidean distances between 1) leaf resource characteristics (stock, C/N
and %N) and 2) ecosystem functions (CK and FK), and 3) Bray-Curtis
dissimilarity between invertebrate communities to assess and visualize
differences in resource, functioning and communities among habitats,
locations and flow regimes. We used non-metric multidimensional scaling
(NMDS) to visualize community variability and PCA ordinations to
visualize resource and ecosystem function variability. We used PERMANOVA
to further assess if leaf stock characteristics, community composition
and decomposition among sites differed depending on their location
(headwater vs. mainstem), habitats (riparian vs. instream) and flow
regime (intermittent vs. perennial). We used “season” as a strata
argument to constrain permutations among seasons. We used post-hoc tests
when interactions between habitat and location or habitat and flow
regime were significant to identify which group combinations differed
from each other.
To further explore H5 and determine how the relationship between
instream invertebrate-driven decomposition (CK) and invertebrate
community changes with flow regime and locations we used a Procrustes
analysis, comparing the similarity between two matrices in
multidimensional space (Peres-Neto and Jackson 2001). For this, we
compared the pair-wise Euclidean distance between CK values and the
Bray-Curtis community dissimilarity matrix among perennial,
intermittent, headwater and mainstem groups, separately.
We used R version 3.5.0 (R Core Team 2018) for all analyses, including
the packages vegan (Oksanen et al. 2019) for PERMANOVA and Procrustes
analyses, lme4 (Bates et al. 2015) and MuMIn (Bartoń 2019) for model
building and averaging, usdm (Naimi et al. 2014) for VIF analyses.