## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Metapopulation stability in branching river networks

Edited by Andrea Rinaldo, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, and approved May 18, 2018 (received for review January 4, 2018)

### This article has a Correction. Please see:

## Significance

Metapopulation stability is a critical ecological property. Although ecosystem size has been considered as a fundamental driver of metapopulation stability, current theories developed in simplified landscapes may not be appropriate for complex branching ecosystems, such as rivers. Here, we show that a scale-independent characteristic of fractal river networks, branching complexity (measured as branching probability), stabilizes watershed metapopulations. We theoretically revealed that a strong association between branching complexity and metapopulation stability is a consequence of purely probabilistic processes. Furthermore, the stabilizing effect of branching complexity was consistently observed in metapopulations of four ecologically distinct riverine fishes. Hence, branching complexity may be a ubiquitous agent of metapopulation stability in branching ecosystems. The loss of such complexity may undermine resilience of metapopulations.

## Abstract

Intraspecific population diversity (specifically, spatial asynchrony of population dynamics) is an essential component of metapopulation stability and persistence in nature. In 2D systems, theory predicts that metapopulation stability should increase with ecosystem size (or habitat network size): Larger ecosystems will harbor more diverse subpopulations with more stable aggregate dynamics. However, current theories developed in simplified landscapes may be inadequate to predict emergent properties of branching ecosystems, an overlooked but widespread habitat geometry. Here, we combine theory and analyses of a unique long-term dataset to show that a scale-invariant characteristic of fractal river networks, branching complexity (measured as branching probability), stabilizes watershed metapopulations. In riverine systems, each branch (i.e., tributary) exhibits distinctive ecological dynamics, and confluences serve as “merging” points of those branches. Hence, increased levels of branching complexity should confer a greater likelihood of integrating asynchronous dynamics over the landscape. We theoretically revealed that the stabilizing effect of branching complexity is a consequence of purely probabilistic processes in natural conditions, where within-branch synchrony exceeds among-branch synchrony. Contrary to current theories developed in 2D systems, metapopulation size (a variable closely related to ecosystem size) had vague effects on metapopulation stability. These theoretical predictions were supported by 18-y observations of fish populations across 31 watersheds: Our cross-watershed comparisons revealed consistent stabilizing effects of branching complexity on metapopulations of very different riverine fishes. A strong association between branching complexity and metapopulation stability is likely to be a pervasive feature of branching networks that strongly affects species persistence during rapid environmental changes.

Over the past few decades, there has been growing interest in understanding how the web of living organisms maintains the temporal stability of ecosystems. Much attention has been paid to diversity-stability relationships of multispecies communities (1). However, recent recognition of intraspecific population diversity, spatial and/or life-history distinctiveness among groups of individuals, has sparked discussion of how this underappreciated component of biological complexity contributes to metapopulation stability and persistence in nature (2⇓–4). The key element of population diversity effects is the spatial asynchrony of population dynamics in metapopulations, which we define here broadly as groups of subpopulations connected via dispersal (5). Integration of multiple asynchronous dynamics over space should dampen stochastic fluctuations occurring at the level of component subpopulations, underpinning the global metapopulation stability (6, 7). Given the anticipated environmental uncertainty due to anthropogenic climate change, it is more imperative than ever to understand underlying processes of emergent stability of metapopulations (8).

Metapopulation stability is thought to be largely determined by ecosystem size (or habitat network size; hereafter, we refer to it collectively as “ecosystem size”), a fundamental driver of various ecosystem properties, such as food chain length (9, 10) and ecosystem production (11). Larger ecosystems will include probabilistically more diverse subpopulations, resulting in more stable metapopulation dynamics through statistical averaging (2, 12⇓–14). In theory, this scaling relation can be expected on statistical grounds alone if all subpopulations are equally and imperfectly synchronized over the landscape and, more specifically, if synchronies for all subpopulation combinations are describable by a single correlation coefficient *ρ* (15). This scale-dependent nature of metapopulation stability has been commonly observed in 2D systems (16⇓–18), where the key theoretical assumption of homogeneous distributions of synchrony may hold true. In those simplified landscapes, two major forces, environmental similarity (i.e., Moran effect) and dispersal, may lead to the emergence of constant synchrony at landscape scales (18⇓–20).

A growing body of evidence suggests, however, that current theories devoted to simplified landscapes may not predict emergent properties of complex branching ecosystems, an overlooked but widespread habitat geometry (21⇓⇓⇓⇓–26). Geomorphological or biological processes form naturally fractal branching structure (e.g., rivers, trees) in which geometric statistical properties (e.g., branching complexity) remain constant across spatial scales (27⇓⇓–30). Thus, branching structure should represent a dimension of system property that is orthogonal to ecosystem size. If synchronization patterns are hierarchically organized by the branching structure (i.e., a scale-invariant characteristic), then metapopulation dynamics may not be describable by existing theories developed in 1D or 2D systems. Fluvial river networks, on which humans rely for many ecosystem services (22, 31), exemplify this alternative landscape (28, 32).

Here, we hypothesize that a scale-invariant characteristic of fractal river networks, branching complexity, should stabilize the watershed metapopulation dynamics. In riverine systems, an individual branch (tributary) may represent a spatial unit of synchronized ecological dynamics because subpopulations living in a same branch should experience similar environmental signals (e.g., flood disturbance) (33, 34) and/or frequent mixing of individuals via dispersal (35, 36). Meanwhile, confluences serve as “merging” points, where branch-specific dynamics aggregate into a larger, ecologically different downstream channel (14, 21). Therefore, increased levels of branching complexity (in terms of branching prevalence) should confer a greater likelihood of integrating asynchronous dynamics over the landscape.

In this study, using a simple theoretical model, we show that this “branching complexity-stability” relationship can be a consequence of purely probabilistic processes. This emergent phenomenon arises only if we consider systematic heterogeneity of population synchrony (i.e., differential levels of synchronization between “within-branch” and “among-branch” subpopulations), the nature of branching networks that has been overlooked in previous theoretical explorations. We further evaluated our theoretical predictions using 18-y observations of fish population abundance across 31 watersheds in Hokkaido, Japan. These small watersheds are separated by the ocean and vary greatly in branching structure due to geological and climatic differences (28). Hence, this unique long-term dataset represents ideal replication of independent river networks (Fig. 1*B* and *SI Appendix*, Fig. S1). We confirmed the stabilizing effect of branching complexity on metapopulations of four fish species involving three major families of freshwater fish in the Northern Hemisphere: *Barbatula toni* (Balitoridae), *Tribolodon hakonensis* (Cyprinidae), *Oncorhynchus masou masou* (Salmonidae), and *Salvelinus leucomaenis* (Salmonidae). Our results demonstrate the fundamental role of branching structure in driving metapopulation stability, an important effect that underpins the long-term persistence of species in the face of increasing environmental variability.

## Results

### Theoretical Prediction.

We described analytical relationships between the coefficient of variation of metapopulations (CV_{m}) and branching probability *P* (Eq. **2**). Our theoretical model predicted that branching complexity should dampen temporal variability of metapopulations if within-branch population synchrony was greater than among-branch synchrony (within-branch correlation *ρ*_{wb} > among-branch correlation *ρ*_{ab}) (lines in Fig. 2). However, these relationships were reversed if *ρ*_{ab} exceeded *ρ*_{wb} (Fig. 2). Metapopulation size *N*, the number of interacting subpopulations in a network (Fig. 1*A*), did not change the relation of branching complexity and CV_{m} (Fig. 2).

Importantly, for both scenarios, stochastic simulations suggested that the uncertainty of expected CV_{m} decreased with increasing branching probability *P* (dots in Fig. 2). This stemmed from low structural variability in networks with high branching probabilities: In our model, stochastic variation in the number of subpopulations within a branch (*n*) fell nonlinearly as branching probability increased [*Materials and Methods*). This facilitated greater structural stability of randomly branching networks, leading to more predictable spatial synchrony and global stability patterns of metapopulations.

Contrary to the effects of branching complexity, effects of metapopulation size were not directly related to metapopulation stability; rather, they were contingent upon branching probability *P* (Fig. 2). As a result, vague relationships between metapopulation size and CV_{m} were observed except when *ρ*_{wb} equaled *ρ*_{ab} (dots in *SI Appendix*, Fig. S2). Specifically, metapopulation size effects were apparent only when branching probability was low (*P* < 0.3; Fig. 2 and *SI Appendix*, Fig. S2). Therefore, the effects may be hardly detectable in real river networks, where both branching probability *P* and metapopulation size *N* vary simultaneously (dots in *SI Appendix*, Fig. S2). The metapopulation stability depended exclusively on *N* and *ρ* when all subpopulations were equally correlated (*ρ*_{wb} = *ρ*_{ab}) (Fig. 2).

The above predictions were also supported in more realistic metapopulations. Even when abundances and CVs of subpopulations were unequal (*SI Appendix*), branching complexity consistently stabilized metapopulations as long as within-branch synchrony was higher than among-branch synchrony (*SI Appendix*, Figs. S3–S8). In contrast, effects of metapopulation size were unclear in most parameter settings (*SI Appendix*, Figs. S9–S14). Furthermore, an individual-based simulation model, which incorporated detailed ecological processes that underlie population synchrony (dispersal, local demography, and environmental stochasticity), reproduced the patterns expected from our analytical model (*SI Appendix*, Figs. S15–S17). Thus, the relationship between branching complexity and metapopulation stability was not sensitive to simplified theoretical assumptions.

### Empirical Evidence.

We first estimated detectability of the study species with a Bayesian multinomial model that accounted for random variation across sites and years; thus, this approach corrects not only for imperfect detection but also for observer-specific errors (*SI Appendix*). On average, detection probabilities were reasonably high across four study species (detection probability = 0.73–0.92; *SI Appendix*, Table S4). The Bayesian state-space model was fitted to the detectability-corrected data. We reconstructed complete 18-y population trends of four freshwater fish species (*SI Appendix*, Figs. S18–S21 and model check fit in Fig. S22) to quantify the stability of metapopulations (CVs of watershed-scale metapopulation dynamics).

Despite the substantial variation in the ecological traits of study species (*SI Appendix*, Table S3), a fixed-response model [no species-specific responses to explanatory variables; widely applicable information criterion (WAIC) = 345.9] outperformed a variable-response model that assumed differential species responses (WAIC = 365.1). Thus, there was no statistical evidence of species-specific responses to the explanatory variables. Furthermore, the Bayesian *P* value for the fixed response model was 0.53, indicating that the model fit was reasonably good.

In the fixed-response model, watershed metapopulation variability (CV) decreased significantly with increasing branching probability (Fig. 3 and Table 1), providing empirical evidence for the stabilizing effect of branching complexity. Moreover, as predicted by stochastic simulations, increased branching probability decreased the uncertainty of watershed metapopulation stability; that is, highly fluctuating metapopulations are less likely to occur in complex river networks (Fig. 3). Thus, branching complexity has dual important impacts on watershed metapopulation dynamics: dampening average and maximum levels of watershed metapopulation variability. Other environmental factors, including watershed area (i.e., an empirical proxy for metapopulation size *N*), had little effect on the watershed metapopulation stability (Table 1).

## Discussion

Ecosystem size, a variable intimately linked to metapopulation size, has been believed to be a general predictor of metapopulation stability. However, in branching networks, size effects were poorly supported by theory (metapopulation size *N*) or empirical analyses (watershed area; an empirical proxy for metapopulation size). Instead, a scale-invariant characteristic of fractal river networks, branching complexity, emerged as an important stabilizing agent across four riverine fishes, in line with our theoretical predictions. The species-invariant effect of branching complexity implies the existence of a general statistical property underpinning watershed metapopulation stability. Our findings are particularly valuable for biodiversity conservation in complex landscapes because species living in branching ecosystems are likely to exhibit classical metapopulation dynamics that are inherently vulnerable to rapid environmental changes (37).

Our results are in general agreement with a few previous studies suggesting that branching complexity enhances metapopulation stability or persistence (24, 25, 38, 39). However, our theoretical exploration provides two insights into metapopulation stability under natural conditions, where within-branch synchrony greatly exceeds among-branch synchrony. First, the stabilizing effect of branching complexity is an emergent phenomenon of purely probabilistic processes. Thus, our theoretical prediction should be applicable to many species inhabiting river networks and other branching ecosystems (e.g., trees, caves, mountain ranges). Second, the form of the “complexity-stability” relationships changed very little along a gradient of metapopulation size *N*. Hence, the impact of branching complexity should appear across a wide range of spatial scales. These two properties, which have not been recognized in previous studies, may explain why stabilizing effects of branching structure were invariant across four ecologically distinct fishes.

Many potential mechanisms could promote within-branch population synchrony, including external forces acting on subpopulations. A common feature of riverine systems is that environmental signals (e.g., flood disturbance) cascade down along branches of stream networks (33). Environmental conditions are highly distinctive among branches due, in part, to variability of local geological and geomorphological features, such as slope, aspect, and soil porosity (14, 40). Such spatial patterning of environmental effects may cause subpopulation responses that are unique to each branch. Alternatively, systematic heterogeneity of population synchrony can stem from internal properties of organisms, such as dispersal. Within-branch mixing of individuals is typically more frequent than among-branch dispersal (35, 36, 41), because cross-branch movements impose longer travel distance along a network (22) or extra costs for overland dispersal [e.g., salamanders (41)] to potential dispersers. This structural constraint, per se, should enhance coherent subpopulation dynamics within a branch. In addition, the limited gene flow across branches, combined with branch-specific environments, may foster local adaptation. Life history variation among streams (e.g., spawning timing) has been found in various aquatic organisms (40, 42, 43) and is thought to bring differential responses to shared environmental fluctuations (2). All of the above mechanisms can act in concert and form the basis of the observed stabilizing effect of branching complexity.

Our theoretical model also suggested that the stabilizing effect of branching complexity should fade out, or even reverse, as among-branch population synchrony (*ρ*_{ab}) increases (Fig. 2). This pattern may arise from homogenization of branch-specific environmental conditions, as well as strong local demographic/environmental stochasticity. Human alterations, such as flow regulations by dams (44, 45), land use change (11, 31), and artificial propagation programs of fishes (3), have homogenized branch distinctiveness of flow, physical conditions, and resource type and availability across the globe. Thus, stabilizing effects of branching complexity may be diminished in severely altered landscapes. Such human-induced processes could be responsible for the apparent lack of stabilizing effects of branching complexity in one notable study, which assessed salmon metapopulation stability across pristine and human-altered watersheds (46). In this sense, our study system offered an excellent opportunity to quantify branching effects since the study watersheds had very little human impact (*SI Appendix* and *SI Appendix*, Table S1). Thus, any influence of unmodeled confounding factors should be minimal. Furthermore, our advanced statistical approach fostered by exhaustive sampling efforts (i.e., two-pass removal across all sites for 18 y) allowed us to account for unavoidable observation errors (e.g., imperfect detection), a critical issue of long-term monitoring programs. These strengths may have helped to illuminate the stabilizing effect of branching complexity that otherwise can be easily missed.

In our empirical analysis, direct comparisons between within-branch and among-branch population synchronies were impossible because they require a greater number of within-network replication than available even in our large dataset. Therefore, the possibility exists that other ecological processes could be responsible for our empirical results. Nevertheless, we are unaware of any general theories that could provide convincing explanations for the observed patterns. Furthermore, it is difficult to envision that detailed ecological traits (e.g., dispersal rate), which are often species-specific, explain the consistent stabilizing effect of branching complexity on very different fish species.

Our proposed theory arises from the physical architecture of spatial structure and is applicable across a wide range of spatial scales. In addition, the wealth of indirect evidence implies the existence of systematic heterogeneity of population synchrony (discussed above). At present, our theory should best explain the observed impacts of branching complexity.

An intriguing avenue for future research is to ask how our findings can be scaled up to metacommunities (47). Branching structure is a primary determinant of species diversity patterns in experimental (48⇓–50) and natural branching landscapes (51, 52); however, it remains elusive how long-term dynamics of metacommunities are mediated by spatial network structure (but ref. 50). We expect that interspecific synchrony may play an important role in determining the relationship between branching complexity and metacommunity dynamics. Although the present study focused on single-species metapopulations separately, an extension of our framework that incorporates interspecific synchrony terms (53) may offer a promising tool to predict metacommunity-level stability in branching ecosystems.

Branching systems are ubiquitous in nature. Biological, geomorphological, and other habitat-forming processes can create fractal branching networks. In those systems, the scale-invariant complexity of ecosystems can drive metapopulation stability via purely probabilistic processes. Hence, there is no reason to believe that our fundamental finding is any less relevant for other aquatic and terrestrial organisms living in branching ecosystems. Recognition of the role of branching structure in driving metapopulation stability should increase our ability to manage and conserve species in complex landscapes.

## Materials and Methods

### Theory.

Branching networks can be depicted by nodes connected via edges (25, 39, 54). In our framework, following Yeakel et al. (39), nodes denote populations living in discretized river sections with a scale length *l*, which include either nonbranching or branching river sections (Fig. 1*A*). The scale length *l* defines the spatial scale of local reproductive interactions (hereafter “subpopulation”) and may depend on attributes of the species of interest, such as dispersal ability. We assumed binary river networks, in which rivers bifurcate at confluences.

We define metapopulation stability as relative fluctuations of metapopulation abundance (or density) through time. The major goal of our theoretical exploration is to describe metapopulation stability by the following properties.

#### Metapopulation size *N* (subpopulation *i* = 1, 2, …, *N*).

Metapopulation size defines the number of interacting subpopulations and is closely related to ecosystem size *L* (total river length of the watershed). The expected metapopulation size can be expressed as *N* = *L*/*l*. For all subpopulations, we assumed an equal temporal mean *μ* and SD *σ* of subpopulation abundance, and thus an identical coefficient of variation [CV of component populations (CV_{p}) = *σ*/*μ*]. Various ecological processes (e.g., environmental fluctuations) impact temporal dynamics of subpopulations *x*(*t*). Network structure influences spatial patterning of these factors (discussed below).

#### Branching probability *P*.

Each node is assigned to be either a branching (or upstream terminal) node with probability *P* or a nonbranching node with probability 1 − *P*. In this formulation, an individual branch represents a series of nonbranching nodes terminated at a branching (or terminal) node (equivalent to “link” in geomorphological terms; Fig. 1*A*). The number of subpopulations involved in branch *j* (*n*_{j}) is a realization of a random variable drawn from a geometric distribution with success (branching) probability *P* [*NP* and *sensu lato*; details are provided in *SI Appendix*). Note that a geometric distribution is a discrete version of an exponential distribution, which has been traditionally used to describe random branch (link) lengths (e.g., ref. 55). The above procedure is equivalent to sampling a river network of metapopulation size *N* from a pool of random river networks with branching probability *P*.

#### Population synchrony *ρ*.

We expressed the degree of population synchrony using the Pearson’s correlation coefficient. Specifically, we considered two levels of synchronization: within-branch (*ρ*_{wb}) and among-branch (*ρ*_{ab}) correlations. Any population (node) combination has a correlation coefficient of either *ρ*_{wb} or *ρ*_{ab} depending on whether the two are in the same branch (Fig. 1*A*). These correlation coefficients were used to characterize temporal subpopulation dynamics *x*(*t*). Population synchrony may emerge as a consequence of ecological processes, including environmental similarity and/or dispersal of individuals, among others. Thus, our model includes implicitly various synchronization forces and spatial interactions.

The global metapopulation stability is a function of component dynamics. Temporal trajectory of metapopulation abundance *X*(*t*) (*X* = *x*_{1} *+ x*_{2} +…+ *x*_{N}) will have a variance of

where *n*_{j} ≥ 2. The parameter φ represents the total number of within-branch combinations of nodes, and its expected value is *SI Appendix*). Eq. **1** shows that the temporal variance of a metapopulation is dependent on the summed variances of each component subpopulation and on the summed covariances among all possible combinations of these subpopulations. By taking the square root and dividing both sides of the equation by *Nμ* (i.e., expected metapopulation abundance), the temporal metapopulation stability (or instability), expressed as relative fluctuations of metapopulation abundance (CV_{m}), can be written as

where CV_{p} is the CV of component subpopulations (*σ*/*μ*) and _{p} does not change the functional relationship of CV_{m} and *P* or *N* because Eq. **2** has the form of a power law. When all populations are correlated equally (*ρ*_{wb} = *ρ*_{ab} = *ρ*), Eq. **2** can be further reduced to

Under this scenario, CV_{m} is independent of branching probability *P*. Eq. **3** is identical to the equation proposed by Doak et al. (15), who formulated the CV for a community aggregate composed of *N* species with equivalent CVs.

We also conducted stochastic simulations to illustrate uncertainty of CV_{m}. We considered 3 × 3 combinations of parameters (*ρ*_{wb} = 0.1, 0.5, 0.9; *ρ*_{ab} = 0.1, 0.5, 0.9). This range of correlation coefficients is comparable to values often observed in nature (e.g., ref. 56). Under each scenario, we replicated 1,000 river networks with varying branching probability *P* and metapopulation size *N*, both of which were drawn randomly from uniform distributions (range: *P* = 0.1–0.9, *N* = 50–150). Metapopulation size *N* was rounded before performing simulations. The number of subpopulations in branch *j* was sampled from a geometric distribution [*n*_{j} ∼ Ge(*P*)] until the sum of *n*_{j} equals *N* without truncation of any nodes in a network. CV_{m} was calculated using Eq. **1** [i.e.,

Finally, we developed supplementary simulation models to confirm robustness of our theoretical predictions in realistic ecological settings. Specifically, we examined effects of the following factors: (*i*) unequal subpopulation abundances and CVs and (*ii*) detailed ecological processes (e.g., dispersal) that underlie population synchrony (*SI Appendix*). All procedures were performed in R 3.3.1 (57).

### Empirical Test.

#### Estimation of long-term trends of stream fish metapopulations.

We assembled time series data of fish abundance across 31 protected watersheds (27 entire watersheds and four subwatersheds) of Hokkaido, Japan. Long-term fish monitoring began in 1971, but effective sampling methods were introduced in 1999 (a combination of cast net sampling and electrofishing). In each watershed, fish abundance was surveyed at two to five permanent sampling sites using a two-pass removal method (*SI Appendix*, Fig. S1), resulting in a total of 106 sites across watersheds (mean survey area per site: 172.9 ± 115.7 m^{2}). Fish abundance data were collected from July to September (typically July and August) with irregular interannual intervals (0- to 3-y intervals in most cases; data availability is provided in *SI Appendix*, Table S2). We quantified metapopulation stability of four fish species (*B. toni*, *T. hakonensis*, *O. m. masou*, and *S. leucomaenis*; their ecological traits are provided in *SI Appendix*, Table S3) during a period from 1999–2016. After filtering the data to remove less reliable information, the number of watersheds analyzed varied from 16 to 31 depending on species. Further details of study locations and data selection are provided in *SI Appendix*.

In our dataset, metapopulation stability was not directly comparable among the watersheds because temporal data availability varied as described above (*SI Appendix*, Table S2). To overcome this difficulty, we employed a state-space modeling approach within a Bayesian framework. Bayesian state-space models can infer complete population trends from sparse time series data while accounting for observation errors (58). Therefore, this approach is one of the best options currently available for analyses of empirical population time series (16). The following Bayesian state-space model was applied to each species separately.

##### Data model.

In our model, detectability-corrected fish abundance *N*_{ti(j)} (year *t*, site *j* nested within watershed *i*; a discussion of detectability correction is provided in *SI Appendix*) was assumed to follow a Poisson distribution with an expected mean of *λ*_{ti(j)}:

where AREA_{ti(j)} denotes areal sampling area (square meters) (i.e., offset term). The parameter log *d*_{obs,ti(j)} represents the log-transformed fish density (individual 100 m^{−2}) that involves site-specific year effects (i.e., observation error). The log-transformed fish density was drawn from a normal distribution:

The inclusion of parameter *σ*_{obs} allows us to account for random variations of site-specific year effects that may stem from ecological and/or artificial factors. For example, *σ*_{obs} may account for spatiotemporal variability imposed by local biotic (e.g., intraspecific competition, immigration) and abiotic processes. Thus, the mean parameter log *d*_{ti(j)} represents “unbiased” fish density at the site scale (hereafter, “local density”).

##### Process model.

The local density log *d*_{ti(j)} was then described as the deviation from the temporal metapopulation trajectory:

where the parameter *σ*_{α} governs the degree of deviation from the metapopulation-level density in year *t* and watershed *i* (log *D*_{ti}). This model assumption is robust and conservative, as no single “litmus test” exists to determine whether the local aggregations are independent subpopulations or parts of a larger panmictic population in a stream continuum (59, 60).

The temporal trajectory of metapopulation-level density log *D*_{ti} [i.e., an empirical proxy for theoretical metapopulation trajectory *X*(*t*)], was modeled as follows:

where log *r*_{ti} is the log-transformed metapopulation growth rate in year *t* and watershed *i*. We assumed density-independent metapopulation growth rates because intraspecific competition is unlikely to operate at a watershed scale. This is a common implicit assumption made in previous empirical and theoretical studies (2, 8, 61). The log-transformed metapopulation growth was assumed to follow a multivariate normal distribution:

where *ζ* is the vector of watershed-specific temporal mean of metapopulation growth rate (i.e., log *r*_{mean,1}, log *r*_{mean,2}, …) and Ʃ_{r} represents the variance-covariance matrix of log *r*_{ti}. The advantage of this hierarchical structure is that parameter estimates of each watershed can be improved over those that would be obtained by fitting separate watershed-specific models (62). This is apparent in our dataset, because data-poor watersheds can borrow information from data-rich watersheds.

Vague priors were assigned to the parameters [i.e., normal distributions for log *r*_{mean} (mean = 0, variance = 10^{3}), an inverse-Wishart distribution for Ʃ_{r} (df = *n* + 1, where *n* is the number of watersheds analyzed), truncated normal distributions for *σ*_{obs} and *σ*_{α} (mean = 0, variance = 10^{3}, range = 0–30)]. The model was fitted to the data using JAGS (version 4.1.0) and the package “runjags” (63) in R 3.3.1 (57). Three Markov chain Monte Carlo (MCMC) chains were run with 75,000 iterations (25,000 burn-in), and 500 samples per chain were used to calculate posterior probabilities. Convergence was assessed by examining whether the R-hat indicator of each parameter approached a value of 1 (64).

Using median estimates, we quantified the metapopulation stability. Specifically, we calculated CVs of watershed metapopulation trajectory *D*_{ti}. Detrending was performed before the calculation to avoid biased estimates of CVs (61): The SD of residuals from a fitted lowess smoother (function “lowess,” a smoothing span of 2/3 of the data) was divided by the original (i.e., before detrending) temporal mean of *D*_{ti}.

#### Watershed characteristics.

We estimated watershed characteristics (branching probability and watershed area) and climatic variables (annual cumulative precipitation and mean air temperature) using QGIS version 2.12.3 (available at https://qgis.org/en/site/). Branching probability was estimated by fitting exponential distributions to histograms of branch length (kilometers, river segment between successive confluences or a confluence and the outlet/upstream terminal) for each watershed as branch length ∼ Exp(*υ*) (55). An exponential distribution is a continuous version of a geometric distribution, and a derived statistical quantity 1 − *e*^{−υ} corresponds to the theoretical branching probability *P*. We used river polylines (defined as >1 km^{2} of watershed area) of the entire watershed for the estimation of branching probability. Watershed area was used as a proxy for metapopulation size and calculated as the entire drainage area upstream of the ocean outlet (27 watersheds) or confluence with another major tributary (four subwatersheds, if all sampling sites were located in a subwatershed). All watershed characteristics were calculated based on a digital elevation map of 90-m resolution (available from the US Geological Survey; https://lta.cr.usgs.gov). Climatic variables were obtained from the Japan Meteorological Agency (available at nlftp.mlit.go.jp/ksj/index.html). We used 10-km mesh data of annual cumulative precipitation and mean air temperature (averaged for 1981–2010). We averaged values of climate meshes for each watershed. Finally, to examine the influence of fish stocking on metapopulation stability, we obtained a record of cumulative number of hatchery fish [masu salmon (*O. m. masou*)] released into each watershed during a period from 1955–2016.

Elevation, watershed land cover, and weir density were not used in the following regression analysis because these variables were less informative and/or risked influences of multicollinearity. The detailed variable selection procedure is described in *SI Appendix*.

#### Stability analysis.

We fitted hierarchical log-normal models to estimated CVs using JAGS (version 4.1.0) to investigate impacts of watershed characteristics on metapopulation stability. Log-transformed CVs met the normality assumption (Shapiro–Wilk normality test for residuals, *P* > 0.58) and had no statistically evident outliers (Grubbs’ test, *P* = 0.42). The detrended CV_{s(i)} (watershed *i* for species *s*) was modeled as follows:

where *BP*_{i}, *WA*_{i}, *Prec*_{i}, *Temp*_{i}, and *Hatchery*_{i} denote branching probability, watershed area (square kilometers), cumulative precipitation (millimeters), air temperature (°C), and the cumulative number of hatchery fish released, respectively. *Nsite*_{i} is the number of sampling sites in watershed *i* and was included as a control variable. Since *Hatchery*_{i} was positively correlated with *WA** _{i}* (Pearson’s

*r*= 0.72), we used residuals of a linear relationship fitted between

*Hatchery*

_{i}and

*WA*

_{i}to avoid multicollinearity. The residuals provide a relative measure of

*Hatchery*

_{i}, independent of

*WA*

_{i}(e.g., positive residuals indicate greater than expected values from a linear regression). All explanatory variables were standardized before the analysis (mean = 0, SD = 1).

We constructed two competing models: fixed-response (random intercept) and variable-response models (random intercept and slope). In a fixed-response model, only the intercept *β*_{0,s} varies randomly by species [*β*_{0,s} ∼ Normal(*β*_{0,global}, *σ*_{β0}^{2})] while keeping slopes constant among species (*β*_{k,constant} = *β*_{k,s} for *k* = 1–6). On the other hand, in a variable-response model, both the intercept and slope vary randomly among species [*β*_{k,s} ∼ Normal(*β*_{k,global}, *σ*_{βk}^{2}) for *k* = 0–6]. To compare the two competing models, we calculated the WAIC (65) using the R package “loo.” The WAIC provides, like other information criteria, a measure of model fit that is penalized by model complexity, but it has a sound theoretical foundation in Bayesian statistics and an applicability to complex hierarchical models like ours. We also checked the model performance with a Bayesian *P* value. The Bayesian *P* value, which was based on a sums-of-square discrepancy, approaches 0.5 if the model perfectly reproduces the data (66).

Vague priors were assigned to the parameters [i.e., normal distributions for *β*_{k,constant} and *β*_{k,global} (mean = 0, variance = 10^{4})] and truncated normal distributions for *σ*_{cv} and *σ*_{βk} (mean = 0, variance = 10^{4}, range = 0–100). Three MCMC chains were run with 7,500 iterations (2,500 burn-in), and 500 samples per chain were used to calculate posterior probabilities. Convergence was assessed as described above.

#### Data accessibility.

The codes (R and JAGS) and data reported in this study are available in *SI Appendix* and Dataset S1.

## Acknowledgments

We thank the numerous people who worked hard to collect and organize the large datasets of protected watersheds. We are grateful to three anonymous reviewers and to Shota Nishijima, Keisuke Atsumi, and Shin-ichiro S. Matsuzaki for their constructive comments on the manuscript. We thank Jesús N. Pinto-Ledezma, Josep Padullés Cubino, and the theory group at the University of Minnesota for helpful discussion. This research was partly supported by Japan Society for the Promotion of Science KAKENHI Grant 18K06404 (to A.T.).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: hanabi0111{at}gmail.com.

Author contributions: A.T. designed research; A.T., N.I., H.U., S.O., J.C.F., and F.N. performed research; A.T. contributed new reagents/analytic tools; A.T. and S.O. analyzed data; and A.T., N.I., H.U., J.C.F., and F.N. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1800060115/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- ↵
- ↵
- Moore JW,
- McClure M,
- Rogers LA,
- Schindler DE

- ↵
- ↵
- ↵
- ↵
- Dey S,
- Joshi A

- ↵
- Anderson SC,
- Moore JW,
- McClure MM,
- Dulvy NK,
- Cooper AB

- ↵
- Sabo JL,
- Finlay JC,
- Kennedy T,
- Post DM

- ↵
- ↵
- Finlay JC

- ↵
- ↵
- ↵
- Moore JW, et al.

- ↵
- ↵
- Matsuzaki SS,
- Kadoya T

- ↵
- Shelton AO, et al.

- ↵
- ↵
- ↵
- Pomara LY,
- Zuckerberg B

- ↵
- ↵
- ↵
- Mari L,
- Casagrandi R,
- Bertuzzo E,
- Rinaldo A,
- Gatto M

- ↵
- Grant EHC

- ↵
- ↵
- Tonkin JD, et al.

- ↵
- Rinaldo A,
- Rigon R,
- Banavar JR,
- Maritan A,
- Rodriguez-Iturbe I

- ↵
- Rodríguez-Iturbe I,
- Rinaldo A

- ↵
- Mandelbrot BB,
- Pignoni R

- ↵
- Veitzer SA,
- Gupta VK

- ↵
- Allan JD,
- Castillo MM

- ↵
- Rinaldo A,
- Gatto M,
- Rodriguez-Iturbe I

- ↵
- Nakamura F,
- Swanson FJ,
- Wondzell SM

- ↵
- ↵
- ↵
- Junker J, et al.

*Cottus gobio*). Conserv Genet 13:545–556. - ↵
- Fronhofer EA,
- Altermatt F

- ↵
- ↵
- Yeakel JD,
- Moore JW,
- Guimarães PR Jr,
- de Aguiar MAM

- ↵
- ↵
- Campbell Grant EH,
- Nichols JD,
- Lowe WH,
- Fagan WF

- ↵
- ↵
- ↵
- Poff NL,
- Olden JD,
- Merritt DM,
- Pepin DM

- ↵
- ↵
- ↵
- ↵
- Carrara F,
- Altermatt F,
- Rodriguez-Iturbe I,
- Rinaldo A

- ↵
- ↵
- ↵
- ↵
- Tonkin JD, et al.

- ↵
- Delsol R,
- Loreau M,
- Haegeman B

- ↵
- Eros T,
- Schmera D,
- Schick RS

- ↵
- Peckham SD,
- Gupta VK

- ↵
- ↵
- R Core Team

- ↵
- Kéry M,
- Schaub M

- ↵
- Quinn TP,
- Rich HB,
- Gosse D,
- Schtickzelle N

*Oncorhynchus nerka*) population structure in Alaska, USA. Can J Fish Aquat Sci 69:297–306. - ↵
- Schtickzelle N,
- Quinn TP

- ↵
- Anderson SC,
- Cooper AB,
- Dulvy NK

- ↵
- Lunn D,
- Jackson C,
- Best N,
- Thomas A,
- Spiegelhalter D

- ↵
- ↵
- Gelman A,
- Hill J

- ↵
- Watanabe S

- ↵
- Kéry M

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Biological Sciences
- Ecology

## Jump to section

## You May Also be Interested in

*Ikaria wariootia*represents one of the oldest organisms with anterior and posterior differentiation.