The geography of metapopulation synchrony in dendritic river networks.

Dendritic habitats, such as river ecosystems, promote the persistence of species by favouring spatial asynchronous dynamics among branches. Yet, our understanding of how network topology influences metapopulation synchrony in these ecosystems remains limited. Here, we introduce the concept of fluvial synchrogram to formulate and test expectations regarding the geography of metapopulation synchrony across watersheds. By combining theoretical simulations and an extensive fish population time-series dataset across Europe, we provide evidence that fish metapopulations can be buffered against synchronous dynamics as a direct consequence of network connectivity and branching complexity. Synchrony was higher between populations connected by direct water flow and decayed faster with distance over the Euclidean than the watercourse dimension. Likewise, synchrony decayed faster with distance in headwater than mainstem populations of the same basin. As network topology and flow directionality generate fundamental spatial patterns of synchrony in fish metapopulations, empirical synchrograms can aid knowledge advancement and inform conservation strategies in complex habitats.

INTRODUCTION

Metapopulation synchrony, the coherent temporal dynamics in the abundance of spatially separated populations, has been observed across a wide range of taxa and is now considered a fundamental property of metapopulations (Liebhold et al. 2004; Wang et al. 2019). Spatial synchrony has well‐recognised implications for the long‐term persistence of species and ecosystem stability (Heino et al. 1997; Wilcox et al. 2017; Erős et al. 2020). Synchrony can decrease opportunities for demographic rescue among populations leading to higher local extinction risks when facing environmental change, while asynchrony may lead to longer‐term stability (Loreau et al. 2003). Spatial synchrony can arise from a combination of intrinsic and extrinsic mechanisms, including dispersal among connected populations, community processes such as predator–prey interactions, and spatially correlated environmental factors, also known as the Moran effect (Moran 1953; Grenfell et al. 1998; Liebhold et al. 2004). Disentangling these mechanisms remains challenging, except in simulation studies (e.g. Wang et al. 2019) or in specific geographical contexts where dispersal between populations is prevented (e.g. Tedesco et al. 2004).

Scientific investigations have recently focused on examining the spatial dimensions of synchrony (Walter et al. 2017), predominantly relying on exploring the decay of synchrony with geographic distance among populations (Hanski & Woiwod 1993; Sutcliffe et al. 1996; Abbott 2007; Jarillo et al. 2018). Despite important insights gained, this approach typically assumes isotropic changes in synchrony across homogeneous landscapes and thus potentially overlooks complex spatial patterns. Such challenges are acutely manifested in river ecosystems where both metapopulations and environmental dynamics reflect the topology, directionality and connectivity of the network. The isotropic assumption is clearly violated in river systems whose network geometry has profound influence on instream physical, ecological and evolutionary processes. In fact, the unique attributes of hierarchical dendritic structures relative to linear or random networks have long been recognised (Campbell Grant et al. 2007; Erős et al. 2012; Filipe et al. 2017; Erős & Lowe 2019), with models linking the branching connectivity of river‐like networks with greater metapopulation persistence (Fagan 2002; Sarhad et al. 2014; Ma et al. 2020). This occurs because river geometry (hereafter ‘network topology’) and the unidirectional water flow promote asynchronous dynamics among populations, thus favouring species persistence (Yeakel et al. 2014; Tonkin et al. 2018).

Riverine metapopulation dynamics have been investigated extensively in recent theoretical and experimental works (e.g. Yeakel et al. 2014; Bertuzzo et al. 2015; Altermatt & Fronhofer 2018; Anderson & Hayes 2018). Yet, assessments of how network topology influences the spatio‐temporal dynamics and synchrony of metapopulations in natural river systems are strikingly rare, despite the conservation implications of such studies (Moore et al. 2015; Terui et al. 2018). Here, we address this issue using freshwater fishes of Europe as an exemplar. Focusing on the influence of dispersal and Moran effect on synchrony within river networks, we first derive theoretical expectations regarding spatial aspects of riverine metapopulation synchrony based on principles of fluvial variography (i.e. geostatistics accounting for spatial dependencies within dendritic networks; Peterson et al. 2013) and network theory (Erős & Lowe 2019). We further support the key expectations using simulations from a spatially explicit metacommunity model applied to river networks. This provides the theoretical basis of the ‘fluvial synchrogram’ concept as a graphical exploratory tool. Next, we confront expectations from theoretical synchrograms using empirical synchrony estimates between > 34 000 pairs of fish species populations from an extensive abundance time‐series data across Europe (Comte et al. 2021). Although not designed to provide an explicit quantification of the different synchrony mechanisms, our analytical framework elucidates the emerging patterns of synchrony manifested from spatial variations in the underlying mechanisms. We thus articulate an empirically driven ‘geography of synchrony’ (Defriez & Reuman 2017; Walter et al. 2017) within river basins. Our results demonstrate the existence of fundamental aspects of population synchrony that are directly predictable from network topology and in‐stream connectivity. Network branching complexity appears to buffer synchronous dynamics, and headwater fish populations display a faster decay of synchrony with distance compared to those in mainstem habitats. By allowing prediction of synchrony patterns even if empirical population time‐series are not available, our findings have implications for the persistence and management of populations in complex habitats, such as streams and rivers that support high levels of diversity and are among the most threatened ecosystems globally (Tickner et al. 2020).

MATERIAL AND METHODS

Theoretical expectations of metapopulation synchrony and the fluvial synchrogram concept

Here, we derive key expectations regarding the geography of fish metapopulations synchrony within river networks (Figs 1 and 2). To visualise these expectations, we propose the concept of fluvial synchrogram (Fig. 1), and present basic theoretical synchrograms using simulated abundance time‐series from a spatially explicit dynamic metacommunity model (Ryser et al. 2019). Importantly, the simulation model was used to illustrate the key expected patterns in synchrony, but not to resolve the role of network complexity per se on fish metapopulation dynamics. As such, we did not perform multiple simulations of river network configurations as this has been covered by previous investigations (e.g. Yeakel et al. 2014; Anderson & Hayes 2018). The metacommunity model corresponds to a food chain composed of a basal resource (e.g. algae) supported by nutrients, an herbivore (e.g. macroinvertebrate) and a fish predator, where fish and herbivores are able to disperse between patches of a river network generated with optimal channel network methods, implemented in R by the OCNet package (Carraro et al. 2020). Using metacommunity simulations allowed us to explore the extent to which spatial patterns of synchrony could emerge directly from dispersal and the unique structure of dendritic river networks (see Appendix S1 in Supplementary Information for description of the metacommunity model).

Figure 1

(a) Hypothetical river network with five populations (labelled a to e) whose geographic separation can be measured as Euclidean (orange dashed lines) and hydrological distance (blue watercourse line). In addition, flow‐connected locations can be identified in which water flows from one to the other (e.g. purple dashed line connecting a to e, but not a to b, c or d). We refer to this distance as flow‐connected. (b) Theoretical ‘Fluvial synchrogram’ derived from simulated metacommunity fish abundance time series (see text and Appendix S1), depicting the decay of pairwise population synchrony over the three types of spatial distances. (c) Theoretical ‘3D synchrogram’ displayed as 2D contour GAM model. The synchrony among pairs of spatially separated populations (small dots) can be plotted as a function of actual Euclidean distance (x‐axis) against the ratio of Euclidean (dE) to Watercourse (dW) distance (dE/dW, y‐axis). Four major types of pairwise distance combinations can be identified on the 3D synchrograms (D1 to D4), as shown in (a) and (c).

Figure 2

Relationship between network branching complexity and the ratio of Euclidean to watercourse distance (dE/dW) between populations, represented by colored dots over the networks (a). For populations distributed over ‘simple’ less branching basins, the mean pairwise dE/dW is expected to be higher; that is closer to the 1:1 line, as indicated by the green dashed lines (b). Conversely, populations distributed in more branching ‘complex’ networks should be separated, on average, by lower dE/dW distances (as the relative dW increases). A geometric demonstration of these patterns is given in Fig. S1

We further derive expectations on how synchrograms are likely to vary within stream networks – from headwater to mainstem reaches – and between networks of different branching complexity.

The fluvial synchrogram concept

In order to capture the inherent spatial complexity of river habitats, dedicated geostatistical approaches are needed to reveal spatial structures over both Euclidean and watercourse dimensions, while concurrently accounting for flow directionality (Peterson et al. 2013; Zimmerman & Ver Hoef 2017). Geographical separation between sampling locations (populations) on a river network can be measured by three types of distance (Fig. 1a): Euclidean, watercourse (hydrological distance) and flow‐connected (hydrological distance between locations where water flows from one to the other). We draw from fluvial geostatistics and introduce the ‘fluvial synchrogram’ as a graphical exploratory tool to depict the decay of pairwise population synchrony with these distances. To illustrate the expected patterns, we derived theoretical synchrograms using simulated fish abundance time‐series on a 100‐segment river network (Fig. 1b; Appendix S1).

Synchrograms among flow‐connected populations can inform on the effect of hydrological connectivity and upstream dependence between populations. Conversely, synchrograms based on watercourse (including locations not directly linked by water flow) and Euclidean distances describe relationships between populations in adjacent tributaries and across the wider landscape context respectively (Ver Hoef et al. 2006; Larsen et al. 2019). Patterns are then expected to differ among geographical distance types (Fig. 1b). In the case of obligate aquatic biota, such as fish, we expect a steeper decay of synchrony over Euclidean than watercourse distance, the latter reflecting the actual connectivity perceived by individuals (Olden et al. 2001). In addition, over relatively short distances, populations directly linked by water flow are likely to experience similar environmental dynamics and higher dispersal rates, leading to higher synchrony than equally distant populations not connected by direct water flow (Tonkin et al. 2018). Over larger distances, however, this pattern may revert, as flow‐connected populations will necessarily inhabit reaches at the opposite ‘margin’ of the network (e.g. headwater to river mouth) and thus display lower synchrony. In fluvial synchrograms, this is depicted as a relatively high synchrony at short distances and steep decay of synchrony for flow‐connected populations (Fig. 1b).

Adding a third dimension to the fluvial synchrogram (3D synchrogram), we can calculate the ratio of Euclidean to watercourse distance (hereafter: d) to describe the relative location of any given pair of spatially separated populations within the network. Thus, the dE/dW metric represents a dimensionless measure of functional (as opposed to absolute) spatial separation as perceived by populations, which can be plotted for each pair‐wise comparison. Values of dE/dW close to 1 imply that populations are likely located on the same branch or segment of the network, whereas small values suggest that populations are located on different branches (Fig. 1a and c). Low values of dE/dW can also indicate separate locations along the same branch of highly meandering rivers when the absolute Euclidean distance is small. According to this rationale, the synchrony between pairs of locations throughout a river network can be examined with respect to their values on the two‐dimensional surface defined by dE/dW and Euclidean distances. We note that by including Euclidean distances on both axes of the 3D synchrogram, a spurious shape may occur due to the non‐independent variable formulation (Pearson, 1897). However, we were not interested in the statistical relationships between dE/dW and Euclidean distances, but rather on how synchrony varies spatially with respect to both. Second, including the Euclidean distances on the x‐axis allows anchoring the dimensionless dE/dW into a spatially explicit context, where the distribution of site pairs is obviously influenced by the inherent structure of river networks (e.g. site pairs are unlikely to be located at the top right corner of the 3D synchrogram where Euclidean distances are large, but dE/dW is small; see Fig. 1c). The spatial synchrony patterns over the 3D synchrogram thus provide an exploratory tool to appraise the underlying mechanisms.

Synchrograms can then help appraise the relative importance of dispersal and Moran effect in determining observed synchrony. Indeed, watercourse distance and flow directionality are inherently related to dispersal probability. In contrast, Euclidean distances are most likely reflecting the probability of a Moran effect, although we note that environmental autocorrelation can occur both along the Euclidean and watercourse dimensions (see Discussion). Four major expectations can be formulated according to combinations of pair‐wise distances on the 3D synchrogram (D1 to D4; Fig. 1a and c). Populations separated by D1 (i.e. small and similar dE and dW) are likely located on the same network branch and expected to display high synchrony, being proximate in terms of both distance types and similarly influenced by dispersal and Moran effect. Conversely, populations separated by D4 (i.e. large dE and much larger dW) are likely positioned on distant and separate branches and thus expected to display the lowest degree of synchrony. Populations separated by D3 (i.e. small dE but much larger dW; located in separated but nearby branches) are expected to display intermediate degree of synchrony, which should be primarily influenced by a Moran effect. Finally, synchrony over D2 (i.e. large and equal dE and dW; distant populations, but likely located on the same branch) is also expected to display intermediate degree of synchrony resulting from a combination of both dispersal and a Moran effect. These expectations were well exemplified by the theoretical 3D synchrogram in Fig. 1c fitted using a tensor‐product Generalised Additive Model (GAM) on simulated metacommunity fish abundance time‐series (described above).

Differences within networks

As the rate of dispersal among populations vary predictably across the network hierarchy, synchrograms are expected to differ between headwaters and mainstem reaches. According to river network theory, mainstem reaches are characterised by higher dispersal rates compared to more isolated headwaters due to their central position in the network, which integrates movements of organisms within and between branches (Brown & Swan 2010; Erős et al. 2012). In addition, local environmental conditions are also likely to vary more rapidly with distance in headwaters (e.g. Clarke et al. 2008). Therefore, we predict a faster decay of synchrony with distance between populations in low‐order compared to high‐order river segments.

Differences between networks

Basin shape and network topology inherently influence the position of population pairs over the 3D synchrogram. Populations in complex branching networks will be, on average, separated by lower dE/dW distances compared to simpler branching networks (Fig. 2). As such, basin‐level dE/dW is a measure of branching complexity that reflects the degree of branching separation experienced by populations, and it is expected to directly influence metapopulation dynamics. Hence, the average dE/dW in a basin varies according to the distribution of populations and effectively captures the realised network complexity (Fig. S1). Moreover, although the data used to present our framework come from monitoring programmes whose sampling locations were not designed to reflect network branching, all available stream orders were well represented across the basins (Fig. S2).

In addition, highly branching networks are often characterised by relatively larger tributaries and higher proportions of geomorphically significant confluences (Benda et al. 2004), which promote physical heterogeneity along the river network (e.g. less predictable sediment size, and variability of reach‐scale habitat features and flow regimes). Given these structural constrains, we expect populations distributed in networks with lower dE/dW distances to display less synchronous dynamics.

Deriving empirical fluvial synchrograms

We compared theoretical synchrograms to the geography of fish metapopulation synchrony within 58 river basins throughout Europe (Fig. 3a) using long‐term stream fish time‐series (≥ 10 years; Comte et al. 2021). Surveys used standardised protocols through time, and sampling occurred during low flow periods (summer‐autumn). The basins were selected to include at least eight sites (i.e. stream reaches; median = 16; range = 8–63). Only species occurring in at least 80% of sampling events (median = 7 species/basin; range = 1–23 species/basin) were included to limit the influence of zeros and low means when calculating synchrony (e.g. Chevalier et al. 2014). Overall, the dataset contained 1150 sites, 48 species and > 34000 pairs of fish population time series.

Figure 3

(a) Location of sites (n = 1150) and basins (58) used in the study. (b) Empirical synchrograms showing the decay of synchrony as separate exponential fits for watercourse (continuous line), Euclidean (dot‐dashed) and flow‐connected (dashed; directly linked by water flow) distances. Confidence intervals were estimated with Monte Carlo uncertainty propagation. To aid visualisation, the mean synchrony values for 30 distance bins is also shown for watercourse (indicated with ‘+’) and Euclidean (‘x’) distances, where each cross includes c.1100 population pairs.

Analytical approach

For each species in a given river basin, we expressed pairwise population synchrony as Spearman correlations between the pairs of abundance time‐series. We extracted the Euclidean and watercourse distance separating each pair of populations, while also distinguishing population pairs directly linked by water flow (flow‐connected), according to HydroRIVERS (Lehner et al. 2008). In the latter case, non‐flow‐connected populations were excluded from the model (Zimmerman & Ver Hoef 2017). These data were combined to construct the empirical fluvial synchrogram depicting the decay in synchrony over the Euclidean, watercourse and flow‐connected distances. The synchrony‐decay was fitted with an exponential model:where y represents the pairwise synchrony, a the intercept, b the decay, and distance is square‐root transformed. Confidence intervals for the different nonlinear fits were estimated using Monte Carlo uncertainty propagation as implemented in the R package ‘propagate’ (Spiess 2018). The exponential model with square‐root transformed distances provided better fits to the data compared to linear and non‐transformed models (delta‐AIC > 10 in all cases; Table S1).

From the synchrograms we derived the decay (b), and the short‐scale synchrony, corresponding to the 1st percentile of the empirical distance distribution (c. 1 km). This 1‐km synchrony was used to indicate the maximum expected synchrony for each distance type and was considered as a more interpretable parameter than the model intercept (synchrony at 0‐distance).

We then used null models to evaluate differences in model parameters between distances. We created null synchrograms by randomly shuffling (999 times) the ‘labels’ of distance types (i.e. Euclidean, watercourse, flow‐connected) within basins, while maintaining the observed distances and synchrony estimates. This procedure broke the correspondence between synchrony and distance types within each basin, but preserved the overall data structure. We compared the observed differences in decay and 1‐km synchrony between distance types with those from null synchrograms. Differences were expressed as z‐score:with |z‐score| > 1.96 indicating significant difference in parameters between the examined distance types (Gotelli 2001).

To display the 3D synchrogram, we deployed a flexible local polynomial smoothing (LOESS) with a span = 0.75, using Euclidean distance and dE/dW as predictors. To account for the complex structure of the data (i.e. multiple observations per species and basins) we also fitted synchrograms including random species and basin effects.

To examine how synchrograms varied across the network hierarchy, we fitted synchrograms based on watercourse distances separately for populations located in low‐order (i.e. headwaters; Strahler order 1‐3, 12 934 population pairs, 37 species), high‐order (i.e. downstream; order 4–7; 11106 population pairs, 40 species) and mixed order (between high‐ and low‐order; 10 452 population pairs, 39 species) river segments. Differences in synchrogram parameters between stream order groups were tested using null‐models that randomly reshuffled synchrony values across stream order within basins, and were then expressed as z‐scores as in eqn (2).

To assess the overall effect of branching complexity on metapopulation synchrony, we calculated the mean dE/dW between all pairs of sites within each basin. To control for the influence of basin size, mean dE/dW was regressed against basin area and residuals were used in the analyses. We then fitted synchrograms (over the watercourse dimension) for each species within each basin separately. For each basin, synchrony decay and 1‐km synchrony were expressed as the mean among species weighted by species‐specific overall abundance (Wang et al. 2019). We also computed basin‐level synchrogram parameters using log‐transformed species abundances to reduce the influence of highly abundant species. The influence of basin‐level dE/dW on the decay of synchrony and its values at 1‐km, was tested using a simple linear regression. A quantile regression was further used to examine whether network branching specifically influenced the lower or upper distribution of synchrogram parameters.

RESULTS

Empirical fluvial synchrograms

In accordance with the theoretical synchrograms, fish metapopulation synchrony decayed more rapidly over Euclidean distance than watercourse distance (Fig. 3b; z‐score difference = 7.35), which also displayed a slightly higher, albeit non‐significant (z‐score = −1.28), short‐scale synchrony (Table 1).

Table 1

Parameters from synchrograms (SE) including synchrony estimates at 1‐km distance (1‐km synch) and decay for watercourse, Euclidean and flow‐connected distances, and between low‐, high‐ and mixed‐order stream pairs (over watercourse distance)

	1‐km synch	Decay
Watercourse	0.21 (0.006)	−0.065 (0.003)
Euclidean	0.23 (0.006)	−0.122 (0.005)
Flow‐connected	0.26 (0.008)	−0.085 (0.006)
Low‐order	0.25 (0.008)	−0.074 (0.004)
High‐order	0.19 (0.007)	−0.028 (0.004)
Mix‐order	0.10 (0.01)	−0.045 (0.01)

Flow‐connected populations displayed the highest short‐scale synchrony (z‐score difference from watercourse and Euclidean = −3.55 and −6.65 respectively) with higher synchrony up to c. 50 km separation (Fig. 3b). However, the decay in synchrony along flow‐connected distances was steeper than along watercourse distances (z‐score = 3.3), but similar to the decay over the Euclidean dimension (z‐score = 1.13). The synchrogram fitted using random basin and species effects displayed similar patterns (Fig. S3).

Agreeing with theoretical expectations, the 3D fish metapopulation synchrogram (Fig. 4) corresponded with the one obtained from the simulations (see Fig. 1c). The highest synchrony was observed between fish populations separated by short and comparable Euclidean and watercourse distances (D1 distance combination in Fig. 4; populations likely located on the same network branch). Conversely, the lowest synchrony was found over D4 distance combination, as expected for populations separated by large Euclidean distances and inhabiting different branches. We also observed high synchrony between fish populations over short Euclidean distances, albeit separated by relatively longer watercourse distances (D3 in Fig. 4), suggesting a Moran effect driving mechanism in adjacent headwaters. Similarly, the comparable level of synchrony predicted over D2 and D3 suggests that as the relative functional separation declines (i.e. ‘straighter’ watercourse distances with dE/dW approaching 1), similar levels of synchrony manifest despite larger physical distances between populations. This implies that, over equal Euclidean separation, fish populations in more branching networks (i.e. with lower dE/dW) are likely to display less synchronous dynamics.

Figure 4

3D synchrogram modelled as 2D contour LOESS illustrating fish metapopulation synchrony over the plane defined by the ratio Euclidean/watercourse distance (dE/dW) against Euclidean distance. The position of four major types of pair‐wise distance combinations over the plane are also shown (D1 to D4, as presented in Fig. 1). See Figs. S4 and S5 for an alternative formulation of the 3D synchrograms using a tensor‐product GAM that include random basin and species effects.

Synchrograms differed when separately modelling fish populations within low‐order, high‐order and mixed‐order river segments. Synchrony decayed faster for populations occupying low‐order streams compared to mainstem populations (Fig. 5; z‐score = 2.98). Also, the short‐scale synchrony for low‐order streams was slightly higher than for high‐order streams, with a marginally significant z‐score = −1.96, reflecting higher synchrony between headwater populations within c. 25 km stream distance. The lowest level of synchrony occurred between populations in low‐ and high‐order reaches (mixed‐order; Table 1) with lower 1‐km synchrony than high‐order populations (z‐score = 5.56), albeit similar decay with distance (z‐score=−1.45).

Figure 5

Synchrograms fitted separately for fish populations within low‐ (headwaters), high‐ (mainstem) and mixed‐order stream reaches. Confidence intervals were estimated with Monte Carlo uncertainty propagation.

Basin‐level 1‐km synchrony (weighted by overall species abundance) increased significantly as network complexity declined (i.e. larger dE/dW, independent of basin size (Fig. 6a; P = 0.03). This occurred despite observing similar basin‐level decay of synchrony with distance (Fig. 6b). This suggests that network branching complexity, as perceived by populations, can influence the maximum levels of synchrony within a network. This is further supported by the significant increase in the upper quantile of 1‐km synchrony in less branching networks (grey line in Fig. 6; quantile regression at q = 0.8; P < 0.001). Qualitatively similar results were observed when weighting basin‐level synchrony by the log‐transformed species abundance (Fig. S6).

Figure 6

Relation of basin‐level short‐scale synchrony (1‐km synch; a) and decay (b) over the watercourse dimension (weighted by overall species abundance) with the mean ratio of Euclidean – Watercourse distance between fish populations (mean dE/dW; controlled for basin area). Dashed line indicates non‐significant relationship. Larger dE/dW values represent populations distributed, on average, over less branching networks. The significant upper quantile regression at q = 0.8 is also shown (grey line; A).

DISCUSSION

In this study, we introduce the fluvial synchrogram as a conceptual framework and exploratory tool to examine the geography of metapopulation synchrony within riverine dendritic networks. Using simulations and empirical data, we showed that fundamental spatial aspects of fish metapopulation synchrony can be predicted from stream network topology and connectivity. Fluvial synchrograms revealed that spatial patterns in fish temporal dynamics resembled those typical of spatial dynamics within dendritic networks. The overall shape of the fluvial synchrogram mirrored the spatial autocorrelation patterns often revealed by empirical variograms in river networks, known as ‘Torgegrams’ (Peterson et al. 2013; Zimmerman & Ver Hoef 2017). For instance, for most in‐stream variables and processes, flow‐connected locations often display higher spatial correlation (analogous to higher synchrony here) than those not directly linked by flow; however, this relationship can often reverse at larger distances (e.g. Fig. 7 in Peterson et al. 2013). The role of spatial processes and network topology in driving synchrony was further supported by the similarity of the empirical synchrograms with patterns derived from simulations that emerged from dispersal within dendritic networks.

The observed empirical patterns also indicate that network structure, such as the relative proportion of low‐ and high‐order segments and network branching complexity, can directly regulate the degree of synchrony among fish populations. In particular, the steeper decline of synchrony among headwater populations implies a stronger effect of isolation by distance. This is consistent with stream network theory, where connectivity and dispersal are predicted to be higher over the central nodes (mainstem) than the marginal nodes (headwater) of the network (Finn et al. 2007; Brown & Swan 2010; Erős & Lowe 2019). This may reflect a general tendency for fish species to disperse less along headwaters than mainstem segments (e.g. Radinger & Wolter 2014), and the fact that environmental conditions could change more rapidly and unpredictably along headwater reaches (Clarke et al. 2008), thus limiting synchrony. However, synchrony remained high between headwater populations up to c. 25 km watercourse distance. It is possible that over such relatively short distances, environmental conditions change less along low‐order than high‐order reaches, where the presence of tributaries and human pressure could have important influences. The empirically derived synchrony patterns across the network hierarchy are also coherent with the theoretical expectations shown in the 3D synchrograms. Fish populations in headwaters were, on average, separated by lower dE/dW (0.45) compared to those in the mainstem (0.57). That is, they would be primarily located on the upper half of the 3D synchrogram in Fig. 4, where synchrony is lower. As expected, the lowest synchrony occurred between headwaters and mainstem populations (mixed‐order), as these are separated by equally low dE/dW (0.45), but are also likely to experience markedly different environmental conditions.

Our results support the growing recognition that riverine metapopulations can be buffered against synchronous dynamics as a direct consequence of network branching (Yeakel et al. 2014; Terui et al. 2018). Interestingly, we found that network branching appeared to limit the maximum expected synchrony (i.e. at short distances) rather than its rate of decay. Our measure of realised network branching – mean dE/dW among populations – is a proxy of dispersal limitation across the basin that reflects the actual distribution of the populations over the network, rather than the overall network shape.

However, basin shape and network branching are also related to the physical diversity of channel and riparian conditions (Benda et al. 2004). This is supported by the correlation between basin level dE/dW and circularity ratio (Fig. S7; P < 0.01; R 2 = 0.15), a network‐wide measure of basin compactness, which reflects confluence‐mediated habitat heterogeneity (Benda et al. 2004). This implies that both dispersal limitation and environmental‐related mechanisms combine to buffer synchrony in branching river networks. Basin shape and branching geometry are controlled by factors such as climate and geology across evolutionary times (Seybold et al. 2017; Yi et al. 2018), suggesting that the overall propensity for synchrony in riverine metapopulations may inherently depend on larger scale processes and biogeographic settings.

This study sought to elucidate the geography of metapopulation synchrony in river networks emerging from the spatial variation in the underlying mechanisms. As such, it was not designed to assess the contribution of species‐specific dispersal or the influence of specific environmental factors (i.e. Moran effect). That is, the modelled synchrograms did not account for species and basin identity; however, patterns based on a mixed‐modelling framework where species and basins were included as random factors were similar (Figs S2 and S4), supporting the generalities of our findings and the utility of synchrograms as exploratory tools. For instance, patterns of metapopulation synchrony could be compared with those obtained from environmental data to appraise the contribution of Moran effects to population dynamics (e.g. Defriez & Reuman 2017). The advantage of using fluvial synchrograms is that they can indicate whether Moran effects may be related to factors operating over the Euclidean (e.g. climate, fire) or the watercourse (e.g. water quality, flow‐regime) dimension. The importance of defining appropriate dispersal pathways and connectivity matrix beyond the Euclidean dimension has clearly emerged from recent studies of population synchrony (e.g. Anderson et al. 2017, 2019; Lopes et al. 2018; Zanon et al. 2018). This is particularly evident in constrained habitats where dispersal is directional or for organisms relying on wind‐borne dispersal (Vindstad et al. 2019). Since fluvial synchrograms rely on multiple types of pairwise distances, several hypothesised drivers of synchrony could be examined using matrix regressions, which can specifically control for alternative mechanisms (Walter et al. 2017).

Similarly, exploring differences in synchrograms fitted for different species (or functional groups) within the same basin represents an exciting research avenue into the relative importance of dispersal (i.e. different patterns over the watercourse and flow‐connected distance) versus environmental factors (i.e. different patterns over the Euclidean dimension) in driving synchrony patterns.

Our understanding of how dendritic habitats influence species dispersal, persistence and diversity, has grown steadily over the last 20 years. Yet, the implications of river network structure for metapopulation synchrony and stability have been explored only recently and primarily through theoretical approaches. Overall, our results are consistent with theoretical work in showing that the geometry of river networks can inherently promote the persistence of species by favouring spatial asynchronous dynamics among localities (Yeakel et al. 2014; Anderson & Hayes 2018). This has important conservation implications. The spatially explicit nature of synchrograms allows estimating the distance at which synchrony falls below any desired thresholds. The ‘scale of synchrony’ is often used as such threshold, indicating the distance at which the exponential function decays at e−1 ≈ 0.37 of its maximum value (e.g. Myers et al. 1997; Jarillo et al. 2018). The scale of synchrony calculated from the overall fluvial synchrogram is comparable to what has been derived from other studies of stream fishes. Over the Euclidean dimensions, the scale of synchrony was 74 km, whereas for watercourse and flow‐connected distance it was 161 and 154 km respectively. For comparison, Myers et al. (1997) reported values of c. 20 km (Euclidean) for species of pike (Esocidae) and perch (Percidae) and up to 180 km for Pacific salmon (Salmonidae) in streams. Tedesco et al. (2004) reported high synchrony values (r = 0.8) up to c. 200 km (Euclidean) for some tropical fish species.

Synchrograms can also help protect the stability of the basin level portfolio by identifying network branches that contribute the most to overall asynchrony (e.g. with shorter scales of synchrony). For instance we show how low‐order segments can contribute to the maintenance of biodiversity in river networks by promoting asynchrony among populations. The importance of headwater ecosystems as hotspots of biodiversity is well recognised (e.g. Richardson 2019) and our results thus provide further theoretical and empirical support for their protection. However, the key role of headwaters and habitat heterogeneity in promoting asynchrony is threatened by modification in catchment land‐use, flow‐regimes and the establishment of non‐native species that homogenise river landscapes globally (Poff et al. 2007; Erős et al. 2020). The consequences of this are already evident, for example in the synchronous dynamics of benthic invertebrates in modified catchments (Ruhi et al. 2018). This calls for more efficient and holistic conservation efforts that explicitly include complex riverine habitats (e.g. Colvin et al. 2019).

Our results indicate that network topology and asymmetric connectivity promote fundamental patterns of synchrony in fish populations. The analytical framework can be further relevant for the study of synchrony in any complex landscape where dispersal is prevalently directional and constrained by a dendritic habitat template (e.g. Rayfield et al. 2011; Fletcher et al. 2013). Expanding the present framework to integrate interspecific synchrony represents a promising avenue for future research, given that interspecific interactions might scale‐up to generate emerging spatial patterns at the metacommunity level (Anderson & Hayes 2018). We therefore hope that our study will stimulate further investigation of the mechanisms underlying the geometry of synchrony in complex habitats.

AUTHORSHIP

All authors contributed to developing the theoretical framework of the study. SL performed the analyses and wrote the first draft of the manuscript, with substantial contribution from LC, AFF, MJF, CJ, PD and JDO as core writing team. CJ, RR and UB developed the metacommunity simulation model. All authors provided critical feedback that helped shape the research and the final draft.

Peer Review

The peer review history for this article is available at https://publons.com/publon/10.1111/ele.13699.

Supplementary Material

Click here for additional data file.