Weakened portfolio effects constrain management effectiveness for population aggregates.

Population diversity can reduce temporal variability in aggregate population abundances in a process known as the portfolio effect. Portfolio effects may weaken, however, due to greater synchrony among component populations. While weakened portfolio effects have been previously documented, the consequences of reduced stability on meeting conservation goals for population aggregates that are harvested (e.g., stock aggregates in fisheries) are rarely quantified. Here, we demonstrate how changes in variability within components, synchrony among components, and population productivity interact to influence the probability of achieving an array of management objectives for Fraser River sockeye salmon: a stock aggregate of high economic, ecological, and cultural value. We first present evidence that component variability and synchrony have increased over the last two decades, consistent with a weakening portfolio effect. We then parameterize a stochastic, closed-loop model that simulates the population dynamics of each stock, the fishery that harvests the stock aggregate, and the management framework used to establish mixed-stock exploitation rates. We find that while median aggregate abundance and catch through time were relatively insensitive to greater aggregate variability, catch stability and performance metrics associated with achieving management targets generally declined as component variability and synchrony increased. A notable exception we observed is that harvest control means that scale exploitation rates based on aggregate abundance may be more effective as synchrony increases. Reductions in productivity led to broad declines in performance, but also moderated the impacts of component variability and synchrony on the proportion of component stocks above management targets and catch stability. Our results suggest that even precautionary management strategies that account for declines in productivity may underestimate risk, particularly to socioeconomic objectives, if they fail to consider changes in aggregate variability. Adequately incorporating changes in portfolio effect strength may be particularly relevant when developing recovery strategies that are robust to climate change, which is likely to increase synchrony and component variability.

Introduction

Metapopulations consist of component populations whose dynamics vary due to life history, local environmental conditions, or stochastic processes. Asynchrony among components can reduce temporal variability of the aggregate resulting in greater productivity, increased availability of ecosystem services, and improved resilience—a process commonly referred to as the portfolio effect (Doak et al. 1998, Tilman 1999). Portfolio effects have been identified as key stabilizing processes across a range of taxa and ecosystems (Anderson et al. 2013) and have been used as justification for systems‐based approaches for management and conservation (Link 2018). A portfolio approach emphasizes monitoring and conserving ecological aggregates, in addition to a focus on individual components (Anderson et al. 2015, Link 2018).

A common way to quantify a portfolio's risk performance is via aggregate variability (e.g., the temporal coefficient of variation, CV) of multiple populations. Due to statistical averaging alone, aggregate variability within an ecological portfolio will decrease as the number of components increases (Doak et al. 1998). For example, the Bristol Bay sockeye salmon (Oncorhynchus nerka) fishery encompasses an aggregate of nine major river systems, each containing tens to hundreds of largely independent spawning populations (Schindler et al. 2010). This population diversity reduces aggregate variability in spawner returns and catches, as well as the probability of fishery closures, relative to a hypothetical fishery containing fewer populations (Hilborn et al. 2003; Schindler et al. 2010). Thus, maintaining population diversity is an effective way to promote resilience, particularly when future environmental conditions are uncertain (Anderson et al. 2015).

Yet even when the number of components is stable, ecological portfolios can exhibit changes in aggregate variability that compromise their ability to provide ecosystem services. For example, although the number of stocks (one or more populations managed as a discrete unit) within California's Central Valley Chinook salmon (O. tshawytscha) aggregate has not declined, aggregate returns to the fishery have collapsed and become increasingly variable in recent decades (Carlson and Satterthwaite 2011, Satterthwaite and Carlson 2015). In this case, decreased stability was associated with greater variability in the returns of individual stocks and greater synchrony among stocks (Satterthwaite and Carlson 2015), and the resulting impacts may have been exacerbated by declines in survival (Lindley et al. 2009). While the stock aggregate still exhibits a portfolio effect (i.e., aggregate variability is reduced relative to that of individual stocks), the buffering conferred by its diversity is substantially weaker than it was historically (Carlson and Satterthwaite 2011). Ultimately, increased aggregate variability has resulted in substantial socio‐economic costs to the region, with the probability of falling below a threshold resulting in severe fishing restrictions or closures increasing by more than 10‐fold relative to historic levels (Yamane et al. 2018).

Aggregate variability is driven by the variance and covariance of component populations. Thus, it can be decomposed into two metrics. the weighted mean coefficient of variation among components (referred to here as “component variability”) and an index of synchrony (Loreau and de Mazancourt 2008, Thibaut and Connolly 2013). While increases in either component variability or synchrony will decrease an aggregate's stability and weaken its portfolio effect, each process can produce unique challenges to systems‐based management. For instance, the dynamics of individual populations become less predictable as component variability increases, raising the risk of overharvest of individual stocks. As synchrony increases, otherwise localized boom‐and‐bust cycles become more widespread and harvesters are less able to shift effort among component stocks, leading to reduced profits (Cline et al. 2017). In contrast, when synchrony is relatively low, divergent dynamics among stocks will reduce the impact of these changes at the aggregate level. Changes in aggregate variability within Central Valley Chinook salmon appear to be linked to increased hatchery contributions (Satterthwaite and Carlson 2015, Yamane et al. 2018); however, weakening portfolio effects may also be associated with anthropogenic disturbance more broadly (Griffiths et al. 2014) or large‐scale environmental processes (Kilduff et al. 2015).

Negative ecological and socioeconomic outcomes associated with relatively high aggregate variability may be exacerbated by changes in other characteristics of population dynamics, particularly declines in productivity. In many fisheries, reduced abundance and lower fishery yields are associated with declines in per capita productivity (Peterman and Dorner 2012, Minto et al. 2014, Britten et al. 2016), due to, for example, shifts in survival or fecundity. Furthermore, the negative impacts of productivity declines on managed systems are likely to increase due to persistent stressors such as climate change (Oliver et al. 2015). While changes in component variability, synchrony, or underlying productivity may independently impact exploited systems, it is unclear how these distinct processes may interact with one another to influence management performance.

In this study, we explore how increased aggregate variability, via changes in either component variability or synchrony, influences the probability of achieving a range of conservation and management objectives using Fraser River sockeye salmon as a case study. Fraser River sockeye salmon are particularly well suited to an evaluation of management actions associated with portfolio effects because of the large number of component stocks that contribute to multiple mixed‐stock fisheries and the high economic, ecological, and cultural value of the aggregate. We first present a retrospective analysis that reveals aggregate temporal variability within the Fraser River aggregate has increased in recent years due to greater component variability, as well as greater synchrony. We then use stochastic, closed‐loop simulations to evaluate how changes in component variability and synchrony influence the probability of meeting a suite of conservation‐ and catch‐based performance metrics under the current management regime. Finally, we repeat the simulations under a range of productivity scenarios to identify how component variability and synchrony interact with changes in productivity to shape dynamics. The closed‐loop simulation framework is analogous to management strategy evaluations (MSE), where the performance of multiple management procedures is evaluated collaboratively by managers, analysts, and stakeholders (Punt et al. 2016). Like an MSE, this multi‐step approach allows us to assess changes in management performance along a gradient of scenarios, in this case, from historically observed levels of low aggregate variability and high productivity to heavily synchronized, unproductive dynamics consistent with degraded systems.

Methods

Sockeye salmon biology, fisheries, and data sources

Sockeye salmon is an anadromous, semelparous fish distributed throughout the northern Pacific. Populations in southern British Columbia typically rear as juveniles in freshwater lakes for 1–2 yr, mature in the Gulf of Alaska, and return to spawn as 2–5 yr olds (Burgner 1991). Under Canada's Wild Salmon Policy (WSP), Pacific salmon status is assessed at the scale of conservation units (CUs): “a group of wild salmon sufficiently isolated from other groups that, if lost, is unlikely to recolonize naturally within an acceptable timeframe” (DFO 2005). During WSP assessments, a range of criteria is used to assign CUs to green, amber, and red zones, which represent increasing degrees of conservation concern and management intervention (DFO 2005). Although CUs are not precisely equivalent to stocks within DFO's management framework, we use the terms interchangeably here. We also describe performance on management units, (MUs), which consist of multiple CUs with a common return migration phenology.

The Fraser River aggregate is Canada's largest sockeye salmon run, but is increasingly vulnerable to a range of threats including anthropogenic development, overexploitation, and climate change (COSEWIC 2017). Productivity declined in the 1990s, resulting in frequent fishery closures and an emergency federal inquiry (Cohen 2012). While there have been signs of recovery in recent years for specific CUs, recruitment continues to be highly variable with numerous CUs at very low productivity levels (Grant et al. 2018) and 9 of 19 assessed in the red zone (DFO 2018). Although status is assessed at the CU level, Fraser River sockeye salmon fisheries are predominantly managed at the scale of MUs and the majority of catch occurs as spawners migrate through nearshore marine areas. Thus, shifting fishery openings to coincide with migration phenology can be used to constrain effort at the MU, but not the CU, scale (though overlap among MUs exists). The Fraser River sockeye salmon CUs included in this analysis are listed in Table 1.

Table 1

Sockeye salmon management units (MU) and component conservation units (CU) within the Fraser River aggregate

MU and CU	Years	Stock–recruit model
Early Stuart
Takla‐Trembleur	1948–2011	Larkin
Early summer
Bowron	1948–2011	Ricker
Shuswap‐ES	1980–2011	Larkin
North Barriere	1967–2011	Ricker
Anderson‐Setona	1968–2011	Ricker
Nadina‐Francoisa	1973–2011	Ricker
Pitt	1948–2011	Ricker
Summer
Takla‐Trembleur‐Stuart	1948–2011	Larkin
Francois‐Fraser	1948–2011	Ricker
Kamloops‐ES	1948–2011	Ricker
Quesnel	1948–2011	Larkin
Chilko	1948–2011	Ricker
Harrison (river‐type)	1948–2011	Ricker
Late summer
Shuswap Complex‐L	1948–2011	Larkin
Lillooet‐Harrison	1948–2011	Ricker
Cultusb	1948–2000	Ricker
Seton	1965–2011	Ricker
Harrison (upstream)a	1966–2011	Ricker

Stock–recruit models reflect whether CUs have cyclic abundance (Larkin) or not (Ricker; see Base operating model for additional details). CUs in boldface type were included in the retrospective aggregate variability analysis. While the majority of the CUs are considered “wild” stocks, a subset have been enhanced (see footnotes for details).

Spawning channels have been used to enhance the productivity of these CUs. Each time series has been truncated to begin immediately post‐construction.

Since 2000, Cultus productivity has been heavily influenced by hatchery supplementation intended to reduce the risk of extirpation. Thus, although spawner abundance and catch estimates extend to 2011, stock–recruit parameters were estimated only using data collected prior to 2000.

We used annual estimates of spawner and recruit abundance for 18 CUs (Grant et al. 2011), with individual time series beginning between 1948 and 1973 (Table 1). We use the term return abundance to refer to the sum of catch and escapement (the number of returning adults that have escaped the fisheries prior to spawning plus an estimate of en route mortality) in a given return year, while recruit abundance refers to catch and escapement produced by a given brood year (see Grant et al. (2011) for details). Ten CUs with long time series were used in the retrospective analysis, while 18 CUs were included in the forward simulation model (these data are provided on Zenodo; see Data Availability).

Retrospective analysis and variability metrics

We calculated temporal changes in ln(recruits per spawner) (a metric of population productivity), aggregate return abundance, and aggregate catch in a subset of 10 CUs with time series extending back to 1948 (bold in Table 1). Since these metrics exhibit considerable interannual variability we present means calculated from four‐year (approximately one generation) moving windows. We then examined temporal changes in annual recruit abundance using three metrics of metapopulation variability adapted from Thibaut and Connolly (2013): component variability, synchrony, and aggregate variability. The first metric is measured as the mean temporal coefficient of variation of components () weighted by each component's mean recruit abundance (i.e., the sum of the components’ CV scaled by each component's mean abundance):

where is the mean recruit abundance (through time) of component CU i, is the mean abundance of the aggregate, and represents the temporal standard deviation of recruit abundance.

The second metric, the synchrony index (), reflects the relative degree of similarity in the dynamics of an ecological aggregate's components. Synchrony is analogous to comparing mean pairwise correlation coefficients, which have been used in similar analyses (e.g., Peterman and Dorner 2012), but makes no distributional assumptions, is normalized (i.e., ranges between 0 and 1 regardless of the number of components), and explicitly accounts for unequal variances among components (Thibaut and Connolly 2013). Synchrony is defined as the total temporal variance of the components (i.e., sum of all elements of the variance–covariance matrix ), divided by the variance of a hypothetical aggregate with the same component variances, but perfect covariance (Loreau and de Mazancourt 2008, Thibaut and Connolly 2013)

where represents the correlation between CUs i and j.

Finally, for the third metric, we calculated the coefficient of variation for the aggregate (CV). Thibaut and Connolly (2013) show that aggregate variability can be expressed a function of component variability and synchrony as

Thus, as synchrony increases, the variability of the aggregate converges on component variability, while aggregate variability is dampened when components vary asynchronously.

To explore historical changes in aggregate variability of Fraser River sockeye salmon, we calculated time series of synchrony, component variability, and aggregate variability using 12‐yr (three‐generation) moving windows of recruit abundance (details in Appendix S1). We used estimates of recruit, rather than spawner, abundance to account for large changes in exploitation rate over the past 70 yr. Since Fraser River CUs vary in the length of their spawner–recruit time series, we generated trends in these metrics using a subset of 10 CUs with data extending back to the 1948 brood year; however, a supplementary analysis analyzing 16 time series beginning in 1975 exhibited similar trends (Appendix S2).

Forward simulation

We used a stochastic, closed‐loop simulation to model the biological dynamics and management system of Fraser River sockeye salmon to explore how interactions between component variability, synchrony, and productivity may influence conservation and fisheries outcomes. The simulation model contains two main components: an operating model that represents a distinct ecological hypothesis about the system's biological components and a management procedure that represents harvest. The closed‐loop simulation analysis is analogous to the quantitative modeling component of an MSE, in which the performance of multiple management procedures (typically in the context of fisheries) is evaluated while accounting for many dimensions of uncertainty (Punt et al. 2016). Here, however, the analysis focuses on a comparison of multiple operating models rather than management procedures and we did not incorporate the collaborative process between managers and stakeholders that is necessary to prioritize objectives in an MSE.

The data and model code necessary to recreate the analyses are provided online (see Data Availability). The closed‐loop simulation model is the recoverySim function within the samSim package, while scripts to run various operating models and generate figures are saved separately.

Operating model

Base operating model

The dynamics of each CU were simulated using an age‐structured Larkin spawner–recruit model (Larkin 1971)

where i represents a CU, y is a given year, R the number of recruits, and S the number of spawners. The parameter represents the number of recruits produced per spawner at low abundance, while is the density‐dependent parameter that represents the reciprocal of the number of spawners that maximizes recruitment. The parameters , and represent delayed density‐dependent effects that occur in cyclic CUs, which produce highly abundant returns every 4 yr (Table 1; DFO 2018). We simulated the dynamics of non‐cyclic CUs with a Ricker model (Ricker 1975), which is equivalent to a Larkin model with delayed‐density dependence parameters set to 0.

Larkin or Ricker models are commonly rearranged so that the error term is additive and normally distributed (Peterman 1981):

where represents the standard deviation of CU‐specific recruitment deviations.

To parameterize each CU's spawner–recruit relationship we used median parameter sets (i.e., median and corresponding and estimates) generated from CU‐specific Bayesian spawner–recruit models (details in Appendix S1). To account for temporal autocorrelation in recruitment and incorporate covariation in recruitment deviations among CUs we simulated year‐ and stock‐specific deviations from the stock–recruit relationship as

where represents the previous year's recruitment deviation, represents an autoregressive lag‐1 yr (AR1) autocorrelation coefficient, and represents multivariate normally distributed errors with mean 0 and standard deviation defined by the variance–covariance matrix V, dimensioned by the number of CUs, n. We adjusted estimates of variance to account for using the transformation, , where is the variance derived from a model without autocorrelation and is the adjusted value (Pestal et al. 2011). represents the pairwise correlation coefficient between CU‐specific time series of recruitment residuals and we made the simplifying assumption that it did not vary among pairs of CUs. We assigned a value of 0.2, consistent with evidence of weak autocorrelation in the residuals of stock–recruit models from preliminary analyses (sensitivity analyses also indicated most performance metrics were not strongly impacted by this parameter; Appendix S2). Note that autocorrelated recruitment deviations were only present in CUs that were forward simulated with Ricker, not Larkin, stock–recruit models (Table 1) because Larkin models already incorporate temporal autocorrelation via delayed‐density dependent effects.

We simulated two sources of mortality in addition to the natural mortality implicitly incorporated in the spawner–recruitment model. The first represented harvest in a mixed‐CU fishery (details in Management procedure). The second represented CU‐specific en route mortality that occurs after fish enter freshwater, but prior to spawning due to a combination of natural causes (thermal stress, pathogen infection, predation, migration impediments), unreported harvest, post‐release fishing mortality, and error associated with estimates of abundance in‐river and on spawning grounds (Grant et al. 2011). We introduced additional stochasticity into the model via interannual variation in en route mortality, age at maturity, and deviations between target and realized exploitation rates (commonly referred to as outcome uncertainty or implementation uncertainty; Holt and Peterman 2006), all of which were parameterized with observations from the Fraser River aggregate (Appendix S1). In addition to the scenarios described in the next subsection, we tested the effect of alternative values that bound the ranges of observed values in a series of local sensitivity analyses (Appendix S2).

Operating model scenarios

The principal drivers of variability in aggregate abundance within the model are deviations from CU‐specific spawner–recruitment relationships (w in Eq. 6a). To explore the consequences of component variability and synchrony on management objectives, we manipulated recruitment deviations to create nine scenarios defined by unique variance–covariance matrices V (Table 2).

Table 2

Parameterization of component variability, , and synchrony, operating models, where is the standard deviation in recruitment estimate derived from CU‐specific stock recruit models (Pestal et al. 2011) and the component variability scenarios represent scalar multiples of this value

Variability	Low CVC	Reference CVC	High CVC
Low φ	0.75σR′;ρ_R=0.05	σR′;ρ_R=0.05	1.25σR′;ρ_R=0.05
Moderate φ	0.75σR′;ρ_R=0.5	σR′;ρ_R=0.5	1.25σR′;ρ_R=0.5
High φ	0.75σR′;ρ_R=0.75	σR′;ρ_R=0.75	1.25σR′;ρ_R=0.75

Values for and were intended to produce relatively large differences in component variability and synchrony, which approximated previously observed values (Fig. 2). We parameterized upper and lower scenarios for component variability as ±25% of the CU‐specific estimates of SD in recruitment , which were within ranges observed for sockeye salmon (Korman et al. 1995, Peterman et al. 2003). We parameterized lower and moderate scenarios for synchrony by adjusting the correlation coefficient to values consistent with 10‐yr moving window estimates of mean pairwise correlations in ln(R/S) among CUs during periods when productivity was weakly ( = 0.05; 1980s and 1990s) or moderately correlated ( = 0.50; 1950s, 1960s, and present). We specified a third high‐correlation scenario ( = 0.75) to represent hypothetical future increases in synchrony associated with a 50% increase in pairwise correlations among CUs. Although synchrony may also arise due to dispersal among populations (Holmes et al. 1994), we did not explicitly model this process because successful dispersal among CUs is assumed, by definition, to be negligible (DFO 2005).

We investigated the impact of reduced productivity, relative to changes in component variability and synchrony, by shrinking relative to its reference value (i.e., the stationary estimate from the longest time series available for a stock; Dorner et al. 2009, Holt and Folkes 2015). Note that is a metric of productivity distinct from ln(recruits per spawner) described in the retrospective analysis with smaller values of resulting in a reduction in the mean recruits‐per‐spawner. To parameterize the low productivity scenario we multiplied each CU's by a scalar (0.65), approximately equal to the mean difference between the median reference estimate of and the minimum estimate of time‐varying for a subset of Fraser River CUs derived from a recursive Bayes stock–recruitment model (time series ranging from brood years 1950–2010, n = 11; C. Michielsens, personal communication, unpublished data). Additionally, we simulated a low productivity regime that incorporated an increased likelihood of extreme “black‐swan” events (Anderson et al. 2017). Specifically, we paired estimates from the low productivity scenario described above with a heavy‐tailed distribution that increases the probability of sampling extreme values. Details and results of the heavy‐tailed scenario are included in Appendix S2.

Management procedure

Total allowable catch (TAC) was calculated each year using a harvest control rule that approximates the total allowable mortality framework used to manage the Fraser River sockeye salmon fishery (Pestal et al. 2011). This framework uses in‐season estimates of return abundance derived from test fisheries to adjust TAC and meet escapement goals specific to each MU (note that MUs include one to six CUs). If in‐season estimates of return abundance exceed escapement goals, TAC is calculated using a fixed maximum target mortality rate of 0.6 (Appendix S1: Fig. S1). Escapement goals vary among years due to cycles and are adjusted upwards to account for en route mortality during upstream migration, as well as spatial overlap between abundant and depleted MUs (Appendix S1). Note that simulated realized exploitation rates deviated from the target due to outcome uncertainty (noted in section Base operating model and described in detail in Appendix S1).

Evaluating model performance

We used a suite of performance metrics to assess how changes in component variability, synchrony, and productivity impact the likelihood of achieving conservation‐ and catch‐based management objectives (Table 3). We present median outputs among simulations for each performance metric, and the 5th and 95th percentiles.

Table 3

Conservation‐ and catch‐based performance metrics, which are presented as median values among trials

Performance metric	Definition
Conservation‐based
Aggregate return abundance	Median aggregate return abundance (millions of individuals) over years and Monte Carlo trials (where aggregate is summation of return abundance across all CUs within a given year).
CU‐specific standardized return abundance	Mean (among CUs) of median (among Monte Carlo trials) CU‐specific returns standardized relative to the low component variability, low synchrony, and reference productivity scenario.
Proportion of MUs meeting escapement goal	The temporal mean proportion of MUs within a return year with return abundance greater than their lower fishery reference point. Under the current harvest control rule, spawner abundances below this level represent a conservation concern that triggers harvest restrictions (Appendix S1 for details).
Proportion of CUs above benchmark delineating sustainable populations	The temporal mean proportion of CUs within a return year with spawner abundance greater than 80% of spawner abundances at maximum sustainable yield (S MSY), defined as an upper biological benchmark for Pacific salmon in Canada (Holt et al. 2009).
Catch‐based
Median catch	Median aggregate catch (i.e., summed across all CUs within a return year; millions of individuals) during the simulation period
Catch stability	The inverse of the temporal coefficient of variation in aggregate catch (i.e., μσ)
Proportion of years above catch threshold	The proportion of years during the simulation period when aggregate TAC was greater than 1,000,000: the level of abundance where managers are able to allocate quota to each quota holder.

We focused our analyses on performance across all CUs and MUs because component variability and synchrony are relevant at the aggregate scale. However, the consequences of increased aggregate variability may vary among components due to intrinsic differences in productivity or carrying capacity, as well as their exposure to harvest in mixed‐CU fisheries. An exhaustive analysis of CU‐specific differences was beyond the scope of this paper; however, to illustrate potential differences, we present simulated changes in median spawner abundance across component variability and synchrony scenarios for two CUs: a red‐status CU, Bowron, and a green‐status CU, Chilko.

We used CU‐specific time series of recruit and spawner abundance to initiate the simulation model (the same data that were used in the retrospective analysis). We ran the simulation over 40 yr, approximately 10 sockeye salmon generations, and simulated each scenario over 1,500 Monte Carlo trials to ensure stable performance metrics. We evaluated each component variability and synchrony scenario under the three productivity scenarios, resulting in 27 distinct operating models.

Results

Retrospective analysis

Generational‐means (i.e., 4‐yr) of Fraser River sockeye salmon ln(recruits/spawner) declined from the late 1980s through 2005, the brood year predominantly responsible for producing the poor return in 2009. Although ln(recruits/spawner) was moderate from 2006 until 2009 (up‐tick in Fig. 1a), it has recently declined again. Generational‐means of aggregate return abundance and aggregate catch increased until the early 1990s before declining (Fig. 1b, c). However, there is substantial interannual variability in raw values due to several abundant cyclic CUs (not shown).

Figure 1

Observed trends in generational‐mean (a) conservation unit (CU)‐specific ln(recruits/spawner), (b) aggregate return abundance, and (c) aggregate catch (millions of fish), as well as estimates of the mean component coefficient of (d) variation, (e) synchrony, and (f) aggregate coefficient of variation for the recruit abundance of 10 CUs. Note that all values are lagged by 4 yr (one generation; a–c) or 12 yr (three generations; d–f) such that the first point in the time series of each panel represents the value for the preceding 4 or 12 yr. In panel a, gray lines represent productivity trends for the 10 individual CUs from which the mean (black line) was calculated. In panels d–f, black lines represent estimates and gray bands represent 95% confidence intervals.

Component coefficient of variation (temporal variability in CU‐specific recruit abundances weighted by the CU's mean abundance) has varied over time, increasing by approximately 50% in the 1950s, 1980s, and 2010s, with considerable uncertainty at the beginning and end of the time series (Fig. 1d). Synchrony in recruit abundances among CUs was relatively high in the 1950s, remained lower and stable from 1960 to 1980, and increased through the 1990s and 2000s, leading to a recent peak for the 2011 brood year (Fig. 1e). Changes in the aggregate CV are the product of component variability and the square root of synchrony and show a gradual increase since the late 1960s followed by a jump to particularly high levels in the first decade of the 2000s (Fig. 1f).

By specifying low, medium, and high values for CU‐specific variability in recruitment residuals, , and the correlation coefficient between CU‐specific recruitment residuals, , (Table 2) we generated a range of component variability and synchrony scenarios (Fig. 2).

Figure 2

Variation in simulated component variability (CVC) and synchrony of recruit abundance (ϕ) under three different component variability scenarios (controlled by scalar; Table 3) and three different synchrony scenarios (controlled by ; Table 3). For the component variability scenarios, is set to its moderate value, while for the synchrony scenarios, the scalar is set to its reference value. Distributions in forward simulations represent median values for each 40‐yr simulation run across 1,500 Monte Carlo trials. Boxes represent the interquartile range (IQR), lower whiskers represent the smallest value greater than 0.5× the lower bound of the IQR, and the upper whisker represent the largest value less than 1.5× the upper bound of the IQR.

The impacts of greater component variability and synchrony varied among performance metrics (Figs. 3 and 4). When productivity was simulated at its reference value (the median retrospective estimate of ) increases in synchrony led to declines in median aggregate return abundance (~7–16% across component variability scenarios; different color symbols in Fig. 3a). Aggregate return abundance declines because the sum of CU‐specific log‐normally distributed recruitments is less right‐skewed, and has a higher median value, when recruitment dynamics are asynchronous than when CU‐specific recruitments are synchronous. Conversely, aggregate return abundance increased by ~9% when component variability was high and synchrony was low (purple symbols Fig. 3a) due to individual CUs experiencing more frequent large positive recruitment deviations. Since high component variability increases the skew of the log‐normal distribution of CU‐specific recruitment, the normalizing effects of asynchronous dynamics are more evident then when component variability is high (Fig. 3a).

Figure 3

Effects of component variability and synchrony on conservation‐based performance metrics for (a, c, e, g) reference and (b, d, f, h) low productivity scenarios on four performance metrics (a, b, aggregate [agg.] return abundance; c, d, CU‐specific standardized [std.] return abundance; e, f, proportion [prop.] of management units [MUs] above their escapement [esc,] goal; and g, h, proportion of CUs above their biological benchmark). Points represent medians and whiskers 90% probability intervals among 1,500 Monte Carlo trials.

Figure 4

Effects of component variability and synchrony on catch‐based performance metrics for (a, c, e) reference and (b, d, f) low productivity scenarios on three performance metrics (a, b, catch; c, d, proportion of years above catch threshold; and e, f, catch stability). Points represent medians and whiskers 90% probability intervals among 1,500 Monte Carlo trials.

Increased component variability had a negative effect on return abundance when it was standardized within each CU (Fig. 3c), rather than summed among CUs (Fig. 3a). Standardized return abundance declined with greater component variability due to density‐dependent effects and serial dependence between recruitment and subsequent spawner abundances, known as time‐series biases (Walters and Martell 2004). Specifically, large spawning events are associated with low recruitments due to overcompensation in the Ricker spawner–recruitment relationship, which results in low spawner abundances in the subsequent year and declines in CU‐specific standardized median return abundances over the time‐series. As component variability increases, these time‐series biases are enhanced and CU‐specific standardized median returns decline as swings from high to low abundances occurred more frequently. Similarly, greater component variability reduced the proportion of MUs above their escapement goal (though declines were modest at low levels of synchrony; Fig. 3e) and the proportion of CUs above their biological benchmark (Fig. 3g).

Increased synchrony also reduced the proportion of MUs exceeding their escapement goals (Fig. 3e), but increased the proportion of CUs exceeding their biological benchmarks (Fig. 3g). The reduction in MU‐level performance with increasing synchrony reflects the reduction in aggregate returns when synchrony is high (Fig. 3a), whereas the increase in CU‐level performance is due to the harvest control rule, which sets exploitation rates based on MU abundance, better tracking the abundance of CUs. The performance of MUs (Fig. 3e) is driven not only by patterns in aggregate abundance (Fig. 3a for synchrony), but also by declines in performance with component variability reflected in CU‐specific abundances (Fig. 3b), since MUs are intermediate in scale between CUs and the total aggregate.

Performance across all conservation‐based performance metrics declined under the low productivity scenario (Fig. 3b, d, f, h). Relative to the reference productivity scenario, variation in return abundance declined (Fig. 3b), while variation in the proportion of MUs above their escapement goal increased (Fig. 3f). Component variability had negligible impacts on CU‐specific standardized returns since abundances were consistently low when productivity was low, and swings between high and low abundances were relatively small (Fig. 3d). Greater component variability increased the probability that biological benchmarks would be exceeded (i.e., the relationship between the proportion of CUs above their benchmark and component variance reversed relative to the reference productivity scenario; Fig. 3g vs. 3h). A similar pattern was also present in the proportion of MUs above their escapement goal, but was restricted to the low and moderate synchrony scenarios (Fig. 3e vs. 3f, purple and green symbols). These increases in performance relative to fixed thresholds reflect the increasing variation in return abundances as component variability increases and a higher proportion of CUs or MUs exceeding those thresholds.

Like conservation‐based performance metrics, the effects of component variability and synchrony varied among catch‐based performance metrics. Under the reference productivity scenario, median aggregate catch declined as synchrony increased (~7–14% across component variability scenarios; different color symbols in Fig. 4a), but increased with greater component variability and low levels of synchrony (purple symbols Fig. 4a). Increasing component variability and synchrony simultaneously led to moderate declines in the proportion of years the target catch threshold was achieved (11% decline between maximum and minimum scenario; Fig. 4c) and severe declines in catch stability (42% decline between maximum and minimum scenario; Fig. 4e).

Median performance declined under the low productivity scenario for all catch‐based performance metrics (Fig. 4b, d, f). Increases in aggregate return abundance at high levels of component variability and low levels of synchrony (Fig. 3b), resulted in increases in aggregate catch abundances and in the proportion of years where the aggregate catch threshold was met (purple symbols Fig. 4b, d).

CU‐specific median spawner abundances declined with greater component variability, although these effects were small relative to declines associated with alternative productivity scenarios (left panels Fig. 5). Since declines with increased component variability are driven by compensatory dynamics and the serial dependence of recruit abundance on prior spawner abundance, the strength of declines varied among CUs due to differences in intrinsic productivity and strength of density dependence. Conversely, increased synchrony resulted in minor increases in CU‐specific median spawner abundance (right panels Fig. 5), due to improved performance of the harvest control rule, which operates at the MU level. The strength of synchrony effects varied among CUs due to differences in sensitivity to fishing mortality within an MU.

Figure 5

Distributions of CU‐specific median spawner abundance (among 1,500 trials) across component variability and synchrony (shading), as well as productivity, scenarios for a depleted (Bowron; orange) and an abundant CU (Chilko; blue). The horizontal dashed line represents each CU's biological benchmark. For the component variability scenarios, is set to its moderate value, while for the synchrony scenarios, the scalar is set to its reference value. Scale of y‐axes differ between CUs.

Discussion

Increases in either the temporal variability of component populations or their synchrony can reduce the stability of an ecological aggregate. In the case of Fraser River sockeye salmon, both component variability and synchrony have increased in recent decades resulting in greater variability in aggregate returns and negative impacts on fisheries. We used a closed‐loop simulation model to demonstrate that simultaneous increases in component variability and synchrony can reduce the probability of achieving a suite of distinct management objectives. However, the magnitude, and even directionality, of such effects vary among performance metrics due to interactions between component variability, synchrony, and productivity at multiple ecological scales. For example, when target exploitation rates are determined by aggregate abundance, greater synchrony can result in exploitation rates that are more likely to be suitable for individual stocks.

Under the reference productivity scenario, aggregate catch stability showed the greatest declines as component variability and synchrony increased. Conversely, the negative impacts of increased synchrony on median aggregate return abundance and catch were relatively modest, while greater component variability actually increased these metrics. Such differences demonstrate how harvests in weakened portfolios will increasingly be driven by boom‐and‐bust cycles even if median catches remain high over longer time horizons. Greater variability in catches may have strong negative impacts on communities that rely on stable fishing opportunities. Indeed, observed declines in aggregate stability have been associated with substantial socioeconomic costs, as documented by the federal inquiry into declines in Fraser River sockeye salmon (Cohen 2012, Peterman and Dorner 2012).

Two threshold‐based performance metrics, the proportion of MUs that exceeded escapement goals and the proportion of years the aggregate catch target was met, also declined as component variability and synchrony increased. Such declines are characteristic of weakened portfolios and arise when asynchronous dynamics no longer buffer aggregates from declines in individual components. Conversely, the proportion of individual CUs above their biological benchmarks exhibited modest increases in high synchrony scenarios due to a lower probability of CUs being depleted by fishing. Specifically, the harvest control rule set exploitation rates based on MU‐specific abundance and when synchrony was high the dynamics of MUs were more representative of their component CUs. An equivalent pattern may occur in systems where “indicator stocks” are used as proxies for the abundance of larger stock complexes; however, this apparent benefit will be sensitive to the management regime in place and should not be considered universal.

While reduced productivity (Peterman and Dorner 2012, Dorner et al. 2018), increased variability (Satterthwaite and Carlson 2015), and increased synchrony (Kilduff et al. 2015, Satterthwaite and Carlson 2015, Freshwater et al. 2018) have been observed in Pacific salmon populations, the cumulative impacts of each process are rarely considered simultaneously. In our simulations, all performance metrics, with the exception of catch stability, were more strongly impacted by reductions in productivity than component variability or synchrony. However, for some performance metrics, declines in productivity reduced or even reversed the negative impacts of greater component variability and synchrony. For example, when productivity and synchrony were low, greater component variability increased the proportion of MUs that exceeded their escapement goals and the proportion of years that the aggregate catch target was met, compared with the reference productivity scenario. Lower productivity also reduced the contrast among levels of component variability or synchrony on catch stability.

The third productivity scenario, which incorporated an increased probability of extreme recruitment events, resulted in particularly low catch stability. Additionally, the heavy tails intermittently generated very large or very small recruitment events, which increased uncertainty across most performance metrics, but resulted in greater median aggregate return and catch abundance at low levels of synchrony.

The results of our simulation analysis have several implications for management strategies that seek to leverage the stability provided by ecological portfolio effects. First, simultaneous declines in productivity and increases in aggregate variability (i.e., the product of component variability and synchrony) will have multiplicative effects that may reduce the probability of achieving management objectives. The magnitude of these impacts, however, depends upon whether changes in aggregate variability are driven by increases in component variability or synchrony and varies among management objectives. Furthermore, ecological scale (i.e., whether the focal unit of an objective is component stocks, the aggregate as a whole, or an intermediate management unit) will affect the degree to which component variability and synchrony affect management objectives. Management actions that fail to account for changes in both component variability and synchrony may underestimate the risk of declines or overestimate the probability of rebuilding ecological aggregates.

Second, we found that even harvest control rules developed to prevent declines below critical thresholds are unlikely to fully buffer aggregates from increases in component variability and synchrony. Although the harvest control rule in our simulation reduced exploitation rates when abundance was low, many performance metrics declined as component variability or synchrony increased. For harvest control rules that are less responsive to declines in abundances (e.g., constant harvest rates common in data‐limited fisheries), high exploitation rates may exacerbate the effects of greater aggregate variability, particularly if productivity is low. In addition, pre‐season forecasts are likely to become less reliable as aggregate variability rises, increasing the probability of overexploitation. A reduction in the apparent effectiveness of a harvest control rule, coupled with declines in catch stability and the probability of meeting catch targets, may result in reduced political support for management strategies that already struggle to balance conservation and socioeconomic objectives.

Natural drivers at various life stages and ecological scales may shape trends in component variability and synchrony in Pacific salmon. Large‐scale processes (i.e., regional or ocean basin scale) may increase component variability if such drivers have themselves become more variable, as well as increase synchrony if their influence has grown relative to local processes (e.g., density‐dependent effects; Walter et al. 2017). Several lines of evidence suggest increases in the component variability and synchrony of Fraser River sockeye salmon stocks may be the result of large‐scale processes. First, the variance of several North Pacific environmental indices associated with Pacific salmon survival has increased since the late 1980s (Sydeman et al. 2013), suggesting interannual variability in habitat quality during marine life stages has also increased. Second, survival during marine residence, when sockeye salmon populations from throughout British Columbia rear in the Gulf of Alaska, has been poor in recent years and has been associated with coherent declines in productivity (Peterman and Dorner 2012, Thomson et al. 2012, Freshwater et al. 2018). Third, environmental forcing at ocean‐basin scales has been linked to greater synchrony among coho salmon (O. kisutch) and Chinook salmon populations from California to southeast Alaska, which have also experienced poor survival in recent years (Kilduff et al. 2015).

Alternatively, anthropogenic disturbance has been identified as a regional process that may lead to greater synchrony in salmon populations. Specifically, stock aggregates that are heavily impacted by hatchery supplementation and hydropower development may experience genetic homogenization that reduces their functional diversity, even if the number of extant stocks is stable (Moore et al. 2010, Satterthwaite and Carlson 2015, Yamane et al. 2018). Although Fraser River sockeye salmon inhabit a heavily developed watershed, hatchery contributions to the aggregate are minimal and hydropower impacts are modest and restricted to four CUs (Grant et al. 2011, Nelitz et al. 2011, COSEWIC 2017). Component variability in Pacific salmon may also be moderated by changes in habitat quality or predator abundance during freshwater residence that influence the reproductive success or mortality of specific stocks (Connor and Pflug 2004, Crozier and Zabel 2006, Crossin et al. 2008, Geist et al. 2008). Local freshwater drivers have likely impacted the component variability of Fraser River sockeye salmon; however it is unlikely that such processes could drive increases in synchrony because trends in freshwater productivity vary among CUs (DFO 2016).

Importantly there may be trade‐offs associated with addressing changes in productivity, variability, and synchrony in real‐world management systems. For example, increasing hatchery supplementation or transporting juveniles may improve the productivity of salmon aggregates, while simultaneously increasing variability and synchrony among populations (Dedrick and Baskett 2018). Conversely, restoring degraded or inaccessible freshwater habitat may reduce synchrony, particularly if colonized by diverse ecotypes, but may provide more modest increases in short‐term productivity. Appropriately parameterized simulation tools can provide a robust methodology for evaluating trade‐offs among management objectives associated with various interventions.

Although the overall patterns we observed are broadly relevant to metapopulations or stock aggregates, the specific effects of greater aggregate variability may be sensitive to several assumptions in our analyses. First, we did not directly model the ecological processes leading to increased component variability and synchrony because there are numerous potential mechanisms, many of which are poorly quantified (e.g., large‐scale environmental drivers). We recognize, however, that the effectiveness of management interventions will depend on their ability to directly address these underlying processes. Second, our conclusions are shaped by the performance metrics selected to assess aggregate status. We principally focused on median outcomes (with the exception of catch stability) over approximately ten generations. Alternative metrics that quantify variability in outcomes will indicate more severe impacts, while those that focus on shorter time horizons may indicate weaker impacts. A critical component of any analysis of the impacts of variability and synchrony will be the careful selection of performance metrics that reflect system‐specific management objectives (Punt et al. 2016). Third, the risks associated with increased variability will vary among components due to differences in initial population size, growth rate, and the relative strength of density dependence. Assemblages that contain a greater proportion of sensitive populations will themselves be more susceptible to weakened portfolio effects.

Systems‐based approaches are increasingly advocated in disciplines such as fisheries biology where accounting for portfolio effects can stabilize the availability of ecosystem services (Link 2018). We demonstrate that systems‐based management efforts may be compromised by increases in either component variability or synchrony, particularly when productivity is depressed, patterns that are consistent with the predicted impacts of climate change. Since stability is maintained by multiple axes of biodiversity (e.g., age structure, species composition; e.g., Schindler et al. 2010, Anderson et al. 2015, Thorson et al. 2018), management efforts may increase long‐term resilience by maintaining intact ecological portfolios.

Click here for additional data file.

Click here for additional data file.