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.