Introduction

Forecasting complex trajectories of marine resources is essential to fishery management and one of the major challenges of our time (Schindler and Hilborn, 2015;Ye et al., 2015). An important factor to be considered in fishery management is the stock-recruitment relationship (SRR), which serves as a basis for any stock assessment procedure to ultimately calculate reference points (Hilborn, 2002;ICES, 2017ICES, , 2018)). SRRs are based on the assumption that recruitment (the number of fish that enter the adult population) is directly related to adult stock size (Kraus et al., 2000;Jennings et al., 2001). Parametric approaches, such as the Ricker model, were developed around the 1950s (Ricker, 1954) and in some cases still represent the method of choice in stock assessments (ICES, 2017). These models are very specific in the type of functional response curve to describe the SRR, and are linear, in the sense that, the relationship between recruitment and biomass can be linearized through log-transformation (Ye et al., 2015). However, they often fail to capture the high variability in recruitment data and this has led to questioning the existence of the relationship itself (Szuwalski et al., 2015;Britten et al., 2016;Perlala et al., 2017). The fit of the SRR is sometimes so poor, that short-term predictions of spawning-stock biomass (SSB) are conducted using an average of recruitment over a particular number of years, instead of a SRR model (Deyle et al., 2018). Both approaches, using average recruitment or a parametric model, assume that natural systems behave in a linear way, which may lead to biased fishery management decisions when stocks show complex dynamics such as aperiodic chaos, non-linearity, or nonstationarity (Ye and Sugihara, 2016;Perlala et al., 2017;Deyle et al., 2018).

Chaos and non-stationary dynamics are pervasive in natural systems and characterize many marine ecosystems and populations (May and Oster, 1976;Scheffer et al., 2001;Mo ¨llmann et al., 2015). These dynamics emerge from the inherent complexity of nature, governed by a multitude of factors (Ye et al., 2015;Deyle et al., 2016;Tu et al., 2018). Assuming linearity and stability in recruitment models can, thus, result in wrong stock predictions (Glaser et al., 2014;Ye et al., 2015). As a consequence, new nonparametric modelling frameworks were developed to predict stock trajectories accounting for state-dependent and chaotic behaviour, such as the empirical dynamic modelling (EDM) framework (Sugihara et al., 2012;Ye et al., 2015;Deyle et al., 2018). EDM is a minimal assumptive approach based on time-series observations, which reconstructs the temporal dynamics of a system by constructing a so-called attractor manifold (Sugihara et al., 2012;Ye et al., 2015). EDM is able to predict the future system trajectory based on its past dynamics (Ye et al., 2015;Deyle et al., 2018), thus accounting for state-dependent dynamics (Sugihara, 1994). This approach, and in particular multivariate simplex projection (MSP) has been applied to predict non-linear fish recruitment dynamics in a range of studies, and has also been applied directly to management, e.g. for the menhaden stocks along the East Coast of the United States (Perretti et al., 2015;Ye et al., 2015;Deyle et al., 2018).

Another non-parametric approach suitable for modelling state-dependent and discontinuous recruitment dynamics is the stochastic cusp model (SCM), which is based on catastrophe theory (Zeeman, 1976;Thom, 1977;Grasman et al., 2009;Petraitis and Dudgeon, 2016;Sguotti et al., 2019). Here, a state variable z (for instance recruitment), depends on two control variables alpha and beta. The model allows z to move from a state A (e.g. high recruitment) to a state B (e.g. low recruitment) following either a continuous or discontinuous path (Diks and Wang, 2016). SCM has been widely applied to economic and behavioural studies (van der Maas et al., 2003;Diks and Wang, 2016), but to a lesser degree to marine ecological studies (Jones and Walters, 1976;Jones, 1977;Petraitis and Dudgeon, 2015;Sguotti et al., 2019).

Another point often neglected in recruitment prediction is the effect of multiple external drivers and potential interactions such as predation, competition, and environmental variables (Myers et al., 1995;Brander, 2005;Ottersen et al., 2006;Stiasny et al., 2016). However, in multiple cases the relationship between recruitment and environment can be spurious, non-linear, or non-stationary, and therefore is often not considered in stock assessments (Myers, 1998;Perlala et al., 2017). Traditional stockrecruitment models, which often are parametric models, assuming fixed parameters, usually fail to correctly incorporate the environmental information, because they often just consider additive effects of SSB and climate variables. Instead, nonparametric models such as MSP and SCM, can model interactions between the different drivers (i.e. biomass and climate variables) and thus may be able to integrate the climate information correctly (Ye et al., 2015;Deyle et al., 2018;Sguotti et al., 2019). This is important because for effectively predicting the status of living marine resources the integration of environmental variables is becoming crucial given the widespread impacts of climate change on ecosystems and marine resources such as commercially important fish (Britten et al., 2016;Gaines et al., 2018).

Atlantic cod (Gadus morhua) is an iconic species from ecological, cultural, and economic points of view (Myers et al., 1996). In recent decades, most North Atlantic cod stocks have collapsed, followed by prolonged periods of no recovery even after the application of strict management measures (e.g. fishing moratoria; Myers et al., 1996;Hutchings, 2000;Hutchings and Rangeley, 2011;Frank et al., 2016;Sguotti et al., 2019). This failed recovery of Atlantic cod stocks suggests the presence of discontinuous dynamics and hysteresis (Frank et al., 2011;Steneck et al., 2011;Sguotti et al., 2019). Eastern and western Atlantic stocks differ in life-history traits, exploitation trajectories and recovery potential (Po ¨rtner et al., 2008;Wang et al., 2014;Frank et al., 2016). Indeed, stocks in the West collapsed more abruptly compared to stocks in the East, which on average show more gradual declines (Frank et al., 2016). Cod recruitment is highly state-dependent, depending on the dimension of the stock and environment conditions. Recruitment is fundamental to Atlantic cod recovery (Myers and Barrowman, 1996;Brander, 2005) and influenced by climate change (Myers and Drinkwater, 1989;Planque and Fre ´dou, 1999;Stige et al., 2006;Po ¨rtner et al., 2008;Pershing et al., 2015). We here used stock assessment data from 20 Atlantic cod stocks to (i) investigate whether cod recruitment can be best described by the parametric Ricker model, by the nonparametric, "discontinuous" SCM, or by the non-parametric, state-dependent MSP approach, and (ii) test whether the model's predictive power can be improved when including environmental variables. We show that the adoption of a multi-model approach should be considered when modelling stocks presenting different dynamics.

Material and methods

Data

We used recruitment (i.e. number of fish for a particular age and stock that recruit to the adult biomass in thousands, R) and SSB (i.e. biomass of mature fish in tonnes) data derived from stock assessments of 20 Atlantic cod stocks (Figure 1, Supplementary Figure S1). Data were provided by the International Council for the Exploration of the Sea (ICES), the National Oceanic and Atmospheric Administration of the United States (NOAA), the northwest Atlantic Fisheries Organization (NAFO), the Department of Fisheries and Ocean in Canada (DFO), and by personal communication (Supplementary Table S1). Recent assessments for cod stocks in the Kattegat, the western Baltic as well as the Norwegian coast have been conducted only for reduced periods. Therefore, we combined recent and older stock assessments after consistency checks of SSB and R time-series, by simply replacing the newer part of the time-series of the older assessments with the newer time-series assessment (see Supplementary Figure S2).

We selected sea surface temperature (SST) and the indices of the North Atlantic oscillation (NAO) and Atlantic multidecadal oscillation (AMO) as climate predictors in our models. SST data were collected from the NOAA Extended Reconstructed Sea Surface Temperature dataset (ERSST, www.ncdc.noaa.gov) version 4. The dataset represents a reconstruction of SST from 1854 to the present and comprises monthly anomalies computed with respect to the period 1971-2000, resolved in a 2 Â 2 grid of spatial resolution. The data were averaged per year and per management unit. SST was chosen because of its importance for recruitment of Atlantic cod and is also a proxy for climate change at a local scale (Planque and Fre ´dou, 1999). NAO and AMO were used as indices of climate variability at the supraregional scale. In particular NAO has been shown to highly correlate with Atlantic cod recruitment (Stige et al., 2006), whereas AMO is a good proxy for climate change at longer timescales in this area. The NAO is a large-scale, high frequency (7-25 years) climatic index depending on the different atmospheric pressure at sea level between Iceland and Azores. The AMO is instead a large-scale, low frequency (60 years) multidecadal index representing climaterelated SST changes in the Atlantic Ocean. The data for both indices were collected from the Earth System Research Laboratory of NOAA (www.esrl.noaa.gov), and the AMO was averaged to annual values, whereas the NAO was averaged annually but just between December and March.

Modelling strategy

We compared multiple stock-recruitment models, the traditional Ricker model, the SCM, and MSPs (from the EDM framework). Recruitment models include either SSB alone or SSB in combination with one of the climate variables (i.e. SST, NAO, and AMO) as predictors. Because recruitment can be influenced by climatic factors at different life-stages (i.e. eggs, larvae, and juveniles), we applied multiple lags on the climate variables depending on recruitment age (Supplementary Table S1). We assessed the predictive power of the different models (three modelling approaches and explanatory variables and corresponding lags selection) on the test data using fivefold cross-validation, which randomly splits the time-series in five parts, fitting the model on four (training data), and using the results to predict the last one (test data). In each of the five iterations, we compared the predicted with the observed test values using Pearson correlation coefficients (q; Ye et al., 2015;Deyle et al., 2018). We repeated this procedure 100 times to increase the robustness and eventually used the median of the 500 values of q for model comparison (Figure 2).

The recruitment models

The Ricker Model fits a curve between recruitment and SSB depending on parameters a and b (Ricker, 1954). These parameters allow the curve to bend in the middle, so that at very high SSB values recruitment will be low because of density-dependent effects. However, this model is log-linear, i.e. the relationship between recruitment and biomass can be linearized through logtransformation, thus, we will refer to it as a linear model throughout the text. Climate effects can be added through a new parameter (c; Figure 2,box1):

where ageR is the age at recruitment, and lags the offset between the effect of a climate variable and R depending on the age of recruitment [i.e. for each stocks the climate variables were lagged from R t to R tÀageR depending on the age at recruitment (Supplementary Table S1)].

The starting values for the parameters were estimated from the linearized version of the function using the FSA (Ogle, 2016) package [Equation ( 2

(2)

Subsequently the Stock-Recruitment function was fitted to the data using a non-linear model with as response variable the logtransformed recruitment. SCM is based on the cusp, one of the seven canonical forms of catastrophe theory that describe sudden changes in a system because of small changes of external drivers (Thom, 1977;van der Maas et al., 2003;Petraitis and Dudgeon, 2016;Sguotti et al., 2019). The cusp model is based on a cubic differential Equation (3) and describes discontinuous transitions in a state variable z t , controlled by two control variables a and b, and thus can be used to describe discontinuous dynamics in recruitment (Figure 2, box 1).

where V z t ; a; b ð Þis a potential function representing the rate of change of the system (z t ), depending on the two control variables (a; bÞ:

Because natural processes and empirical data often include stochasticity, Equation (3) was reformulated as a stochastic differential equation, adding the Wiener process (r z dW t ) with variance r 2 :

where the first part of the equation is the drift term, r z is the diffusion parameter, and W t represents the Wiener process. The state variable, z t ; and the parameters, a and b [Equations (3) and (4)] are estimated as a linear function of one or more exogenous variables using a likelihood approach [Equations (5a)-(5c)]. In SCM the state variable R depends on two control variables (i.e. SSB and temperature). Although SSB controls the dimension of R, i.e. if R is found at the upper or lower shield, temperature controls the type of path that R will follow, either linear or discontinuous (shaped, i.e. the two red paths). The three-dimensional landscape can be projected on the two-dimensional plane in which the folded area, i.e. the area of discontinuous dynamics, is shown in grey. In the attractor reconstruction of R depending on SSB and climate made with multivariate simplex projection (MSP; right) every point in the attractor correspond to a time-step of the system. MSP allows for the prediction of the future step of the system based on Euclidean Distance Calculations and thus is a state-dependent approach. All methods were fit as baseline models using just SSB as control variable then adding the environmental variable. Box 2 shows the model evaluation procedure using fivefold cross-validation. From this procedure, the predictive power of all the models was calculated and finally compared (Box 3) to assess performance among models of the same type and between the three different methods.

where a 0 ; b 0 ; and z 0 are the intercepts and a 1 ; b 1 ; and z 1 the slopes of the models.

In our study z t , the state variable, is a linear function of recruitment. a is the so-called asymmetry parameter and controls the size of z t , thus in our study is a function of SSB. b is called the bifurcation parameter because it controls whether the state variable follows a continuous or discontinuous path (Petraitis and Dudgeon, 2015;Diks and Wang, 2016;Sguotti et al., 2019), and in our study is a function of the environmental variables [climate, see below Equations (7a)-( 7c)].

The system presents multiple equilibria if it follows a discontinuous path (i.e. two stable and one unstable) and just one if it follows a continuous path. The number of equilibria of the system depends on the solution of Equation ( 4), from which the Cardan's discriminant (d) is derived:

If d > 0, the system has one equilibrium, indicating a continuous path. Whereas if d 0 the system has three equilibria, indicating a discontinuous path (Diks and Wang, 2016). Therefore, SCM allows the detection of interactive effects of the two control variables on the state variable. Any changes in the bifurcation parameter b, will lead to changes in the relationship between a and z t and consequently dramatic changes of the state variable (Figure 2, box1, Supplementary Figure S3).

The model is represented by:

In order to test the predictive power of the model, we first produced the linear predictors of the parameters and the state variable. These were then fit into Equation (8) to predict the new points on the surface.

MSP is based on the EDM framework. The cornerstone of this framework is the Simplex Projection method. The principle of EDM is to reconstruct the dynamics of one or multiple time-series in a multidimensional space, i.e. an attractor manifold, and predict the future trajectory of the system based on these past dynamics (Figure 2, box1; Sugihara et al., 2012;Ye et al., 2015;Chang et al., 2017). Reconstructing the past dynamics of a system (in our case recruitment) is possible either using multiple variables (i.e. SSB or climate indices) or just time-lags of the system itself (i.e. recruitment; Sugihara et al., 2012). We here used differentiated recruitment time-series to build the attractor for each cod stock, and Simplex Projection [Equations ( 9) and ( 10)] to approximate the attractor dynamics of the system (Sugihara et al., 2012;Ye et al., 2015;Deyle et al., 2018). The time-series is transformed in a set of time-delayed coordinate vectors:

x ¼ fx t ; x tÀs ; x tÀ2s ; x tÀ3s ; . . . x tÀðEÀ1Þs g;

(9) where x is the system, in our study recruitment, t is time, s is the time-lag, and E the embedding dimension. E represents the dimensionality of the attractor (Ye et al., 2015). E is selected by predicting the attractor manifold one step ahead into the future (using leave-one-out cross-validation) then comparing the predictive power of models with a varying E. In order to predict the system into the future, x tþ1 , Euclidean distance is used and the system is predicted using nearest neighbourhood estimations

where w i represents the weights [Equation ( 11)], which are the Euclidean distance to the neighbour vector i relative to the nearest neighbour d .

MSP uses Equation ( 9) but with multiple variables. In our study, the attractor reconstruction of recruitment was based on SSB alone or together with climate variables [climate; Equations (12a) and (12b)]:

Performing MSP requires two preliminary tests, the S-Map and the convergent cross mapping (CCM), to unravel recruitment dynamics and the relationship between recruitment and explanatory variables, respectively.

EDM-specific preliminary tests S-Map and CCM

The S-Map, was performed after the attractor reconstruction with Simplex Projection. This test includes a tuning parameter h that controls the weights w i from Equation (10), and, if bigger than 0 indicates non-linearity (Sugihara, 1994;Klein et al., 2016;Dakos et al., 2017). Significance of non-linearity was assessed using a null distribution generated from 500 surrogate time-series for each S-Map model. The surrogate time-series were created following Deyle et al. (2018) and were phase-randomized, which preserves the basic statistical properties of the original time-series (Ebisuzaki, 1997). We averaged the S-Map results for all Atlantic cod recruitment time-series to understand the overall dynamics.

We performed CCM between R and SSB and the climate variables (SST, NAO, and AMO), a technique to understand causality between time-series without assuming any distribution (Sugihara et al., 2012;Deyle et al., 2016;Pierre et al., 2018). CCM is based on the principle that, if SSB or climate variables have an influence, then the R time-series will contain information about the past state of these variables. CCM is performed using Equation (10; see Deyle et al., 2018).

All analyses were performed in the programming environment R (R Core Team, 2016, version 3.3.1) using the packages FSA (Ogle, 2016), cusp (Grasman et al., 2009), rEDM (Ye et al., 2016).

Results

In our multi-model approach, we compared the parametric, linear Ricker model with two non-parametric, state-dependent approaches, i.e. the SCM and the state-dependent MSP, with or without environmental variables as predictors (Figure 2). The two preliminary tests of the EDM, necessary to perform the MSP, revealed on average significantly non-linear dynamics in recruitment of Atlantic cod stocks, and an appropriate choice of explanatory variables (Supplementary Figures S4 andS5), thus allowing us to proceed with the analyses. For most of the Atlantic cod stocks, the best performing models produced high correlations between observed and predicted values (0.7 < q < 0.8). An exception were northeast Arctic, Iceland, and Gulf of Maine cod stocks where the predictive power was lower compared to the other stocks (about q ¼ 0.4). Differences between the three model types were in general low (Figure 3). The Ricker model performed best for seven stocks, the SCM for eight stocks, and the MSP for five stocks (Figure 3, Supplementary Table S2). For stocks where SCM was the best, the MSP generally showed also a high predictive power, indicating that both models can well describe abrupt dynamics [e.g. Figure 3 (12,13,14,17)]. The addition of climate variables as explanatory variables to the baseline SSB models generally increased the predictive power, independently of the model type, although SSB was often the most correlated explanatory variable (Figure 3, Supplementary Table S2 and as shown in CCM, Supplementary Figure S5). SST and AMO were selected, based on the predictive power of the model, in respectively eight stocks and NAO in the remaining four stocks, generally agreeing with CCM results (Figure 4, Supplementary Figure S5). However, adding a climate variable had only a weak or even no additional effect when the baseline SSB model performed already poorly [e.g. Figure 3 (8), Ricker model].

The Ricker model best represented more gradual declines in recruitment, typical for cod stocks around the British Isles (i.e. North Sea, West of Scotland, and Irish Sea), those closer to the Arctic (i.e. Faroe Plateau, northeast Arctic, and Iceland cod) and Georges Bank cod [Figure 4 (4,5,6,10,9,11,19), Supplementary Table S3], as illustrated by their individual time-series (Supplementary Figure S1). All of these stocks, except Georges Bank, displayed strong density-dependence in recruitment (Supplementary Table S3), which is characteristic for the Ricker model. Furthermore, Ricker models clearly revealed that recruitment in warmer years is usually lower for the same level of SSB when compared to colder conditions (as indicated by low SST, NAO,or AMO in Figure 4 (4,5,6). The only exception with the reverse pattern of higher recruitment values at warmer conditions was northeast Arctic, hence the only cod stock that really profited from climate warming [Figure 4 (9)].

SCM instead is an approach from catastrophe theory, which models best discontinuous dynamics characterized by abrupt sudden shifts and hysteresis (i.e. in this case delayed recovery). The recruitment and SSB time-series of Canadian stocks on the western Atlantic side, but also Greenland and eastern Baltic cod (Supplementary Figure S1) show this type of dynamics, and hence SCM was the best approach for these stocks. SCM identified discontinuous stock-recruitment dynamics caused by the interaction of SSB and the climate variable. Moreover, SCM can identify catastrophic collapse, which occurs when SSB is found in the "folded" area, or area of instability (see blue shaded areas in Figure 4 (12,13,14,15,16,17,18,1), Supplementary Figure S3). Recruitment collapsed in these stocks, when in the instability area, in response to only small reductions in stock size 1)]. Consequently, SSB was a significant predictor in all SCMs, controlling the stocks dimension, while the climate variables modified the relationship between recruitment and SSB rendering it discontinuous, and thus inducing hysteresis (Supplementary Table S4). These two factors lead to the presence of stable low recruitment levels towards the end of the time-series. Low SSB coupled with warming (as indicated by climate variables SST, NAO, and AMO, Supplementary Table S4) had the potential to stabilize low recruitment. This is indicated by values outside the bifurcation area as best demonstrated by northern and Grand Banks cod [Figure 4 (13,17)]. Other cod stocks such as those from the Gulf of St Lawrence, on the eastern Scotian Shelf and off Greenland were at the boarder of stable low recruitment levels [Figure 4 (12,[14][15][16]].

Eventually, we found MSP to be the best model for recruitment of stocks that did not show collapses, but mostly fluctuating dynamics such as cod in the western Baltic, the Kattegat (because, even if the SCM was the best the model, the fit was invalid), the Celtic Sea, the Norwegian coast and in the Gulf of Maine [Figure 4 (3,2,7,8,20), Supplementary Figure S1]. The MSP however, seemed also appropriate to model catastrophic dynamics, but less effectively than the SCM. In contrast to the stocks best modelled with SCM and Ricker, stocks best modelled with MSP showed a mixed response to recent warming with a clear negative effect on recruitment in the western Baltic only [Figure 3 (2)].

Discussion

Short-term predictions of the size of an incoming year-class is essential to modern assessments of commercial fish species, but often suffers from the accuracy of available models predicting recruitment based on continuous, linear relationship with SSB. In our study, we investigated (i) whether recruitment dynamics in Atlantic cod stocks are better predicted by non-parametric, statedependent, or catastrophic statistical methodology compared to traditional parametric, linear approaches such as the Ricker stock-recruitment model and (ii) whether using climate variables as predictors in addition to SSB improves the predictive performance of such models.

The main result of our study is that predicting fish stock-recruitment can be improved by tailoring the modelling approach to the dynamical properties of each individual stock. We found cod stocks with more gradual and mostly linear dynamics to be best predicted by the traditional linear Ricker model, whereas stocks that experienced sudden abrupt changes in recruitment and stock size are best described by the SCM. SCM, based on catastrophe theory, is well suited to represent such discontinuous regime shift dynamics (Thom, 1972;Grasman et al., 2009;Diks and Wang, 2016;Sguotti et al., 2019). SCM allows for the identification of drivers and how their interactions result in unstable recruitment dynamics and hence provides a form of vulnerability assessment that can be instrumental in management (Petraitis and Dudgeon, 2015;Diks and Wang, 2016;Sguotti et al., 2019). Eventually, MSP was most appropriate for stocks that displayed The comparison between the predictive power of the best models resulting from the model selection between the Ricker, stochastic cusp model (SCM), and multivariate simplex projection (MSP) models. The median of the predictive power, derived from the cross-validation is shown for the three models without (blue, darker) and with (green, lighter) the inclusion of climate variables. The best model among the three, i.e. the model presenting the highest Pearson q between observed and predicted values of the test dataset, is indicated by a yellow (lighter) star for each stock. Black stars indicate the best models which however had a poor fit to recruitment and thus were substituted by the second-best model. The environmental variable that resulted in the best predictions can be found in Figure 4 and Supplementary Table S2. The number of years in the time-series is indicated for each stock. The colours underlying the names of the stocks correspond to the geographical location of the stock in pink (or from 1 to 11) in the East Atlantic, and orange (from 12 to 20) in the West. The numbers correspond to the stocks number in Figure1. (4,5,6,10,9,11,19). Results for cod stocks best represented by the Ricker model. The colour of the dots corresponds to the state of the climate variable, red (lighter) above the mean, blue (darker) below the mean. The two lines indicate separate SRRs for the two climate states (above and below the mean. (12-18, 1). Results for cod stocks best represented by the stochastic cusp model (SCM), dots are scaled to the size of R. The blue area (or grey in black and white) corresponds to the instability area, thus the fold in the three-dimensional visualization (Figure 2, Supplementary Figure S3) where three equilibria are possible. (3,2,7,8,20). Results for cod stocks best represented by multivariate simplex projection (MSP); dots are scaled to the size of SSB. Colours correspond to the state of the climate variable indicated, red (lighter) above the mean, blue below the mean. The lines show the predicted trends of R over time. The colours in the boxes correspond to the geographical location of the stock in pink (or from 1 to 11) in the East Atlantic, and orange (from 12 to 20) in the West. The numbers correspond to stock numbers in Figure 1. The comparison between observed and predicted values for each model can be seen in Supplementary Figure S6. more chaotic and fluctuating behaviours (Sugihara et al., 2012;Ye et al., 2015). Indeed, being a minimally assumptive model the most complex dynamics are better captured by it. MSP as part of the EDM suite of methods is based on attractor reconstruction and accounts for state-dependent dynamics (Ye et al., 2015), which makes it a suitable approach to model also discontinuous dynamics (Ye et al., 2015;Deyle et al., 2018). Mostly, both SCM and MSP models performed similarly in our analysis and their relatively high predictive power indicated the importance of using state-dependent and/or discontinuous approaches to model recruitment (Ye et al., 2015;Deyle et al., 2018;Munch et al., 2018).

Our study highlights that important differences exist between cod stocks in the eastern and western areas of the North Atlantic (Frank et al., 2016). Stocks from the western Atlantic and in particular off Canada and Greenland often experienced pronounced catastrophic dynamics, i.e. abrupt and sudden changes in stock size and recruitment. Eastern Atlantic stocks instead showed more continuous dynamics and thus a higher degree of stability. In general western Atlantic cod stocks seemed to be less resilient to abrupt collapses because of more fragile life-history traits, an overall more extreme and difficult environment, and different exploitation histories (Ra ¨tz and Lloret, 2003;Po ¨rtner et al., 2008;Wang et al., 2014;Frank et al., 2016). Moreover, SST was selected in eastern Atlantic cod stocks models, whereas for Western stocks the climate indices explained better the recruitment variability. This difference might indicate that the eastern cod stocks are more influenced by local processes, whereas in the western Atlantic large-scale climatic fluctuations are more important. Nevertheless, the addition of the climate factors in the best stockrecruitment models almost always increased its predictive power and thus highlights the importance of using environmental information also in stock assessment and management considerations to consider broader ecosystem dynamics (Punt et al., 2013;Skern-Mauritzen et al., 2016).

These results highlight the presence of multiple dynamics in cod stocks, which are also supported by the results of the preliminary S-Maps tests revealing a significant level of non-linearity in recruitment time-series of Atlantic cod stock. However, the nonlinearity signal is lower than expected, which we assume is because of the nature of the stock assessment data we used, and thus could be an underestimation (Brooks and Deroba, 2015). Such model output tends to be smoother and more linear than survey data (Storch et al., 2017), which are unfortunately not available for all cod stocks and longer periods needed for our study.

Finally, the different models allow us to draw conclusions about the recovery potential for collapsed Atlantic cod stocks. Most of the stocks are negatively influenced by warming and climate variability, because the lowest recruitment and SSB coincide with the highest temperature (Brander, 2005;Drinkwater, 2005;Po ¨rtner et al., 2008). The only exception is northeast Arctic cod where a warming environment positively influences recruitment, because this stocks resides at the northern distribution limit of the species (Po ¨rtner et al., 2008). Apart from northeast Arctic and Iceland cod where SSB has recently reached high levels, the stocks for which the traditional Ricker model performed best, such as the ones from the North Sea and around the British Isles, show continuously low recruitment and SSB in recent years and a continuous relationship between these parameters. These imply that, with low exploitation pressure these stocks have a higher recovery potential, but with climate change the productivity will likely remain low (Drinkwater, 2005). The situation is even worse for stocks that are best described by the SCM such as the western Atlantic stocks where the relationship between recruitment and SSB is discontinuous and thus the stocks display a strong hysteresis effect. Most of them are at present in a stable low state, suggesting that recovery might be even further delayed and productivity will remain low.

Conclusions

We demonstrated that discontinuous, state-dependent dynamics are pervasive in at least half of Atlantic cod stocks and need to be considered when predicting year-class strength. Indeed, even if our study does not necessarily reflect the goodness of the models to predict future recruitment, because the cross-validation included years after those predicted, we highlight the presence of different dynamics between stocks. Furthermore, we show the importance of accounting for environmental factors in recruitment predictions. Our findings indicate the need for more flexibility in the stock assessment process and highlight the importance for an adaptive multi-model approach that accounts for the inherent dynamics of living marine resource populations (Punt et al., 2016). Flexible models and adaptive management are fundamental to move towards an ecosystem-based management approach, especially in the face of climate change. To achieve this, we need to move away from fixed and established model procedures and explore options, to be ready to adapt to the new challenges that climate change will pose (King et al., 2015).