Barriers decouple population dynamics of riverine fish, and asynchrony of subpopulations promotes stability within fragments

The spatial synchrony framework suggests that asynchrony among subpopulations in different branches of a river network should stabilise the metapopulation. However, how barriers affect this framework remains poorly understood. This is a significant knowledge gap given that population synchrony arises from dispersal and environmental similarity, both of which are influenced by barriers. Here, we empirically evaluated how barriers impact fish population synchrony and, subsequently, the associations between synchrony and metapopulation persistence, productivity, stability, and trajectory within fragments. We found that barriers demographically decouple populations by decreasing synchrony in brown trout (Salmo trutta) and Eurasian minnow (Phoxinus phoxinus) but not northern pike (Esox lucius), suggesting species-specific responses to fragmentation. Additionally, asynchrony had a stabilising portfolio effect on metapopulation stability at the fragment level that was statistically significant for S. trutta. Higher fragment synchrony also made S. trutta and P. phoxinus populations less stable. The impact of barriers on riverine fish population synchrony emphasises the need to include barriers in future studies on the causes and consequences of synchrony in rivers. That asynchrony stabilises populations in some riverine fishes suggests that conservation prioritisations should lie in restoring or retaining larger fragment sizes and higher branching complexity with intact connectivity.

Methods
Study species, data acquisition and river structure

Fish density data was acquired from Swedish Electrofishing RegiStry (SERS) at Swedish University of Agricultural Sciences (SLU). Electrofishing is a non-lethal fish sampling method mainly conducted in streams where it is possible to wade. Electric current (DC) is used to attract fish to swim towards a handheld anode where they are caught with a dipping net. It is an established and reliable method for quantifying fish density[1] in rivers, see the Swedish and European Standard[2] for detailed descriptions on the method. To evaluate generality and improve inference space, we used three model species – Salmo trutta (brown trout), Phoxinus phoxinus (Eurasian minnow), and Esox lucius (northern pike). These species represent three families of freshwater fish (Salmonidae, Esocidae and Cyprinidae) that are widely distributed in the northern hemisphere. Although coexisting in riverine habitats, these species require different habitats to complete their life cycles and sustain local populations. S. trutta is an obligate lotic spawner laying their eggs on gravel[3], E. lucius has a strong preference for spawning in lentic habitats with vegetated soft bottoms[4], whereas the P. phoxinus is intermediate, preferably spawning on hard bottoms in slow flowing lotic habitats[5].

Because of the extensive requirement on spatial and temporal coverage of sampling sites to answer our questions, an initial filtering was done where we removed all electrofishing sites with less than 6 sampling occasions. This allowed us to visually select 24 catchments with sufficient spatial and temporal coverage throughout Sweden (Fig. 1) with sites located between 56.62 to 66.83 degrees latitude and 12.10 to 21.82 degrees longitude. Temporal span of sampling data was 1951 to 2021, but 90% of the sampling occasions occurred after 1988. Vectored polylines of the rivers in these 24 catchments were downloaded from Lantmäteriet (Land Survey). The spatial distribution of barriers was downloaded from Swedish Meteorological and Hydrological Institute (SMHI). These three geographic information layers were integrated in R package riverdist[6], which was used to generate tables with the segment and vertex location of all barriers and electrofishing sites. Watercourse routes between all pairwise sites were recorded, and if the route crossed a segment and vertex of a dam, the two sample sites were assigned to different fragments. A matrix was then constructed with this information which clustered connected sites into fragments. The riverdist package also recorded watercourse distance, Euclidean distance, and flow-connectivity between all sites. See Fig. 1b for a visualisation of a catchment with sites split into fragments by barriers.

Of the total 1505 dam occurrences within in the study area, 224 were classified as removed and were thus filtered away before the spatial data extraction process. The metainformation in the registry for barrier occurrences was very patchy. Of the remaining 1281 dam occurrences, only 301 specify a dam height: the average height of those dam constructions was 6.4 (± 19 SD) meters. Information on whether barriers had the capacity to regulate the flow, or if they were barriers for fish, was also rarely specified. Rather than designing questions dependent on such patchy data, all remaining barriers were assumed to be barriers for fish.

Quantifying synchrony

Between sites

Synchrony between all sites in each catchment was quantified using Spearman’s correlations with at least 6 matching pairwise years for each species (n = 9998, n = 4653 and n = 2033 for S. trutta, P. phoxinus and E. lucius, respectively). Correlation analyses have been used to analyse spatial synchrony in several previous studies (Spearman’s correlations[7–9] and Pearson’s correlations[10]). We found that site-pairs that did not have any matching years with >0 density lead to negatively biased correlations, so to improve data quality all such site-pairs were excluded.

To evaluate whether the results and conclusions based on estimates of synchrony were robust or influenced to any important degree by the length of the time period, we also performed a sensitivity analysis using only site-pairs with at least 10 matching years (instead of at least 6 matching years). The results were qualitatively similar as no significances nor direction of effects were altered; these results are available in the Supplementary Material but we report the results based on analyses of at least 6 matching years in the main text below.

Because distances between site-pairs across fragments have the potential to be longer than those within fragments, the distance of between-fragment site-pairs was truncated to the maximum distance of site-pairs within fragments (for each catchment respectively). This resolves the statistical issue of unequal variances (and ranges) between groups, which can be problematic in models with interactions.

Linear mixed models were performed with R package lme4[11]. Because each “observation” in this dataset is based on information obtained from two sites, the data has a “pairwise” correlation structure (one row equals the synchrony, Euclidean distance, watercourse distance, etc. between two sites). Thus, each site (e.g., A, B, C, D) is used for the calculation of multiple synchrony values (e.g., A-B, A-C, A-D, etc) and that the same site sometimes occurs in different columns (e.g., A-B, B-C). To account for this non-independency issue, we used a multimembership random effect added using lme4-wrapper LmerMultiMember[12]. With this approach, a synchrony value can be seen as a group with the two sites as “members”, and some part of the variation can be partitioned to the sites across groups.

In light of recent studies[7], accounting for dendritic structure is critical in analyses of synchrony in river networks. To evaluate the influence of spatial distribution, flow-connectedness and fragmentation on site-pair synchrony, we constructed a linear mixed model with site-pair synchrony as the response variable, and the main effects of, and the three-way interaction between, ‘water course distance’, ‘flow-connectedness’ and ‘fragment border passing’ as explanatory variables. Random intercept for catchment identity and multimembership effect for site were included as random effects in the models. We analysed the data for each of the three species separately. To ensure that our model maintains the expected Type I error rate under the null hypothesis, we conducted a simulation-based evaluation outlined in the Supplementary Materials. Reassuringly, the proportion of p-values below the alpha level (0.05) was approximately 5% across all model terms except the intercept (Supplementary Fig. S1), confirming that our approach maintains appropriate Type I error rates.

To improve interpretability of results, we also performed models with flow-unconnected (n = 6632, n = 2561 and n = 1243 for S. trutta, P. phoxinus and E. lucius; Fig. 2a-c) and flow-connected site-pairs (n = 3366, n = 2092 and n = 790 for S. trutta, P. phoxinus and E. lucius; Fig. 2d-f) separately. Linear mixed models (LMMs) with ‘site-pair synchrony’ as response and ‘watercourse distance’, ‘fragment border passing’, and their interaction, as explanatory variables were performed for each species. Random intercept for catchment identity and multimembership effect for site were included as random effects. For flow-connected sites, the interaction effect between ‘watercourse distance’ and ‘fragment border passing’ was not significant (interaction term: t = 0.028) for pike so an additional model with only main effects was performed (Table 2; Fig. 2).

Scaling up site-synchrony to fragment-level

To quantify and compare synchrony at the level of fragments all site-pairs with different fragment identities were filtered away (leaving n = 4747, n = 2412 and n = 730 for S. trutta, P. phoxinus and E. lucius). The level of synchrony for each fragment (i.e., the metapopulation of the fragment) was then estimated by calculating the arithmetic mean[13,14]. This metric disregards spatial sampling density or dendritic structure so it is likely that fragments that contain sites that are more spaced out, or that are located in different branches, will have lower synchrony values, for example. While being dependent on the spatial sampling distribution, this metric does give an accurate view of the synchrony of the sampled sites.

Consequences of synchrony for populations

The portfolio effect

To quantify the stability (i.e., the portfolio effect response variable; PE) of each fragment metapopulation we chose the metric and calculated the PE as outlined by the practical guide by Anderson and colleagues[15]. The PE is defined as the ratio of the observed metapopulation CV to the CV the population as if it was one uniform population[15]. A value of 1.5 in this metric signifies a metapopulation that is 1.5 times more stable than a theoretically uniform population. Noteworthy is that “This metric does not address the benefit of increases in portfolio size (e.g., metapopulation size) itself”[15], implying that unequal sampling in fragments ought not to be a problem.

To choose the appropriate PE metric, we first plotted the mean and variance of each subpopulation time series across all metapopulations on log-log axes, which showed a linear (as opposed to non-linear) relationship, and retrieved species-specific slopes of z_strutta = 1.61, z_pphoxinus = 1.68, and z_elucius = 1.43. Since all of these z are different from z = 2 (an assumption of the average-CV PE), we used the more conservative PE metric, the mean-variance CV[15]. To make the estimation of the PE as accurate as possible, the data inclusion criteria were set at similar levels of the salmon example in Anderson et al. (2013)[15]: Only fragments with at least 4 subpopulations where each of the subpopulations had at least 10 sampled occasions were included. This resulted in a dataset of n = 40 fragments for trout, n = 33 fragments for minnow, and n = 26 fragments for pike. The mean ± SD number of sites per fragment was 13.7 ± 11.6 for S. trutta, 14.8 ± 12.4 for P. phoxinus, and 16.6 ± 13.32 for E. lucius.

We used linear mixed models to evaluate the relationship between fragment mean synchrony and the mean-variance portfolio effect response variable. Because some catchments might be represented by several fragments, we added catchment identity as a random effect to the model to account for any residual geographical non-independence.

Other population performance metrics

We also quantified five other population performances; (i) occurrence rate, (ii) mean log-density + 1, (iii) standard deviation in log-density, (iv) residual standard deviation (RSD), and (v) population trajectory slope, for all respective sites using all the available data (i.e., the entire time series of each site) (Table 1). These metrics thus have a larger sample size than for the portfolio effect above. The first three measurements were calculated at site level and then aggregated to fragment level by averaging the values from all sites with the same fragment identity. The last two measurements (residual standard deviation and population trajectory slope) were calculated through performing a model per fragment: the residual standard deviation is calculated, per fragment, from regressing all containing sites’ time series of log-density on years in interaction with site identity and using the residuals to calculate a population variability (Please see Supplementary Fig. S2 for visual explanation). An important advantage of this last way to estimate variability of sites in fragments is that it also detrends the data (i.e., preventing longer-term trends in data from being included in estimations of year-to-year variability). The mean fragment trajectory slope was used to evaluate if level of synchrony is associated with population trajectory. This resulted in a dataset of n = 101 fragments for trout, n = 70 fragments for minnow, and n = 62 fragments for pike. The mean ± SD number of sites per fragment was 7.35 ± 9 for S. trutta, 8.9 ± 10.2 for P. phoxinus, and 8.8 ± 10.8 for E. lucius.

A justification for using a range of different metrics is that they all capture slightly different ecological aspects of population performance[13]. Some, but not all, of the metrics of population performance were positively correlated, and the pairwise correlations differ somewhat between the three species (Supplementary Figs. S3, S4, S5)

To evaluate the association between synchrony and these population performance metrics, we performed linear mixed models with mean fragment synchrony as explanatory variable, and the five different measures of performance as responses, respectively, for each species. The test statistics were recorded in Table 3.

The Benjamini-Hochberg method was used to adjust the p-values for all population metrics, including the portfolio effect in the previous section.

Does adding information on barriers improve a previous framework?

To investigate if adding barriers adds significantly to the predictive power of the recent development of a “geography of spatial synchrony in dendritic river networks”[7], we replicated the fluvial synchrogram variables defined therein in our dataset. The synchrograms developed by Larsen and colleagues in 2021[7] offer insights into the effects of both hydrological connectivity and upstream dependence among populations, while also capturing relationships among populations across adjacent tributaries and the broader landscape context. They calculate Euclidean distance (dE) and watercourse distance (dW) and their quotient (dE/dW) between all pairwise sites and subsequently recognise that sites can be separated by four functionally different distances: D1, D2, D3 and D4. The authors present four major expectations based on these four combinations of pair-wise distances:

D1: Populations with small and similar dE and dW values are likely on the same network branch, displaying high synchrony. They are close both in terms of physical distance and water flow, suggesting a strong potential for interaction and shared dynamics.
D2: Populations with large and equal dE and dW values are likely located on the same branch but are more distant. They are still expected to show intermediate synchrony, influenced by a combination of dispersal and a Moran effect.
D3: Populations with small dE but much larger dW values are likely situated in separate but nearby branches. They should display intermediate synchrony, primarily driven by a Moran effect. This suggests that even though they may not be directly connected by water flow, there might still be shared environmental influences.
D4: Populations with large dE and much larger dW values are positioned on distant and separate branches. They are expected to show the lowest degree of synchrony due to their significant spatial and hydrological separation.

We set out to explore whether and how the inclusion of barriers enhances the predictive capacity of this framework within our dataset. Euclidean distance, watercourse distance and their quotient (dE/dW) were calculated for all sites. The classification was then done river-catchment-wise, so that the delineation into short or long dE, and high or low dE/dW, was set at the catchment-wise median values, and all pairwise observations in the present study were subsequently categorised into the four classes of distance (D1, D2, D3, D4). Two models were then performed per species, one that modeled the synchrony main effect of the synchrogram category, and another model that included main effect of synchrogram and its interaction effect with fragment border crossing (1/0). Akaike Information Criterion (AIC) was then used to compare the parsimony of these two (per species) models[16]. We hypothesised that closely located sites (those separated by D1) would show the largest divergence in synchrony when either separated or not separated by barriers. Because of the nature of synchrony decay with distance, we hypothesised that differences between models would be smaller for D2 through D4.

Softwares and packages used

All data manipulations, calculations and analyses were done in R language[17]. Widely used packages include those in the tidyverse ecosystem[18], lme4[11], PerformanceAnalytics[19], and riverdist[6]. We used the R package car[20] with the Anova function to estimate p-values. We used the R package performance[21] with the r2 function to estimate conditional and marginal R2 values. Correction for multiple testing was done using Benjamini-Hochberg with the p.adjust function from stats[17] package.

References

1.                   Bohlin T, Hamrin S, Heggberget TG, Rasmussen G, Saltveit SJ. 1989 Electrofishing — Theory and practice with special emphasis on salmonids. Hydrobiologia 173, 9–43. (doi:10.1007/BF00008596)

2.                   2003 SS‐EN 14011:2003 Water quality – Sampling of fish with electricity.

3.                   Jonsson B. 1989 Life history and habitat use of Norwegian brown trout (Salmo trutta). Freshwater Biology 21, 71–86.

4.                   Nilsson J, Engstedt O, Larsson P. 2014 Wetlands for northern pike (Esox lucius L.) recruitment in the Baltic Sea. Hydrobiologia 721, 145–154.

5.                   Frost WE. 1943 The natural history of the minnow, Phoxinus phoxinus. The Journal of Animal Ecology , 139–162.

6.                   Tyers M. 2022 riverdist: River Network Distance Computation and Applications. See https://CRAN.R-project.org/package=riverdist.

7.                   Larsen S et al. 2021 The geography of metapopulation synchrony in dendritic river networks. Ecology Letters 24, 791–801. (doi:10.1111/ele.13699)

8.                   Vendrametto Granzotti R, Agostinho AA, Bini LM. 2022 Drivers and spatial patterns of population synchrony of fish species in a floodplain. Freshwater Biology 67, 857–872. (doi:10.1111/fwb.13886)

9.                   Chevalier M, Laffaille P, Grenouillet G. 2014 Spatial synchrony in stream fish populations: influence of species traits. Ecography 37, 960–968. (doi:10.1111/ecog.00662)

10.                 Terui A, Ishiyama N, Urabe H, Ono S, Finlay JC, Nakamura F. 2018 Metapopulation stability in branching river networks. Proc. Natl. Acad. Sci. U.S.A. 115. (doi:10.1073/pnas.1800060115)

11.                 Bates D, Mächler M, Bolker B, Walker S. 2015 Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software 67, 1–48. (doi:10.18637/jss.v067.i01)

12.                 Paridon J van, Bolker B, Alday P. 2022 lmerMultiMember: Multimember Random Effects Wrapper for lme4.

13.                 Yang Q, Hong P, Luo M, Jiang L, Wang S. 2022 Dispersal Increases Spatial Synchrony of Populations but Has Weak Effects on Population Variability: A Meta-analysis. The American Naturalist 200, 544–555. (doi:10.1086/720715)

14.                 Abbott KC. 2011 A dispersal-induced paradox: synchrony and stability in stochastic metapopulations. Ecology Letters 14, 1158–1169. (doi:10.1111/j.1461-0248.2011.01670.x)

15.                 Anderson SC, Cooper AB, Dulvy NK. 2013 Ecological prophets: quantifying metapopulation portfolio effects. Methods Ecol Evol 4, 971–981. (doi:10.1111/2041-210X.12093)

16.                 Zuur AF, Ieno EN, Smith GM, others. 2007 Analysing ecological data. Springer.

17.                 R Core Team. 2022 R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. See https://www.R-project.org/.

18.                 Wickham H et al. 2019 Welcome to the tidyverse. Journal of Open Source Software 4, 1686. (doi:10.21105/joss.01686)

19.                 Peterson BG, Carl P. 2024 PerformanceAnalytics: Econometric Tools for Performance and Risk Analysis. See https://CRAN.R-project.org/package=PerformanceAnalytics.

20.                 Fox J, Weisberg S. 2019 An R Companion to Applied Regression. Third. Thousand Oaks CA: Sage. See https://www.john-fox.ca/Companion/.

21.                 Lüdecke D, Ben-Shachar MS, Patil I, Waggoner P, Makowski D. 2021 performance: An R Package for Assessment, Comparison and Testing of Statistical Models. Journal of Open Source Software 6, 3139. (doi:10.21105/joss.03139)