Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

On signals of phase transitions in salmon population dynamics Martin Krkosek and John M. Drake Proc. R. Soc. B 2014 281, 20133221, published 23 April 2014

Supplementary data

"Data Supplement" http://rspb.royalsocietypublishing.org/content/suppl/2014/04/17/rspb.2013.3221.DC1.h tml

References

This article cites 39 articles, 13 of which can be accessed free

Subject collections

Articles on similar topics can be found in the following collections

http://rspb.royalsocietypublishing.org/content/281/1784/20133221.full.html#ref-list-1

ecology (1804 articles) theoretical biology (103 articles)

Email alerting service

Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here

To subscribe to Proc. R. Soc. B go to: http://rspb.royalsocietypublishing.org/subscriptions

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

On signals of phase transitions in salmon population dynamics Martin Krkosˇek1,2 and John M. Drake3,4

rspb.royalsocietypublishing.org

Research Cite this article: Krkosˇek M, Drake JM. 2014 On signals of phase transitions in salmon population dynamics. Proc. R. Soc. B 281: 20133221. http://dx.doi.org/10.1098/rspb.2013.3221

Received: 9 December 2013 Accepted: 25 March 2014

Subject Areas: ecology, theoretical biology Keywords: early warning signals, critical slowing down, phase transitions, portfolio effects, stochasticity, salmon

Author for correspondence: Martin Krkosˇek e-mail: [email protected]

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspb.2013.3221 or via http://rspb.royalsocietypublishing.org.

1 Department of Ecology and Evolutionary Biology, University of Toronto, 25 Willcocks St., Toronto, Ontario, Canada, M5S 3B2 2 Department of Zoology, University of Otago, 340 Great King Street, Dunedin 9016, New Zealand 3 Odum School of Ecology, University of Georgia, 140 East Green Street, Athens, GA 30602-2202, USA 4 Department of Zoology, University of Oxford, South Parks Road, Oxford OX1 3PS, UK

Critical slowing down (CSD) reflects the decline in resilience of equilibria near a bifurcation and may reveal early warning signals (EWS) of ecological phase transitions. We studied CSD in the recruitment dynamics of 120 stocks of three Pacific salmon (Oncorhynchus spp.) species in relation to critical transitions in fishery models. Pink salmon (Oncorhynchus gorbuscha) exhibited increased variability and autocorrelation in populations that had a growth parameter, r, close to zero, consistent with EWS of extinction. However, models and data for sockeye salmon (Oncorhynchus nerka) indicate that portfolio effects from heterogeneity in age-at-maturity may obscure EWS. Chum salmon (Oncorhynchus keta) show intermediate results. The data do not reveal EWS of Ricker-type bifurcations that cause oscillations and chaos at high r. These results not only provide empirical support for CSD in some ecological systems, but also indicate that portfolio effects of age structure may conceal EWS of some critical transitions.

1. Introduction Critical transitions in ecosystems include catastrophic change to a new stable state [1,2], destabilization of populations into cyclical or chaotic dynamics [3,4], and extinction [5]. All are examples of bifurcations [6], where a small change in conditions induces large changes in the steady state and dynamical behaviours of a complex system. Such critical transitions may be foreshadowed by early warning signals (EWS) arising from the interaction of stochastic and nonlinear processes [5,7,8]. Typically, signs of an approaching transition are due to critical slowing down (CSD): the resilience of equilibria (defined as the rate of recovery from perturbation) declines near bifurcations, which allows noise to increase the magnitude, shape and autocorrelation (AC) of system variability [7]. However, some critical transitions do not express EWS [9,10], and the intrinsic heterogeneity of ecological systems may play a role. In population dynamics, portfolio effects describe how heterogeneity in habitat, ecotypes and life history may dampen the influence of environmental stochasticity on population fluctuations [11– 13]. In some cases, however, life-history heterogeneity, via age – or stage–structure, may also amplify population fluctuations owing to resonance arising from interactions between random perturbations and population structure or nonlinear processes [14 –16]. It is therefore likely that heterogeneity in ecological systems modulates EWS of critical transitions by dampening or amplifying the effects of stochastic perturbation. Simple population models can exhibit a rich array of bifurcations as the population growth parameter (r) changes [3]. For the common transcritical bifurcation where r crosses zero, CSD has been shown to provide EWS of extinction in experimental populations [5]. However, some models, such as the Ricker equation [17], Ntþ 1 ¼ Ntexp(r 2 bNt), where r is the population growth parameter, b controls density-dependent mortality and N is abundance [17–19], also have destabilizing bifurcations as r increases [3]. Such instability underlies a leading explanation for how fishing magnifies fluctuations in marine fish populations [20], owing to

& 2014 The Author(s) Published by the Royal Society. All rights reserved.

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

2. Methods (a) Critical slowing down in Beverton– Holt and Ricker models We studied CSD in the stochastic Ricker Ntþ1 ¼ Nt exp(r  bNt þ 1t ), and stochastic Beverton –Holt   Nt , Ntþ1 ¼ exp(r þ 1t ) 1 þ bNt

(2:1)

(2:2)

models via computer simulation. In these models, r is the population growth parameter, b controls density-dependent mortality, N is abundance and 1t is a normally distributed random variable with mean of zero and standard deviation, s, that represents environmental stochasticity [17 – 19,34]. For demonstration of concept, we set a carrying capacity of k ¼ r/b ¼ 1 000 000 and used s ¼ 0.9, and simulated the model over a range of values of r. The qualitative patterns are similar for other choices of s in the range of 0.4 – 1.2 that is commonly seen in the analyses of salmon population dynamics. Other choices of s do not change the predicted relationships between CSD indicators and r. For each simulation, we allowed an initialization period of 450 iterations (this appeared sufficient based on visual inspection) and then analysed the final 50 iterations. CSD indicators were calculated from the final 50 iterations in the same manner as they were calculated from the salmon recruitment data by calculating the CV and lag-1 AC from ln(Nt), the log-transformed abundance. For each value of r (ranging from 0.01 to 3.0), this was repeated 1000 times to generate a distribution of CV and AC estimates from which we estimated the mean and 95% confidence intervals on CV and AC for each value of r. The Ricker model was also used to evaluate the effects of fishing

2

Proc. R. Soc. B 281: 20133221

large—the positive equilibrium in the Beverton–Holt model remains stable for all positive r owing to its eigenvalue l ¼ er , but the Ricker model becomes unstable at r ¼ 2, where its eigenvalue becomes l ¼ 1 2 r ¼ 21, and therefore a flip ( period-doubling) bifurcation occurs [6] where the equilibrium becomes unstable, and there is instead a period-two oscillation. CSD also precedes the flip bifurcation, because the magnitude of l becomes larger as r approaches two, and hence rate of recovery is weaker, and in this case, is mixed with an oscillation, because l is negative. The Ricker model exhibits further period-doubling bifurcations and chaos as r is increased further [3]. However, compared with the simple population models above, some salmon species exhibit heterogeneity in ageat-maturity [31], and such age– or stage –structure may have important consequences for fish population dynamics via nonlinear dynamics or portfolio effects [11,12,32,33]. To analyse this biocomplexity of salmon life histories in relation to CSD indicators, we first analyse simple Ricker and Beverton–Holt population models for the existence of CSD near bifurcations. Then, we study how stage-structured elaborations of Ricker and Beverton– Holt models may affect CSD indicators. Finally, we then examine historical timeseries of population abundance of Pacific salmon to characterize CSD indicators over an among-population gradient in average population growth rates, r. Results of these analyses revealed evidence of EWS of extinction in pink salmon (Oncorhynchus gorbuscha) and chum salmon (Oncorhynchus keta), but not sockeye salmon (Oncorhynchus nerka), where EWS may be obscured by heterogeneity in age-at-maturity.

rspb.royalsocietypublishing.org

interaction between nonlinearity and stochasticity [21–23]. Collectively, these properties suggest that CSD may occur in the vicinity of multiple, but not necessarily all, critical transitions in the parameter space of a complex system. Particularly, we predict that variability and AC in populations with Ricker-type density-dependence would be a U-shaped function of r. Because the Beverton–Holt model, Ntþ 1 ¼ exp(r)[Nt/(1 þ bNt)], does not exhibit this destabilizing bifurcation, variability and AC in populations with density dependence of this type are predicted to decline monotonically as r increases. Ecological studies that examine CSD mostly have been theoretical analyses of simulated dynamics obtained from a stochastic dynamical system under continuous change in a forcing parameter [8,9,24,25]. Data from laboratory populations [5,26,27] and a whole-lake experiment [28] support the theory of EWS, but empirical characterizations of CSD indicators across a parameter range from wild populations are lacking (but see [29]). Such an investigation is data-intensive, because it requires replicated timeseries of ecological dynamics in which to study CSD indicators via at least two different scales of contrast: (i) variation in proximity to bifurcations among replicated timeseries of stationary systems; and (ii) as non-stationary systems transit (or approach) bifurcations within their individual timeseries. While the latter has been empirically investigated in the context of a marine ecosystem regime shift in the Baltic Sea [29], studies of the dynamics of replicated natural ecological systems are particularly rare. Our study addresses this empirical gap using Pacific salmon (Oncorhynchus spp.) as a model system. Populations of Pacific salmon (Oncorhynchus spp.) are known to vary among each other in their average population growth rates and therefore provide an opportunity to investigate if indicators of CSD are present among Pacific salmon populations. To theoretically characterize CSD in salmon populations, we first studied the influence of r on coefficient of variation (CV) and temporal AC in the Ricker and Beverton– Holt fisheries models. Both models contain a transcritical bifurcation at r ¼ 0, but differ at high values of r where the Ricker model has unstable dynamics, but the dynamics of the Beverton–Holt model remain stable. To briefly review these bifurcations, the positive equilibria for Beverton–Holt and Ricker models are N* ¼ (er 2 1)/b and N* ¼ r/b, respectively, and both models also have extinction equilibria at N* ¼ 0. Local stability analysis reveals the stability of the equilibria based on their corresponding eigenvalues, l ¼ f 0 (N  ), where f(N) is the recruitment function (i.e. Ntþ 1 ¼ f(N)), and which measures how a small perturbation near an equilibrium, dt ¼ Nt – N* subsequently changes according to dtþ 1 ¼ ldt [30]. Equilibria are therefore stable if jlj , 1 and unstable if jlj . 1 [18]. The eigenvalues for the positive equilibria are l ¼ er for the Beverton–Holt and l ¼ 1  r for the Ricker model, and both models have an eigenvalue at the extinction equilibrium that is l ¼ f 0 (0) ¼ er . Thus, in both models, when r crosses zero from above, the positive equilibrium is initially stable, and the extinction equilibrium is unstable, the two equilibria collide at r ¼ 0 where they exchange stability, so that the extinction equilibrium becomes stable (i.e. the definition of a transcritical bifurcation [6]). In addition, as r approaches zero, the magnitude of l becomes larger, and therefore the rate of recovery from the perturbation modelled as dtþ 1 ¼ ldt becomes slower, and hence there is CSD. The models differ, however, when r becomes

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

and pre-spawn mortality on the relationship between CV and r (electronic supplementary material, figure S1).

which written in matrix form is Ntþ1 ¼ Bt Nt ,

(2:4)

where B is the transition matrix, F is the abundance of age-1 offspring in a river (fry), that survive to the next year with probability u to become age-2 offspring that have left the river (J; post-smoults). After their first winter at sea, the female smoults either remain at sea as subadults (U1) with probability g1 or return to their river to spawn as 1-sea winter adults (A1) with probability m1. Of the S1 females that spend their second winter at sea, they either remain at sea as subadults (U2) with probability g2 or return to spawn as 2-sea winter adults (A2) with probability m2. All remaining females that spend three winters at sea then return to spawn (A3) with probability m3. All females die after spawning. The fecundity function we used P was w j,t ¼ exp ( fj  A j,t þ 1t ) for Ricker-type dynamics and P w j,t ¼ exp ( fj þ 1t )=(1 þ A j,t ) for Beverton – Holt-type dynamics to represent density-dependent mortality from egg-laying to first year of life. In these functions, fi is the production of age-1 offspring by spawning salmon in class A1, A2 or A3. The parameter 1t is a normally distributed random variable with mean of zero and standard deviation, s, that represents environmental stochasticity. Note that 1t is drawn randomly once for each year and applied to fecundity for each spawning age group, so fecundity among age groups is stochastic but covaries among ages. We simulated the age-structured model (equation (2.3)) with just environmental stochasticity entering through 1t and also through a full-stochastic model where each parameter was randomly selected in each year according to uniform distributions (electronic supplementary material, table S1). In each simulation, we generated a timeseries of 500 recruit observations from which we extracted the last 50 data for further analysis. Recruitment that corresponds to a particular brood year t is Rt ¼ A1,tþ 3 þ A2,tþ 4 þ A3,tþ 5. From the last 50 recruit data, we calculated the CV and lag-1 as well as lag-4 AC of the simulated ln(Rt) recruitment timeseries (where recruitment was log-transformed as in the empirical analysis). Lag-1 AC is expected to have the highest power, but lag-4 was also studied to correspond to the 4 year age at maturity of chum and sockeye salmon. Note that lag-1 for pink salmon is 2 years owing to their data structure. Simulations were repeated 1000 times each over a range of values of f to compare CSD indicators across a range of the population growth parameter. The population growth parameter for the stagestructured model (which relates spawners and the recruits that they produce) is r ¼ ln(R0), where R0 is the dominant eigenvalue of the next-generation matrix F(I 2 T)21, F contains the linearized fecundity elements of B, and T contains the transition probabilities of B so that F þ T ¼ B when N ! 0 [35,36]. The distribution of age-at-maturity was approximately 20% A1, 60% A2 and 20% A3.

(d) Empirical analysis As a preliminary analysis, we compared the magnitude of variation in r at among-population scales relative to nonstationarity in r within individual timeseries to evaluate if there was a relatively large amount of variation among populations in r to investigate variation in CSD indicators. This was done using the tabulated results of Kalman-filtered estimates of within-timeseries variation in r published by Dorner et al. [38] from which we calculated and compared the average standard deviation in r within timeseries to the standard deviation of average r among populations. For each stock, we estimated the CV as well as the lag-1 and lag-4 temporal AC of log-recruitment directly from the timeseries data. To estimate the population growth parameter of each stock, we fitted a Ricker model for spawner – recruit data for each population in the database using the log-transformed version of the Ricker model:   Ri,t ln ¼ ri  bi Si,t þ 1i,t , (2:5) Si,t where Ri,t is the number of recruits in stock i produced by spawners, Si,t, in year t. The population growth parameter of population i is measured as ri, density-dependent mortality within each stock is determined by bi, leaving normally distributed residuals 1j,t, which represent environmental stochasticity, with mean of zero and variance, si 2, to be estimated. The use of the Ricker model here does not imply that the data generating processes were in fact governed by the Ricker equation, only that the population growth parameter can be reasonably estimated from the logarithm of recruits per spawner [40]. Estimates of CV, AC, s and r for each salmon stock then formed the basis of our analysis of CSD indicators for each species. For each species, we looked for a declining or Ushaped relationship between CV and r as well as between AC and r using generalized additive models (GAMs) with Gaussian error and identity link function. First, to place CSD indicators and r on the same scale, we log-transformed the recruitment data prior to calculating the CV and AC values. To account for conventional (non-CSD) effects of environmental variability that might confound the analysis, we compared GAM models containing only s as an explanatory variable versus GAM models that contained s and r as explanatory variables. The models were fit using maximum-likelihood in the mgcv package in R (v. 2.11.1, www.R-project.org) [41] and compared using Akaike’s information criterion (AIC) [42]. For AC, we followed the same procedure. To account for uncertainty in the

3

Proc. R. Soc. B 281: 20133221

To evaluate the influence that heterogeneity in age-at-maturity, common in chum and sockeye salmon, has on CSD indicators, we studied the dynamics of the following matrix population model of female salmon with census time at spawning in autumn: 10 1 1 0 0 0 0 0 0 w1,t w2,t w3,t F F C Bu 0 B J C B 0 0 0 0 0 C CB J C C B B CB U1 C B 0 g1 0 B U1 C 0 0 0 0 CB C C B B C B U2 C ¼ B 0 0 B g2 0 0 0 0 C (2:3) CB U2 C, C B B C C B 0 m1 0 B A1 C B 0 0 0 0 A 1 CB C C B B @0 0 @ A2 A m2 0 0 0 0 A@ A2 A A3 tþ1 m3 0 0 0 A3 t 0 0 0

The data on the population dynamics of Pacific salmon consist of timeseries of spawner – recruit data for 120 stocks: n ¼ 43 for pink salmon (O. gorbuscha); n ¼ 40 for chum salmon (O. keta) and n ¼ 37 for sockeye salmon (O. nerka). These data give the abundance of spawners in each year for each stock and the corresponding abundance of recruits, which is the abundance of adult offspring generated by the spawners and which returned to spawn or were caught in fisheries. The average lengths of the timeseries are 28 years for pink salmon, 27 years for chum salmon and 39 years for sockeye salmon, all distributed among regions of Washington, British Columbia and Alaska. Each pink salmon stock was further divided into odd and even year runs, and analysed separately (pink salmon have a 2-year life cycle, which creates odd and even year populations that are reproductively isolated in each river) [37,38]. A full description of the data as well as the full spawner–recruit dataset itself are available online in reference [38], where among-population variation in average r and within-timeseries non-stationarity in r (associated with a climatic regime shift [39]) is quantified for each population.

rspb.royalsocietypublishing.org

(b) Critical slowing down in age-structured models

(c) Salmon data

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

Ricker

Beverton–Holt

coefficient of variation

0.4 0.3 0.2 0.1

Proc. R. Soc. B 281: 20133221

0 1.0 autocorrelation coefficient

4

rspb.royalsocietypublishing.org

0.5

0.5

0

–0.5

–1.0 0

0.5 1.0 1.5 2.0 2.5 3.0 population growth parameter (r)

0

0.5 1.0 1.5 2.0 2.5 3.0 population growth parameter (r)

Figure 1. Predicted indicators of CSD in unstructured populations show evidence of early warning at the transcritical bifurcation leading to extinction at r ¼ 0. Shown are theoretical predictions from unstructured Ricker and Beverton – Holt models (black lines are the mean and grey regions are 95% confidence regions from simulations of Ricker and Beverton – Holt models) The vertical dotted line is the location of the flip bifurcation (r ¼ 2) in the Ricker model. To demonstrate that the autocorrelation (AC) observed here is due to CSD, the environmental stochasticity used in model simulations was constrained to have no temporal AC (white noise). estimates of r, we weighted the data in the GAM models inversely relative to the standard errors of the r estimates. We report results from weighted and unweighted analyses. To control for potential non-stationarity within the timeseries associated with a 1976– 1977 shift in climate, we repeated the analysis with the data constrained to brood years after 1977. Timeseries prior to the regime shift were frequently too short for estimating r, and so analysis of timeseries prior to the regime shift was not pursued. To summarize, the steps involved in the analysis are as follows: — for each salmon stock, estimate the CV and temporal AC of the log-recruitment timeseries; — for each salmon stock, estimate the population growth parameter (r) and environmental stochasticity (s) by fitting a Ricker model to the spawner– recruit data pairs; — assemble the CV, AC, r and s estimates from each stock into a dataset, one for each species of salmon—pink, chum and sockeye; and — fit and compare GAMs for the relationship between the CSD indicator and s versus the CSD indicator and s and r using maximum-likelihood and AIC.

3. Results First, to describe the among-population and within-population variation in r, we compared among-population variation in average r-values with non-stationarity in r within individual timeseries from a Kalman filter analysis of these data published in Dorner et al. [38]. Among-stock variation in the

average growth parameter (SDr ) was larger than the average of non-stationary variation in r within timeseries (measured as the average among stocks in their within timeseries standard deviation of the growth parameter, SDr ): pink salmon SDr ¼ 0.57 versus SDr ¼ 0.24; chum salmon SDr ¼ 0.43 versus SDr ¼ 0.35; sockeye salmon SDr ¼ 0.44 versus SDr ¼ 0.31. Thus, the key ingredient we need to study CSD indicators—variation in r that spans parameter space near and far from bifurcation—is present in the timeseries. In Ricker and Beverton–Holt models with no age structure, variability and temporal AC increased with proximity to the transcritical bifurcation at r ¼ 0 (figure 1). Further, instability in the Ricker model at the flip bifurcation at r ¼ 2 [3,6] also induced CSD (figure 1), giving rise in simulations to the hypothesized U-shaped relationship between variance and r. The effect of fishing, or other pre-spawn mortality, which occur in the salmon populations we studied, can affect the location of the minimum and curvature of the overall U-shaped pattern (electronic supplementary material, figure S1). Reddening the noise term (i.e. increasing AC) did not change the overall predicted relationships between CSD indicators and r, but it did cause the U-shaped relationship between CV and r to shift rightward and, unsurprisingly, also increased the overall level of AC (electronic supplementary material, figure S2). Adding stage-structure to the population model to account for variation in age-at-maturity—common to chum and sockeye salmon [31]—can reduce EWS near r ¼ 0 (where r ¼ ln(R0) and R0 is the dominant eigenvalue of the

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

Ricker

Beverton–Holt

5

coefficient of variation

rspb.royalsocietypublishing.org

0.5 0.4 0.3 0.2 0.1

Proc. R. Soc. B 281: 20133221

0

autocorrelation coefficient

1.0

0.5

0

–0.5 lag-1

lag-1

lag-4

lag-4

–1.0

autocorrelation coefficient

1.0

0.5

0

–0.5

–1.0 0

0.5

1.0

1.5

2.0

2.5

population growth parameter (r)

3.0

0

0.5

1.0

1.5

2.0

2.5

3.0

population growth parameter (r)

Figure 2. Predicted indicators of CSD in stage-structured populations show that variable age-at-maturity dampens the early warning of extinction. Shown are theoretical predictions from stage-structured Ricker and Beverton –Holt models that contain variation in age-at-maturity (black lines are the mean and grey regions are 95% confidence regions from simulations of Ricker and Beverton– Holt models). The vertical dotted line is the location of the bifurcation (near r ¼ 2). Simulated data are from the full stochastic model (all parameters randomly drawn per time step), though results of simulations where stochasticity entered only through 1t yielded very similar results. To demonstrate that the autocorrelation (AC) observed here is due to CSD, the environmental stochasticity used in model simulations was constrained to have no temporal AC (white noise). linearized next-generation matrix; see Methods; figure 2) although CSD owing to instability of Ricker-type dynamics at higher r persists. The effect of heterogeneity in age-atmaturity occurred in both the fully stochastic model (where all parameters were stochastic) as well as the model where stochasticity entered only through fecundity and fecundity among age classes was 100% correlated (electronic supplementary material, figure S3). Reddening the noise in the stage-structured model had the effect of slightly recovering patterns of CSD indicators and r observed in the unstructured model (electronic supplementary material, figures S4 and S5). For the empirical analyses, the most reliable empirical results are those based on post-1977 data (after the regime

shift) with data weighted by uncertainty in r estimates (column D of table 1), because these minimize influences of non-stationarity within timeseries and account for uncertainty. For the bifurcation at r ¼ 0, pink salmon, which do not exhibit variable age-at-maturity, provided the most data near the bifurcation and revealed a negative relationship in amongpopulation variation between CV and r and between AC and r, as predicted (table 1 and figure 3; electronic supplementary material, figure S6). However, sockeye salmon, which exhibit life-history heterogeneity via variable age-at-maturity [31], did not reveal evidence of CSD near r ¼ 0 at the among-population level consistent with predictions from the stage-structured models (figure 3). Chum salmon, which also

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

pink

chum

sockeye

–1 0 1 2 3 4 population growth parameter (r)

–1 0 1 2 3 4 population growth parameter (r)

–1 0 1 2 3 4 population growth parameter (r)

6

coefficient of variation

rspb.royalsocietypublishing.org

0.5 0.4 0.3 0.2 0.1 0

Proc. R. Soc. B 281: 20133221

autocorrelation coefficient

1.0 0.5 0 –0.5 –1.0

Figure 3. Empirical estimates for CSD indicators—coefficient of variation and autocorrelation (AC) coefficient—across the population growth rate from stocks of pink, chum and sockeye salmon for brood years after the 1977 regime shift. Note that pink salmon do not have variable age-at-maturity, and hence no age structure, whereas chum and sockeye salmon have variation in age at maturity with a typical mode of 4 years. AC estimates are lag-1 for pink salmon (equates to 2 years) and 4 years for chum and sockeye. Solid lines are plots of the fitted GAM models (+2 s.e. in dashed lines), but note that these are not the expected least-squares fits, because the model is marginalized over a third axis corresponding to s.

Table 1. Model selection statistics of CSD indicators in salmon recruitment—the coefficient of variation (CV) and temporal autocorrelation (AC, with lag 1 (2 years) for pink salmon and lag 1 and 4 years for chum and sockeye). (The DAIC values are the differences in AIC values for generalized additive models of the indicator as a function of the population growth parameter, r, and conventional environmental stochasticity, s. Four analyses are shown: (A) full datasets with no weighting; (B) full dataset but weighted inversely with the standard errors of r; (C) data for brood years after the 1977 regime shift with no weighting; and (D) data for brood years after the 1977 regime shift and data weighted as in (B).) DAIC species

indicator

model

pink

CV

s sþr s sþr s sþr s sþr s sþr s sþr s sþr s sþr

AC chum

CV AC-1 AC-4

sockeye

CV AC-1 AC-4

A

B

C

D

3.65

3.92

6.85

6.24

0 20.99

0 23.03

0 21.96

0 22.40

0

0

0

0

2.11 0

2.43 0

2.11 0

6.71 0

0.30 0

0.42 0

0.29 0

21.99 0

8.42 0

11.22 0

3.49 0

3.70 0

5.14

1.89

2.95

21.36

0 1.58

0 0.23

0 21.99

0 6.5

0 2.21

0 1.21

0

0

0

0 22 0 21.49 0

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

4. Discussion

7

Proc. R. Soc. B 281: 20133221

Our theoretical predictions of EWS of ecological phase transitions for salmon population dynamics indicate that heterogeneity in age-at-maturity may reduce EWS of extinction for both Ricker and Beverton–Holt-type dynamics. The mechanism is that, within one recruitment year, the effects of environmental stochasticity are averaged over several age classes. This averaging effectively dampens interannual variability caused by environmental stochasticity akin to how a diversified set of assets dampens the overall fluctuations of an investment portfolio. This so-called portfolio effect is relatively new to population dynamics [12] and is similar also to the statistical inevitability of diversity– stability relationships in community ecology [43]. Such portfolio effects also dampen the variance in recruitment during Ricker-type destabilizing bifurcations when the population growth parameter, r, is large and increasing, although the qualitative nature of increasing CV as r is increased through the bifurcation is preserved for the age-structured Ricker model. Thus, there is a variety of qualitative forms that CSD indicators in salmon population dynamics may theoretically take in relation to r that is mediated by the heterogeneity that is inherent to the species life history. Stage –structure or multidimensional systems, in general, may also affect the existence of EWS in other ways. In particular, EWS owing to CSD near a catastrophic bifurcation in a simple predator –prey model is detectable only in the juvenile component of the prey population, but not in the adult [44]. Most fishery datasets only track abundance or catch of adults, and so in cases where the model in reference [44] is representative, the bifurcation may not be presaged, because the relevant variable has not be measured. Note that the mechanism in reference [44] is not a consequence of statistical averaging of asynchronous variability among population components (i.e. portfolio effects) but rather via a different mechanism: the orientation of eigenvectors and the value of their corresponding eigenvalues that characterize the dynamical system near the bifurcation point allow some variables to not express CSD [44]. Thus, our results provide a new mechanism that adds to EWS theory in that EWS in multidimensional systems are not universal, but rather dependent on the structure of the system. Comparing CSD indicators across populations in relation to their average r-values revealed patterns of EWS of extinction that were consistent with theoretical predictions in two of the three species. For pink salmon, which have minimal variation in age-at-maturity, the data revealed evidence for CSD, both for CV and AC indicators. For chum salmon, we also observed evidence for EWS of extinction in CV as well as lag-4 AC, even though chum salmon have variability in age-at-maturity. By contrast, for sockeye salmon, which also have variable age-at-maturity, there was no evidence for EWS of extinction consistent with theoretical predictions.

Our simulation study suggests that AC rises sooner than CV when r approaches zero, suggesting this may be a more reliable EWS of extinction. The empirical results on CSD were also stronger for AC than CV, suggesting that AC may be a more reliable EWS of extinction. Evidence that heterogeneity in age-at-maturity may diminish EWS of extinction mostly rests on our theoretical results; the power to detect effects of heterogeneity in age-at-maturity in data is limited by few replicates of species that vary in heterogeneity in age-at-maturity as well as variation among species in the quantity of data at low r-values. The results did not reveal evidence for CSD at high values of r as was predicted by the Ricker model, but not by the Beverton– Holt model. However, this result is not strong evidence that instability of fish population dynamics at high r-values does not occur. For salmon population dynamics, potential overcompensation that would lead to Ricker-type dynamics could occur through overcrowding of spawning habitat in natural populations. Such dynamics may rarely occur in the data we studied, because fisheries harvests are likely to often reduce abundance of spawners, so that overcompensation occurs infrequently. This idea is consistent with our simulations where fishing makes CSD indicators flatten out at high r-values making CSD less detectible there, but also makes CSD at low r-values more detectable. More generally, instability of fish stocks owing to increased r is a main explanation for increased amplitude of fluctuations of exploited fish species [20,45]. The main objective of our study was to test for predicted directional relationships between CSD indicators and population growth—a qualitative test. Although this qualitative objective was achieved, we failed to observe clear quantitative correspondence between the observed and predicted values of CSD indicators. Why? One possibility is that other factors, such as variation in the local magnitude of environmental stochasticity may cause quantitative estimates of CSD indicators from the theoretical models to differ from observed estimates. Additionally, parameter values used in the theoretical models were chosen to reflect known variation in those parameters and to generate qualitatively testable predictions on the relationships between CSD indicators and population growth rates. The observed directional relationships between CSD indicators and population growth rates were indeed consistent with predictions, but not in all cases. Finally, because the empirical results are based on correlations from field data and not experimental manipulations, it is possible that there are other mechanisms that underlie the results. We did not analyse the data for within-timeseries nonstationarity in r for EWS of the 1977 climate-associated regime shift [39]. We did not expect to see EWS of the regime shift in the salmon data, because the regime shift is an external environmental change that occurred abruptly relative to salmon population dynamics. Indeed, even if this regime shift caused r to approach or cross a bifurcation, then we would not expect to observe EWS of the transition, because the change would have been abrupt. Furthermore, according to theory one only expects to see EWS of the regime shift expressed via changing r-values in salmon-recruitment dynamics if the shift is a smooth transition over multiple spawning periods. These arguments highlight the need to better understand how the rate of environmental change and the associated speed of approach to bifurcation affect the expression of CSD.

rspb.royalsocietypublishing.org

exhibit variation in age-at-maturity did, however, reveal evidence for a negative relationship between CV and r at low r-values as well as between lag-4 AC and r at low r-values (table 1 and figure 3). There was no evidence for a positive relationship between CV and r nor AC and r for any species (table 1 and figure 3), as would be expected owing to Ricker-type instability at high values of r.

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

Acknowledgements. We thank Randall Peterman for contributing the Pacific salmon spawner–recruit dataset as well as RM and referees for helpful comments that improved the work. Funding statement. This work was supported by funding from the Natural Sciences and Engineering Research Council of Canada (M.K.) and the Odum School of Ecology, the Sarah Moss Fellowship at the University of Georgia, and a grant from the Leverhulme Trust (J.M.D.).

References 1.

2.

3.

Holling CS. 1973 Resilience and stability of ecological systems. Annu. Rev. Ecol. Syst. 4, 1– 23. (doi:10.1146/annurev.es.04.110173.000245) Scheffer M, Carpenter S, Foley JA, Folke C, Walker B. 2001 Catastrophic shifts in ecosystems. Nature 413, 591–596. (doi:10.1038/35098000) May RM. 1974 Biological populations with nonoverlapping generations: stable points, stable cycles, and chaos. Science 186, 645 –647. (doi:10.1126/science.186.4164.645)

4.

5.

6.

Fussmann GF, Ellner SP, Shertzer KW, Hairston NG. 2000 Crossing the Hopf bifurcation in a live predator–prey system. Science 290, 1358–1360. (doi:10.1126/science.290.5495.1358) Drake JM, Griffen BD. 2010 Early warning signals of extinction in deteriorating environments. Nature 467, 456 –459. (doi:10.1038/nature09389) Strogatz SH. 1994 Nonlinear dynamics and chaos. Cambridge, UK: Perseus Books Publishing.

7.

8.

9.

Scheffer M et al. 2009 Early-warning signals for critical transitions. Nature 461, 53 –59. (doi:10.1038/nature08227) Carpenter SR, Brock WA. 2006 Rising variance: a leading indicator of ecological transition. Ecol. Lett. 9, 308–315. (doi:10.1111/j.1461-0248.2005.00877.x) Hastings A, Wysham DB. 2010 Regime shifts in ecological systems can occur with no warning. Ecol. Lett. 13, 464–472. (doi:10.1111/j.1461-0248.2010. 01439.x)

8

Proc. R. Soc. B 281: 20133221

chum salmon (but not sockeye) that were consistent with predictions of CSD at low values of r in theoretical models. However, we did not directly test for changes in CSD indicators within individual timeseries during a bifurcation, and as such it departs from other studies that typically change a parameter, so that the system transits a bifurcation and then analyses data for EWS of the bifurcation [8,9,24,25]. By contrast, we empirically characterized parameter space of salmon population dynamics using variation among distinct populations in their average distance from locations in parameter space where bifurcations are predicted by theory. EWS theory, according to our simulation results, yields testable predictions at this among-population scale that were evident in two of the three species we studied. EWS of critical transitions may be useful for anticipating changes in system state or dynamics, including extinction [5]. Of course, if one has adequate data and well-validated models, then the approach to bifurcations may be anticipated directly by trends in parameter estimates alone—for example through state–space modelling [38,47]. In most situations, however, sufficient models and data will not be available, so that any inferences about impending system change must be based on more limited information. Our results indicate that partial data (i.e. recruitment data alone as opposed to full spawner–recruit data) in the absence of detailed dynamical models may be helpful, particularly for anticipating extinction in species with simple life cycles. We have also shown, however, that species with intact heterogeneity owing to life-history variation may conceal EWS of extinction. Thus, extinction may be less predictable in ecological systems with intact heterogeneity than in systems that have been simplified by human activities such as exploitation or habitat modification. However, although intact heterogeneity may conceal EWS of extinction it may also reduce the probability of extinction owing to reduced variability at low r. Furthermore, populations in degraded habitats may be more variable owing to the erosion of portfolio effects [12] as well as CSD associated with a decline in r. In conclusion, our results may not only be the first observations of CSD across a parameter range from natural populations, but also caution that age–structure effects may conceal EWS of some critical transitions.

rspb.royalsocietypublishing.org

While we compared empirical patterns of CSD with theoretical predictions from Beverton–Holt and Ricker models (and their stage-structured elaborations), our analysis does not require that these models were the data generating processes of salmon population dynamics. The true data generating processes are probably far more complex. Rather, these models illustrate and contrast how CSD indicators may be qualitatively related to the population growth parameter for unstructured and stage-structured populations, and for compensatory and overcompensatory density dependence. Similarly, the parameter values used in simulations were not empirically derived, but rather chosen to give representative values for fecundity, smoult survival and distribution of age at maturity [31]. Nevertheless, these models are routinely used in contemporary fisheries science and population ecology [16,18,19,46]. Indeed, a key feature of the theory of CSD is that its predictions are qualitative— a full understanding of the nonlinear dynamical mechanisms is not required for the theory to be applied. While these results indicate that a rise in CSD indicators may be useful as EWS, it is important to note that the practicality of this approach also depends on the pace of change in r, which must be slow compared with the timescale of recruitment. Furthermore, our results pertain to bifurcations away from a stable equilibrium and not among unstable equilibria, where dynamics are more complex and the utility of EWS is far from clear. For EWS theory to be applied practically, one needs to look within non-stationary timeseries using a movingwindow analysis to detect temporal increases in CSD indicators. Such an approach is more data intensive than the salmon data can support. This is because many recruitment observations are needed to calculate CSD indicators, which takes time, and it is possible that fast environmental changes could cause bifurcations to occur without the accumulation of data necessary to reveal EWS. This is partly the basis of why we do not expect salmon population dynamics to reveal EWS of shifts in the Pacific decadal oscillation. As such, age structure is not the only mechanism by which EWS may be eroded; fast environmental change is another. Overall, the practical use of early warning systems will depend on the speed at which the critical transition is approached, the severity of the nonlinearity in the approach to bifurcation, the magnitude of intrinsic noise, and the precision of the data. However, analyses from laboratory populations indicated that CSD could be detected as early as eight generations prior to the critical transition in an experimentally deteriorating model ecosystem [5]. Our main empirical results revealed patterns of CV and AC among timeseries from replicated wild populations of pink and

Downloaded from rspb.royalsocietypublishing.org on November 6, 2014

37.

38.

39.

40.

41.

42. 43.

44.

45.

46.

47.

output. Am. Nat. 172, 128–139. (doi:10.1086/ 588060) Pyper BJ, Mueter FJ, Peterman RM, Blackbourn DJ, Wood CC. 2001 Spatial covariation in survival rates of Northeast Pacific pink salmon (Oncorhynchus gorbuscha). Can. J. Fish. Aquat. Sci. 58, 1501–1515. (doi:10.1139/f01-096) Dorner B, Peterman RM, Haeseker SL. 2008 Historical trends in productivity of 120 Pacific pink, chum, and sockeye salmon stocks reconstructed by using a Kalman filter. Can. J. Fish. Aquat. Sci. 65, 1842– 1866. (doi:10.1139/F08-094) Mantua NJ, Hare SR, Zhang Y, Wallace JM, Francis RC. 1997 A Pacific interdecadal climate oscillation with impacts on salmon production. Bull. Am. Meteorol. Soc. 78, 1069 –1079. (doi:10.1175/15200477(1997)078,1069:APICOW.2.0.CO;2) Myers RA, Bowen KG, Barrowman NJ. 1999 Maximum reproductive rate of fish at low population sizes. Can. J. Fish. Aquat. Sci. 56, 2404– 2419. Wood SN. 2006 Generalized additive models: an introduction with R. Boca Raton, FL: Chapman and Hall/CRC. Burnham KP, Anderson DR. 2002 Model selection and multimodel inference. New York, NY: Springer. Doak DF, Bigger D, Harding EK, Marvier MA, O’Malley RE, Thomson D. 1998 The statistical inevitability of stability-diversity relationships in community ecology. Am. Nat. 151, 264–276. (doi:10.1086/286117) Boerlijst MC, Oudman T, De Roos AM. 2013 Catastrophic collapse can occur without early warning: examples of silent catastrophes in structured ecological models. PLoS ONE 8, e62033. (doi:10.1371/journal.pone.0062033) Hsieh CH, Reiss CS, Hunter JR, Beddington J, May RM, Sugihara G. 2006 Fishing elevates variability in the abundance of exploited species. Nature 443, 859–862. (doi:10.1038/nature05232) Melbourne BA, Hastings A. 2008 Extinction risk depends strongly on factors contributing to stochasticity. Nature 454, 100–103. (doi:10.1038/ nature06922) Peterman RM et al. 2010 Synthesis of evidence from a workshop on the decline of Fraser River sockeye. June 15 –17, 2010. Vancouver, A report to the Pacific Salmon Commission.

9

Proc. R. Soc. B 281: 20133221

23. Sugihara G et al. 2011 Are exploited fish populations stable? Proc. Natl Acad. Sci. USA 108, E1224 –E1225. (doi:10.1073/pnas.1112033108) 24. Biggs R, Carpenter SR, Brock WA. 2009 Turning back from the brink: detecting an impending regime shift in time to avert it. Proc. Natl Acad. Sci. USA 106, 826 –831. (doi:10.1073/pnas.0811729106) 25. Carpenter SR, Brock WA, Cole JJ, Kitchell JF, Pace ML. 2008 Leading indicators of trophic cascades. Ecol. Lett. 11, 128–138. 26. Veraart AJ, Faassen EJ, Dakos V, van Nes EH, Lurling M, Scheffer M. 2012 Recovery rates reflect distance to a tipping point in a living system. Nature 481, 357. 27. Dai L, Vorselen D, Korolev KS, Gore J. 2012 Generic Indicators for loss of resilience before a tipping point leading to population collapse. Science 336, 1175 –1177. (doi:10.1126/science.1219805) 28. Carpenter SR et al. 2011 Early warnings of regime shifts: a whole-ecosystem experiment. Science 332, 1079 –1082. (doi:10.1126/science.1203672) 29. Lindergren M, Dakos V, Gro¨ger JP, Ga˚rdmark A, Kornilovs G, Otto SA, Mo¨llmann C. 2012 Early detection of ecosystem regime shifts: a multiple method evaluation for management application. PLoS ONE 7, e38410. (doi:10.1371/journal.pone. 0038410) 30. Murray JD. 1989 Mathematical biology. Berlin, Germany: Springer. 31. Groot C, Margolis L. 1991 Pacific salmon life histories. Vancouver, Canada: UBC Press. 32. Botsford LW, Holland MD, Samhouri JF, White JW, Hastings A. 2011 Importance of age structure in models of the response of upper trophic levels to fishing and climate change. ICES J. Mar. Sci. 68, 1270 –1283. (doi:10.1093/icesjms/fsr042) 33. Hsieh CH, Yamauchi A, Nakazawa T, Wang WF. 2010 Fishing effects on age and spatial structures undermine population stability of fishes. Aquat. Sci. 72, 165– 178. (doi:10.1007/s00027-009-0122-2) 34. Beverton RJH, Holt SJ. 1993 On the dynamics of exploited fish populations. Fish Fish. Ser. 11, 1–537. 35. de-Camino-Beck T, Lewis MA. 2007 A new method for calculating net reproductive rate from graph reduction with applications to the control of invasive species. Bull. Math. Biol. 69, 1341–1354. (doi:10.1007/s11538-006-9162-0) 36. de-Camino-Beck T, Lewis MA. 2008 On net reproductive rate and the timing of reproductive

rspb.royalsocietypublishing.org

10. Brock WA, Carpenter SR. 2010 Interacting regime shifts in ecosystems: implication for early warnings. Ecol. Monogr. 80, 353–367. (doi:10.1890/09-1824.1) 11. Hilborn R, Quinn TP, Schindler DE, Rogers DE. 2003 Biocomplexity and fisheries sustainability. Proc. Natl Acad. Sci. USA 100, 6564 –6568. (doi:10.1073/pnas. 1037274100) 12. Schindler DE, Hilborn R, Chasco B, Boatright CP, Quinn TP, Rogers LA, Webster MS. 2010 Population diversity and the portfolio effect in an exploited species. Nature 465, 609–612. (doi:10.1038/ nature09060) 13. Greene CM, Hall JE, Guilbault KR, Quinn TP. 2010 Improved viability of populations with diverse lifehistory portfolios. Biol. Lett. 6, 382 –386. (doi:10. 1098/rsbl.2009.0780) 14. Stenseth NC, Bjornstad ON, Falck W, Fromentin JM, Gjosaeter J, Gray JS. 1999 Dynamics of coastal cod populations: intra- and intercohort density dependence and stochastic processes. Proc. R. Soc. Lond. B 266, 1645–1654. (doi:10.1098/rspb.1999.0827) 15. Guill C, Drossel B, Just W, Carmack E. 2011 A threespecies model explaining cyclic dominance of Pacific salmon. J. Theor. Biol. 276, 16 –21. (doi:10.1016/j. jtbi.2011.01.036) 16. Krkosek M, Hilborn R, Peterman RM, Quinn TP. 2011 Cycles, stochasticity and density dependence in pink salmon population dynamics. Proc. R. Soc. B 278, 2060–2068. (doi:10.1098/rspb.2010.2335) 17. Ricker WE. 1954 Stock and recruitment. J. Fish. Res. Bd Can. 11, 559– 623. (doi:10.1139/f54-039) 18. Kot M. 2001 Elements of mathematical ecology. Cambridge, UK: Cambridge University Press. 19. Hilborn R, Walters CJ. 1992 Quantitative fisheries stock assessment. New York, NY: Chapman and Hall. 20. Anderson CNK, Hsieh CH, Sandin SA, Hewitt R, Hollowed A, Beddington J, May RM, Sugihara G. 2008 Why fishing magnifies fluctuations in fish abundance. Nature 452, 835–839. (doi:10.1038/ nature06851) 21. Shelton AO, Mangel M. 2011 Reply to Sugihara et al.: the biology of variability in fish populations. Proc. Natl Acad. Sci. USA 108, E1226. (doi:10.1073/ pnas.1115765108) 22. Shelton AO, Mangel M. 2011 Fluctuations of fish populations and the magnifying effects of fishing. Proc. Natl Acad. Sci. USA 108, 7075– 7080. (doi:10. 1073/pnas.1100334108)

On signals of phase transitions in salmon population dynamics.

Critical slowing down (CSD) reflects the decline in resilience of equilibria near a bifurcation and may reveal early warning signals (EWS) of ecologic...
680KB Sizes 2 Downloads 4 Views