Abstract

Mantzouni, I., Sørensen, H., O'Hara, R. B., and MacKenzie, B. R. 2010. Hierarchical modelling of temperature and habitat size effects on population dynamics of North Atlantic cod. – ICES Journal of Marine Science, 67: 833–855.

Understanding how temperature affects cod (Gadus morhua) ecology is important for forecasting how populations will develop as climate changes in future. The effects of spawning-season temperature and habitat size on cod recruitment dynamics have been investigated across the North Atlantic. Ricker and Beverton and Holt stock–recruitment (SR) models were extended by applying hierarchical methods, mixed-effects models, and Bayesian inference to incorporate the influence of these ecosystem factors on model parameters representing cod maximum reproductive rate and carrying capacity. We identified the pattern of temperature effects on cod productivity at the species level and estimated SR model parameters with increased precision. Temperature impacts vary geographically, being positive in areas where temperatures are <5°C, and negative for higher temperatures. Using the relationship derived, it is possible to predict expected changes in population-specific reproductive rates and carrying capacities resulting from temperature increases. Further, carrying capacity covaries with available habitat size, explaining at least half its variability across stocks. These patterns improve our understanding of environmental impacts on key population parameters, which is required for an ecosystem approach to cod management, particularly under ocean-warming scenarios.

Introduction

Evaluation of stock status and recovery policies in fisheries management relies heavily on biological reference points estimated from the parameters of single-stock, stock–recruit (SR) models (Hilborn and Walters, 1992; Mace, 1994; FAO, 1995; Myers and Mertz, 1998a; Quinn and Deriso, 1999; Needle, 2002; ICES, 2003). In this context, the key parameters of SR models are the maximum reproductive rate at low stock size, and the maximum long-term equilibrium spawner biomass (i.e. the carrying capacity, K, of the specific ecosystem for the stock). Together, these parameters determine, for example, how quickly suppressed populations might recover and to what level. Further, Martell et al. (2008) showed recently that SR models can be parameterized in terms of key management quantities, e.g. maximum sustainable yield (MSY) and fishing mortality required to achieve MSY (FMSY). Therefore, it is now possible to allow for more transparent application of SR models in fisheries management, especially when informative priors for these reference points are used.

SR models are often characterized by non-stationarity, causing the parameters to vary within stocks (Walters, 1987; Myers and Mertz, 1998b; Needle, 2002). This variability is caused mainly by changes in the factors shaping the reproductive success of stocks: ecosystem conditions (e.g. trends in key biotic or abiotic variables, including physicochemical conditions, food availability, and abundance of predators and/or competitors) and fisheries exploitation (Walters, 1987; Hilborn and Walters, 1992; Kuparinen and Merilä, 2008; Rijnsdorp et al., 2009). These factors, also acting interactively (Perry et al., in press), directly affect survival probabilities, especially of the early life stages, and the reproductive output of the stock (through changes in age and size composition) and can also induce phenotypic or evolutionary changes in biological traits controlling productivity (e.g. behaviour, fecundity, growth rates, age and size at maturity, habitat selection). As a result, inferences based on the SR parameters, such as estimated recovery rates and target stock recovery levels, will have additional uncertainty and, in some cases, could give overly optimistic (or pessimistic) estimates of the resilience of populations to perturbations such as exploitation and ecosystem variability.

Usually, ecosystem conditions are not explicitly included in the derivation of SR models. As a consequence, estimated parameters are insensitive to non-random changes in ecosystem properties that affect population dynamics. This situation has at least two consequences. First, the ability of parameter values themselves to reflect true levels of future stock productivity deteriorates when new ecosystem conditions (e.g. foodweb configurations, temperatures) develop in the region occupied by the stock. In other words, parameter values may no longer be reliable under different ecosystem circumstances. There are many examples of situations where stock productivity has changed (for instance through regime shifts; ICES, 2007a), and there will be more in future (for instance, as climate change progresses). Second, variability of parameters (a measure of their uncertainty) is larger than if ecosystem factors that affect stock productivity were identifiable, quantifiable, and incorporated directly in SR models.

Here, we apply and develop analytical methods, with the aim to identify potential environmental patterns in recruitment dynamics and also to improve the estimation of biology-based SR model parameters across North Atlantic stocks of cod (Gadus morhua). First, given the nature of the cod multistock SR datasets, we use hierarchical modelling approaches to describe SR dynamics both within and across stocks. Hierarchical modelling is especially suitable in the present case, where multiparameter SR models are fitted to a wide range of cod stocks, which are expected to be related to some extent in their responses. Moreover, combining information across stocks allows borrowing strength for the estimation of individual parameters from the broader dataset (Ntzoufras, 2009). Therefore, these methods are advantageous, especially for stocks with shorter time-series, because they can reduce the uncertainty in estimates of SR model parameters (Myers and Mertz, 1998b; Myers, 2001). Second, we incorporate ecosystem properties in the analytical process to investigate whether accounting for such effects can provide an improved representation of SR dynamics. The hierarchical methodology implemented in this context allows one to quantify the magnitude and shape of the functional response of the parameters derived to variations in ecosystem properties at different ecological levels of organization. Most importantly, we attempt to maintain a balance between statistical and biological reasoning; SR models (Ricker and Beverton and Holt, BH) developed based on fish biological processes, rather than solely on the statistical goodness-of-fit, are employed. The approach therefore stands in the middle ground between assessment and management objectives, aiming to enhance understanding of natural processes or implement the precautionary approach, respectively.

The study has two related objectives. First, hierarchical approaches are employed to estimate SR model parameters and their uncertainty for all 21 cod stocks in the North Atlantic. Second, the outputs from the analysis are used to increase our understanding of how the productivity of individual cod stocks varies geographically, over time, and in relation to ecosystem properties. In the latter context, we focus on temperature because that variable has been shown to have important impacts on cod biology (e.g. Planque and Frédou, 1999; Brander, 2000; Drinkwater, 2005). Hierarchical approaches are used to quantify how local variations in temperature affect productivity and carrying capacity within a single stock, as well as the aggregate species-level functional response (both mean and uncertainty), across the range of cod stocks, as well as whether uncertainty in SR model parameters can be reduced by incorporating ecosystem properties.

Material and methods

Data

Our database consists of spawning-stock biomass (S), recruitment (R), and spring surface layer (0–100 m) temperature (T) time-series (presented in Figure S1), and habitat size (between 40 and 300 m deep) for each of 21 North Atlantic cod stocks (Table 1, Figure 1).

Figure 1.

Fisheries statistical areas in (a) the western, and (b, next page) the eastern North Atlantic. The locations of the cod stocks are listed in Table 1. The maps are reproduced with the permission of FAO.

Figure 1.

Fisheries statistical areas in (a) the western, and (b, next page) the eastern North Atlantic. The locations of the cod stocks are listed in Table 1. The maps are reproduced with the permission of FAO.

Table 1.

Summary of codes, geographic locations, time-series length, recruitment, and S (spawning-stock biomass) data sources and SPRF=0 estimates for the stocks used in the analyses.

Stock Area SPRF=0 (kg) Habitat size (km2Year range Reference for stock data 
codgb Georges Bank (5Z) 23.8a 85 306 1978–2004 O'Brien et al. (2005) 
codgom Gulf of Maine (5Y) 27.9a 47 629 1982–2003 Mayo and Col (2005) 
cod4x Western Scotian Shelf (4X) 14.7a 64 331 1983–2000 Clark et al. (2002) 
cod4vsw Eastern Scotian Shelf (4VsW) 11.7a 91 656 1970–2001 Fanning et al. (2003) 
cod4tvn* Southern Gulf of St Lawrence (4TVn) 7.0a 80 338 1950–2004 Chouinard et al. (2006) 
cod3pn4rs* Northern Gulf of St Lawrence (3Pn4RS) 4.1a 89 734 1974–2003 Fréchet et al. (2005) 
cod3ps* Southern Newfoundland (3Ps) 5.9a 55 509 1977–2002 Brattey et al. (2004) 
cod3no* Grand Banks (3NO) 6.3a 114 967 1959–2004 Power et al. (2005) 
cod3m Flemish Cap (3M) 9.8a 17 398 1972–1993 Vázquez and Cerviño (2002) 
cod2j3kl* Northern Newfoundland (2J3KL) 7.4a 263 893 1962–1989 Bishop et al. (1993) 
cod-iceg* Iceland (Va) 18.9b 63 158 1955–2003 ICES DB 2006 
cod-farp* Faroe Plateau (Vb) 11.6b 12 883 1961–2004 ICES DB 2006 
cod-arct* Northeast Arctic (I, II) 12.1b 858 979 1946–2003 ICES DB 2006 
cod-coas* Norwegian Coastal (IIa) 6.2c 234 615 1984–2004 ICES DB 2006 
cod-2224* Western Baltic (IIId—west) 5.3b 41 823 1970–2005 ICES DB 2006 
cod-2532* Eastern Baltic (IIId—east) 3.3c 131 019 1966–2003 ICES DB 2006 
cod-kat Kattegat (IIIa—east) 7.8b 19 926 1971–2004 ICES DB 2006 
cod-347d North Sea (IIIa, IV, VIId) 18.2b 418 821 1963–2004 ICES (2007b) 
codviia Irish Sea (VIIa) 12.7b 23 940 1968–2005 ICES (2006) 
cod-7ek Celtic Sea (VIIe–k) 19.9b 177 620 1971–2005 ICES DB 2006 
codvia West of Scotland (VIa) 12.9b 81 016 1978–2004 ICES (2006) 
Stock Area SPRF=0 (kg) Habitat size (km2Year range Reference for stock data 
codgb Georges Bank (5Z) 23.8a 85 306 1978–2004 O'Brien et al. (2005) 
codgom Gulf of Maine (5Y) 27.9a 47 629 1982–2003 Mayo and Col (2005) 
cod4x Western Scotian Shelf (4X) 14.7a 64 331 1983–2000 Clark et al. (2002) 
cod4vsw Eastern Scotian Shelf (4VsW) 11.7a 91 656 1970–2001 Fanning et al. (2003) 
cod4tvn* Southern Gulf of St Lawrence (4TVn) 7.0a 80 338 1950–2004 Chouinard et al. (2006) 
cod3pn4rs* Northern Gulf of St Lawrence (3Pn4RS) 4.1a 89 734 1974–2003 Fréchet et al. (2005) 
cod3ps* Southern Newfoundland (3Ps) 5.9a 55 509 1977–2002 Brattey et al. (2004) 
cod3no* Grand Banks (3NO) 6.3a 114 967 1959–2004 Power et al. (2005) 
cod3m Flemish Cap (3M) 9.8a 17 398 1972–1993 Vázquez and Cerviño (2002) 
cod2j3kl* Northern Newfoundland (2J3KL) 7.4a 263 893 1962–1989 Bishop et al. (1993) 
cod-iceg* Iceland (Va) 18.9b 63 158 1955–2003 ICES DB 2006 
cod-farp* Faroe Plateau (Vb) 11.6b 12 883 1961–2004 ICES DB 2006 
cod-arct* Northeast Arctic (I, II) 12.1b 858 979 1946–2003 ICES DB 2006 
cod-coas* Norwegian Coastal (IIa) 6.2c 234 615 1984–2004 ICES DB 2006 
cod-2224* Western Baltic (IIId—west) 5.3b 41 823 1970–2005 ICES DB 2006 
cod-2532* Eastern Baltic (IIId—east) 3.3c 131 019 1966–2003 ICES DB 2006 
cod-kat Kattegat (IIIa—east) 7.8b 19 926 1971–2004 ICES DB 2006 
cod-347d North Sea (IIIa, IV, VIId) 18.2b 418 821 1963–2004 ICES (2007b) 
codviia Irish Sea (VIIa) 12.7b 23 940 1968–2005 ICES (2006) 
cod-7ek Celtic Sea (VIIe–k) 19.9b 177 620 1971–2005 ICES DB 2006 
codvia West of Scotland (VIa) 12.9b 81 016 1978–2004 ICES (2006) 

Data for most Northeast Atlantic stocks were extracted from the 2006 ICES Stock Assessment Summary Database (ICES DB 2006; http://www.ices.dk/datacentre/StdGraphDB.asp), unless otherwise stated. The stocks marked with an asterisk display significant autocorrelation at lag 1 in the log(R/S) time-series. SPRF =0 refers to spawners produced per recruit (SPR) in the absence of fishing mortality (F). The areas are shown in Figure 1.

cEstimated by the authors.

Population data for the eastern and western stocks were extracted from ICES and NAFO (Northwest Atlantic Fisheries Organization) published reports, respectively (Table 1). The estimates are based on sequential population analyses, standardized usually with fisheries-independent (such as research trawl survey) data.

Cod prerecruits (eggs, larvae, and pelagic juveniles; approximately ages 0 to ca. 3–4 months) are the most sensitive to temperature, so it is during that period that possible effects can emerge (Brander, 2000). We assumed that the upper water layer between 0 and 100 m was representative of the pelagic environment experienced by these stages, so spring (March–May) temperature data in the 0–100 m water layer were used. The time-series for the eastern North Atlantic were extracted from the ICES oceanographic database (http://www.ices.dk/datacentre/), and data for the western areas were compiled using the DFO (Fisheries and Oceans Canada) oceanographic databases (Gregory, 2004). Note, however, that special considerations apply to certain areas; for the eastern Baltic stock (ICES Subdivisions 25–32), temperature estimates were obtained by averaging over ICES Subdivisions 25–29, because the low salinity in Subdivisions 30–32 is unfavourable for cod reproduction (Nissling and Westin, 1991). For the Barents Sea (arct), temperature was estimated for the area south of 78°N because cold water can limit cod distribution (Ottersen et al., 1998). For similar reasons, temperature for the Icelandic cod (iceg) corresponds to the region south of 62°N in Division Va (ICES, 2005). Spatial restrictions were also applied in four Northwest Atlantic areas [Flemish Cap (NAFO Division 3M), Grand Bank (3NO), and western and eastern Scotian Shelf (4VsW and 4X)], which extend offshore into areas that could be affected by the Gulf Stream; for those areas, temperature observations south of 42°N were excluded. For area 3M, temperature data from the Flemish Cap only were used.

Estimates of juvenile habitat size (40–300 m) were also compiled for every area because available space at that stage is assumed to limit cod recruitment (Myers and Cadigan, 1993). A similar index is associated with a significant proportion (36%) of the long-term mean recruitment of 20 cod stocks in the North Atlantic (MacKenzie et al., 2003). It is assumed that habitat area is a reasonable proxy for the large number of processes and mechanisms that influence juvenile cod survival, but which remain poorly quantified in most areas of the North Atlantic. Habitat size was estimated by calculating the area between 40 and 300 m in the fisheries division occupied by the stocks (Figure 1). Note that the same spatial considerations regarding the exclusion of certain subareas for the estimation of the T time-series were also taken into account. Therefore, by available habitat size, we define the continental shelf area between the 40- and 300-m bathymetric contours within the entire statistical, ICES or NAFO, division(s) occupied by each stock, but excluding certain subareas, as described in the paragraph above.

To combine data on recruitment across cod stocks, the time-series have to be standardized for differences in life-history characteristics affecting reproductive output. For example, recruitment in some stocks is estimated at age 1, but at age 3 in others. Standardization was performed following the method used in previous meta-analyses of cod recruitment (e.g. Myers et al., 1999, 2001). According to this method, recruit numbers are multiplied by the stock-specific SPRF=0 parameter (spawners produced per recruit in the absence of fishing mortality), estimated as a function of natural mortality, weight-at-age, and age-at-maturity (Mace, 1994; Myers and Mertz, 1998b). In this way, recruits are transformed to total spawner biomass produced per spawner over its lifetime in the absence of fishing, and this also affects the interpretation of the SR model parameters, as discussed later. For most stocks, estimates of the SPRF=0 parameter are provided by Goodwin et al. (2006) and Shelton et al. (2006); for the others, they were estimated using data extracted from the assessment reports (Table 1).

Modelling and statistical analyses

Hierarchical models

Hierarchical or multilevel modelling is a rigorous probabilistic framework appropriate for combining data and making inference across various independent sources with similar characteristics and expected to exhibit comparable patterns in their dynamics (Berliner, 1996; Hilborn and Liermann, 1998; Wikle, 2003; Gelman et al., 2004; Gelman and Hill, 2007). The method makes use of these similarities to define a common population distribution, from which the unit-specific (here, stock) parameters are sampled (Gelman et al., 2004). The common distribution not only allows an improved fit of the model to the data, but also structures dependence among parameters, thereby avoiding overfitting (Gelman et al., 2004). In this context, the exchangeability assumption plays a key role (Gelman et al., 2004; Ntzoufras, 2009). The assumption is used when no relevant information, apart from the modelled data, is available to distinguish among the stock-specific parameters, so no ordering or grouping can be applied (Gelman et al., 2004). As will be shown later, however, it is possible to relax the exchangeability assumption by introducing stock-specific characteristics, besides the SR data, potentially influencing the parameters (Su et al., 2004). An additional advantage of hierarchical inference is that, especially when combined with state–space modelling methods, it allows for the explicit incorporation, and hence isolation, of the observation and systematic model error (Meyer and Millar, 2001), usually characterizing fisheries models (Hilborn and Walters, 1992; MacKenzie et al., 2003). This point is further considered in the Discussion section.

Implementation of the models is based on decomposing them into three stages or levels (Berliner, 1996; Wikle, 2003; Clark, 2007; McCarthy, 2007). The first level describes the probability of the data, given the explanatory variables and the parameters describing the corresponding effects (i.e. the functional form of the SR model). The second level describes the variation in the parameters that determine the dynamics in the first level. Its importance is twofold because at that level we define (i) the functional form for how ecosystem factors affect the parameters, and (ii) the distribution of SR model parameters across cod stocks. The latter can be extended to account for mechanisms generating among-stocks differences, so this stage is also referred to as the stock-level model. The third level is the parameter model, and it concerns the hyperparameters used to define the probability distributions of the parameters in the previous stages. These last two levels are based on the assumption that certain SR model parameters are connected across stocks and also on the exchangeability assumption; hence, they lie in the core of the hierarchical inference. The common probability distribution and the process generating the parameters, or describing the differences among them, are both described by the hyperparameters of the third stage and, therefore, form the interface for the combination of the individual datasets and, consequently, for exchange of estimation strength across stocks (Gelman et al., 2004).

Owing to their probabilistic background, most hierarchical applications have been implemented under the Bayesian paradigm (Gelman et al., 2004; Clark, 2007). Hierarchical Bayesian inference averages over the above levels of uncertainty and variability by likelihood (data model), priors (process or stock-level models), and hyperpriors (parameter model) to produce the posterior distribution of the SR model parameters (see also Supplementary material). The posterior distributions are interpreted as descriptions of the uncertainty about the parameters. Another popular framework in this context is frequentist mixed (or variance components) modelling (Searle et al., 1992; Snijders and Bosker, 1999; Demidenko, 2004; West et al., 2006; Clark, 2007; Gelman and Hill, 2007). In particular, the estimated variance components describe variability across stocks. Mixed models can be regarded as a combination of Bayesian and frequentist approaches (Demidenko, 2004), with results parallel to empirical Bayesian inference (Robinson, 1991; Snijders and Bosker, 1999). Both Bayesian and frequentist mixed models treat some parameters as random variables, with common (estimated) variances. For frequentists, the hyperparameters and additional unknown parameters used to specify the variance of the random effects and other parts of the model are considered fixed, but unknown, and are estimated from the data by maximizing their likelihood (or a similar measure). The distributions of the estimates describe sampling variability (uncertainty of estimates under resampling). In the Bayesian inference, however, these parameters are given hyper-prior distributions, and the full (posterior) distribution of the parameter is estimated (Demidenko, 2004).

Both approaches have some advantages for implementing multilevel modelling (Clark, 2007; Gelman and Hill, 2007). Mixed models are usually quick and easy to fit, but estimation may fail under certain circumstances. For instance, the non-linear, mixed-models approach does not have a closed-form solution, leading to more computationally intensive estimation algorithms and to less reliable inference results (Pinheiro and Bates, 2000). Bayesian hierarchical models, on the other hand, are more flexible, allowing estimation for more complex model structures, as well as easier inference about the uncertainty of variance components.

In this study, the hierarchical SR models were developed mainly in the Bayesian framework and simulated in BUGS (Lunn et al., 2000). Mixed modelling approaches, implemented in the R language and environment for statistical computing using the nlme library (Pinheiro and Bates, 2000), were also employed in a complementary manner to identify the best model structure for the linear SR models. The parameter estimates provided by the two frameworks were compared and, as will be shown, there was sufficient comparability, justifying the choice of combining both methods. Subsequently, the more flexible Bayesian framework was used to implement complex, non-linear, SR model formulations, to estimate inferential uncertainty about all model parameters, and to present key results. Both mixed modelling and Bayesian terminology were used to describe the SR model development below to provide a more complete overview of the principles and methods of hierarchical models.

SR models

The Ricker SR model (Ricker, 1954) is one of the standard models used in fisheries science (Hilborn and Walters, 1992):  

1
formula
where R is recruitment, standardized as described previously, S the spawner biomass, t denotes year for S and year class for R, and the process errors εt are assumed to be lognormal and multiplicative. Parameter ARIC represents the slope of the curve near the origin and is related to stock productivity and the density-independent survival rate. In the present case, where standardized recruitment is used, it can be interpreted as the average rate at which replacement spawners are produced per spawner over its lifetime at low spawner abundance and in the absence of fishing mortality. This rate is standardized across stocks for differences in weight-, maturity-, and natural mortality-at-age that are incorporated in the standardization parameter SPRF=0. Therefore, as will be shown in detail below, ARIC can depend on a number of time-varying factors, such as temperature, accounting for effects on prerecruit survival rates and also for potential fluctuations in the SPRF=0 components.

Parameter B is related to carrying capacity because 1/B) equals spawner biomass when recruitment reaches the maximum (Smax, also termed beta), and −B represents the density (stock)-dependent mortality dominating after this point. The parameter therefore depends on habitat size, which differs across regions and can cause, and explain, at least some of the across-stocks variability in beta. Moreover, the suitable habitat available can also be influenced by ecosystem variables and vary in time within each stock (Kell et al., 2005a). The possibility for within-stock temperature effects on beta, as well as the among-stocks relationship between beta and habitat size, is investigated below.

The K indices for a given stock i can be derived using the following formulation of the Ricker SR model, expressed in terms of the curve's maximum point Rmax and Smax (Brander and Mohn, 2004a, b): forumla. According to that formula, the parameters of the Ricker model can be expressed as ARIC= exp(1)(Rmax/Smax) and beta = 1/B = Smax. Therefore, Kmax, representing the maximum number of recruits produced by the maximum number of spawners sustained by the ecosystem (i.e. Rmax), is estimated as forumla. Keq is the carrying capacity when S and R are at equilibrium (i.e. when R = S) and hence is a useful parameter for management, representing the minimum spawner biomass required to produce replacement recruitment: forumla. The Ricker model can be linearized by natural log-transformation and assuming lognormal errors: log(Rt/St) = log(ARIC) − BSt + εt. To simplify notation, the Ricker model for stock i is written as  

2
formula
where yit = log(Rit/Sit), forumlaforumla, and i denotes the stock. For simplicity, we denote forumla as alpha, understanding that alpha = log(ARIC). The errors in this and the following models are assumed to be stock-specific.

The BH (Beverton and Holt, 1957) model is also broadly used for the study of SR dynamics:  

3
formula
Ricker and BH SR models display similar behaviour at the limit of low S (i.e. compensation); therefore, parameter A (denoted as ABH for the BH model) has common interpretation and estimation in both models (Myers et al., 1999). Parameter KBH has the same dimensions as S and can be interpreted as the “threshold biomass” resulting in half the maximum recruitment, which equals ABH × KBH.

The BH model cannot be linearized and, as discussed above, was developed in the Bayesian framework. In this context, a useful reformulation is the following:  

4
formula
where yit = log(Rit/Sit), xit = Sit, forumla and i denotes the stock. In this form, the model is linear to forumla, which has the same interpretation as forumla of the Ricker model in Equation (2) and will also be denoted as alpha, and forumla is comparable with the Ricker model forumla, because it corresponds to S, resulting in half the maximum recruitment given by forumla. For simplicity, forumla is also termed beta (distinction between the Ricker and BH model parameters is made on the basis of context, or reference is made to both parameters jointly). We can also estimate K at equilibrium, Keq, for the BH model as forumla. Recruitment K indices under the two SR models are estimated as a function of alpha and beta, so depend on both density-independent and -dependent processes represented by each parameter.

The SR models were, subsequently, developed in the “multilevel” framework. The Ricker model is used to present this section, simplifying notation by dropping the SR model identifiers, i.e. using αi and βi in the place of forumla and forumla, respectively. However, the same approaches apply to the BH model and its parameters forumla and forumla. Estimation for this model is described in the Supplementary material.

A convenient way to conceptualize the multilevel SR framework is by starting with a simple regression model fitted to all stocks (Gelman and Hill, 2007). This type of model is referred to as a complete-pooling model and is based on the “extreme” assumption that each parameter is fixed to a certain value, common across stocks. For the Ricker model, the complete-pooling model can be written simply as yit = α + βxit + εit, with forumla

The previous model can be written in another generalized way, as in Equation (2). In this form, it corresponds to the no-pooling model, a classic regression model that can be estimated for each stock separately, using indicators and assuming that parameters are completely independent across stocks. In the Bayesian framework, the data-level models represent the likelihood, describing the distribution of the data given the model. In other words, they convey information about the range of parameter values that are most consistent (likely) with the data for each stock (Gelman and Hill, 2007).

Hierarchical modelling is a compromise between these two extremes, imposing a soft constraint on the stock-specific parameters by assuming that they are derived from a common probability distribution (Gelman and Hill, 2007). The lowest-level model is extended by introducing the next level of complexity, i.e. specifying the distribution of the data model (i.e. the Ricker SR model) parameters across stocks. These across-stock distributions of the parameters represent the stock-level models and are estimated from the full dataset, so strength is borrowed across stocks. Usually, Gaussian distributions are assumed:  

5a
formula
and  
5b
formula
where μα and forumla (or μβ and forumla) are the mean and variance of the parameter alpha (or β = −1/beta) distribution, respectively. The distribution means, μα and μβ, are referred to as fixed effects in mixed-models terminology, and they represent the average value of the corresponding parameter across all stocks, whereas forumla and forumla are the variance components (Searle et al., 1992; Pinheiro and Bates, 2000). The previous stock-level models [Equations (5a) and (5b)] can also be written as:  
6a
formula
and  
6b
formula
where forumla and forumla are the stock-level model errors referred to as random effects and represent the deviation of stock i parameter αi or βi from the corresponding across-stocks mean. Jointly, they are represented by a multivariate normal distribution:

 

formula
where ρ is the correlation between the random effects. In the Bayesian context, the stock-level models represent the priors, which usually convey the existing (i.e. before the data under study are seen) knowledge of the parameters and are hence used to update or weight the likelihood. In the present context, however, whereby priors are common to all stocks, they are used to incorporate information on the distribution of parameters across stocks, which is relevant to obtaining individual estimates.

Consequently, the Ricker data-level model incorporating the stock-level models is  

7
formula
It is assumed that the residuals forumla associated with each stock are independent and that the residuals and random effects are independent of each other. The parameter forumla is also a variance component (Searle et al., 1992).

The stock-level models partially pool the parameters towards the mean of the distribution, so the estimates are referred to as shrinkage estimates (Gelman and Hill, 2007). For instance, the estimate of αi can be expressed as a weighted average of the no-pooling, stock-specific model forumla, corresponding to forumla in Equation (2) and its mean across stocks (the fixed effect) μα:

 

8
formula
where forumla is the stock-specific model variance and ni the number of observations for stock i (Gelman and Hill, 2007). Pooling is stronger when forumla is small and for stocks with fewer observations and greater variability (forumla). The stock-specific weight forumla is defined as the reliability of the corresponding stock mean, and the ratio of the two weights is the ratio of the true variance forumla to the error variance forumla. As the true values are unknown, they are substituted by their estimates in Equation (8). The above estimate forumla is also known as the posterior mean or the empirical Bayes estimate, because it is obtained by combining the stock with the population (representing the prior, in this context) information (Snijders and Bosker, 1999).

Incorporating ecosystem effects

Having introduced the hierarchical model as composed of uncertainty/variation levels [within (data-level) and across (process- and parameter-level) stocks], it is straightforward to extend the previous model to incorporate ecosystem (log-transformed habitat size, H, and temperature, T) effects on the parameters. The stock-level models [Equations (6a) and (6b)], describing the across-stocks variation in the parameters, are simply modified to account for the effects of these predictors. Therefore, stock-specific parameters are standardized for differences in known characteristics of the ecosystems occupied by cod stocks.

Temperature

Initially, both parameters alpha and beta are allowed to be temperature-dependent, so the stock-specific SR parameters are now also time-varying. Different functional forms of these relationships can be assumed and tested. To facilitate model comparison using the linear, mixed modelling framework, the relationships are initially assumed to be quadratic to allow for dome-shaped impacts, whereby the size and the magnitude of the effect depend on the temperature range. Therefore,  

9a
formula
and  
9b
formula

Those relationships, however, can impose structural constraints on the effect because they imply equal rates of increase and decrease in the parameters with temperature beyond the optimum point. Therefore, more flexible relationships, such as the inverse polynomials (Nelder, 1966), were also tested using the Bayesian models. For example, the dependence of parameter alpha on temperature can be written as:  

9c
formula

Note that the temperature time-series in Equations (9a) and (9b) were centred to the overall (across stocks) mean (i.e. average of all Tit observations) to remove the correlation between the first- and the second-order term estimates. Therefore, the intercepts coi and doi in the above relationships represent the stock-specific values of alpha and beta at across-stocks mean temperature. These among-stocks differences remain after accounting for temperature effects and can be attributed to additional ecosystem factors affecting the SR-model parameters in a relatively stable manner. The intercepts can be modelled by assuming across-stocks distributions (stock-level models):  

10a
formula
and  
10b
formula
These models describe the (constant in time) divergence between the stock-specific alpha or beta estimate and the across-stocks grand mean forumla or forumla, respectively. The coefficients describing the temperature effect on alpha and beta, (forumla, forumla) and (forumla, forumla), respectively, are assumed to be common across stocks. This assumption can be relaxed by imposing stock-level models on the coefficients. For instance, forumla in Equation (9a) can be substituted by the stock-specific terms forumla. Consequently, the Ricker SR model, updated to incorporate temperature effects on alpha and beta, is  
11
formula
where αit and βit are given by Equations (9a) [or (9c)] and (9b), respectively.

Habitat size

As discussed previously, parameter β = −1/beta of the Ricker model is a proxy of the carrying capacity of the ecosystem for stock i. K is expected to be positively related to the size of the stock-specific habitat, defined as the juvenile feeding ground, where space can be a limiting factor causing density-dependent effects (Myers et al., 2001; Myers, 2002). Initially, a quadratic relationship is assumed to allow for the possibility of non-linear effects. Differences in habitat size can, therefore, explain part of the across-stocks variability in β = −1/beta, which is represented by the distribution of intercepts, doi, in the stock-level model [Equation (10b)]. Therefore, H (the natural log of the habitat area, centred to the mean to eliminate correlation between the coefficients) is included as a predictor in Equation (10b) and the β stock-level model becomes  

12
formula

Combining Equations (9a), (9b), (10a), (11), and (12), the full model is  

13
formula
with random effects forumlaforumlaforumlaforumlaforumla and forumla Also, note that Equation (9c) can be used instead of Equation (9a).

The model can be conceptualized in the following ways (Gelman and Hill, 2007). The first way is with a two-level model, with the first level describing the dependence of yi on xi for each stock i, through the Ricker SR model updated to include temperature-related processes in its functional form, as in Equation (11). The second level includes the across-stocks distributions of the first-level model parameters [i.e. stock-level models in Equations (10a) and (12)]. In the second way, in the Bayesian inference framework, the second-level model [Equation (12)] imposes different prior distributions for the βi values, which are now non-exchangeable, because they depend on the stock-specific Hi.

Identification of best model structure

For the mixed models, likelihood ratio tests (LRT) were used to compare the different formulations fitted using REML (restricted maximum likelihood) or ML (maximum likelihood), depending on whether the test applied to random or fixed components, respectively (Pinheiro and Bates, 2000). The LRT results generally agree with the AIC (Akaike Information Criterion) comparisons.

The Deviance Information Criterion (DIC; Spiegelhalter et al., 2002), a generalization of the AIC, was used to compare the different Bayesian model formulations. The DIC is estimated as DIC = mean deviance + 2pd. The mean deviance is estimated as −2 × the log-likelihood averaged over the number of simulations, hence quantifying the lack of model fit (smaller is better). The addition of pd, the effective number of parameters, is a penalty for complex models. In other words, using the DIC criterion, we seek models with good technical fit (high likelihood) that are not overly complex or overfitted.

Assessments of model suitability

The uncertainty associated with the estimated model parameters was compared between the hierarchical Bayesian Ricker (HBR) and the single-stock Ricker (SSR) models to assess the degree of “borrowing strength” achieved with the former models. The HBR goodness-of-fit was also evaluated in terms of the coherence between the trends in the fitted and observed data.

Comparison of uncertainty in SR model parameters between HBR and SSR

The precision in the estimates of alpha and beta = −1/β were compared between the hierarchical and the single-stock Ricker models to demonstrate the reduction in uncertainty that could be achieved by applying a hierarchical approach incorporating ecosystem properties (e.g. temperature and habitat size effects). As the parameter estimates differ in both their means (μ) and standard errors (s.e.), the latter were expressed in units of the mean by using the s.e./μ ratio, here defined as an estimate of uncertainty indicator (i.e. analogous to the more familiar coefficient of variability). For alpha, the SSR estimate corresponds to an estimate under the long-term, stock-specific temperature conditions experienced by the stock during the period for which data were available. Therefore, when comparing the values of alpha from the HBR and SSR models, the hierarchical alpha values corresponding to stock-specific mean temperatures were used. The comparison itself was conducted by calculating the ratio of the single-stock to hierarchical-uncertainty indicators. Values >1, therefore, quantify the reductions in uncertainty achieved using hierarchical approaches.

Evaluation of pattern coherence between data and model predictions

The ability of the Ricker-Bayesian model to reproduce the main trends and fluctuations in the observed recruitment survival over time was evaluated by plotting the observed SR data and comparing the patterns with the predictions of the hierarchical SR models. Also, the coherence between recruitment survival and productivity patterns with temperature was assessed.

Results

Ricker SR model formulation

Both mixed modelling and Bayesian approaches were used to identify the final structure of the hierarchical linear Ricker SR model. Therefore, the two methods, based on different principles, were first checked for consistency in the parameter estimates, and the results were satisfactory (Figure S2).

The simplest model form [Equation (7)] was used as the basis to build the final model by identifying (i) whether the random effects associated with alpha and beta should be assumed independent or correlated, (ii) whether the residual variance should be assumed common or stock-specific, and (iii) the significance and the patterns of the temperature and the habitat-size effects. In particular, we were interested in exploring the possible temperature effect, so this effect was tested as the final one, where other non-significant effects were eliminated. The mixed-Ricker (MR) model codes and structure and the AIC, statistics, and p-values of LRT are summarized in Table 2.

Table 2.

Summary of Ricker mixed and Bayesian models developed and compared in analyses of temperature, habitat size, and spawner-biomass effects on cod recruitment in 21 populations in the North Atlantic (see Table 1 and Figure 1 for locations of populations).

Model code Model structure
 
AIC (fitting method) logLik (fitting method) [d.f.] Hypothesis tested (model comparison) LRT p-value (L. ratio) 
Fixed effects Random effects 
Mixed models 
 MR.1 forumla forumla 1 662.5 (REML) −825.2 (REML) [6] – – 
 MR.2 As in MR.1 forumla 1 661.6 (REML); −825.8 (REML) [5] Random effects are independent; ρ = 0 (MR.1 vs. MR.20.30 (1.1) 
 MR.3 forumla As in MR.2 1 595.5 (REML); 1 578.9 (ML) −772.7 (REML); −764.4 (ML) [25] Residual errors are stock-specific (MR.3 vs. MR.2) <0.01 (106.2) 
 MR3.HL forumla forumla, forumla 1 576.3 (ML) −762.2 (ML) [26] Linear dependence of βi on habitat size (H); fH1 = 0 (MR3.HL vs. MR.3) 0.03 (4.6) 
 MR3.HQ forumla As in MR3.HL 1 561.2 (ML) −753.6 (ML) [27] Quadratic dependence of βi on habitat size (H); fH2 = 0 (MR3.HQ vs. MR3.HL) <0.01 (17.2) 
 MR3.HQ.TAL forumla forumla, forumla1 562.4 (ML) −753.2(ML) [28] Linear dependence of αi on temperature (T); forumla (MR3.HQ.TAL vs. MR3.HQ.TAQ<0.01 (14.3) 
 MR3.HQ.TAQ forumla as in MR3.HL.TAL 1 550.2 (ML); 1 606.8 (REML) −746.1 (ML); −774.4 (REML) [29] Quadratic dependence of αi on temperature (T); forumla (MR3.HQ.TAQ vs. MR3.HQ) <0.01 (15.0) 
 MR3.HQ.TBL forumla forumla, forumla 1 562.6 (ML) −753.3 (ML) [28] Linear dependence of βi on temperature (T); forumla (MR3.HQ.TBL vs. MR3.HQ0.46 (0.6) 
 MR3.HQ.TBQ forumla As in MR3.HQ.TBL 1 562.2 (ML) −752.3 (ML) [29] Quadratic dependence of βi on temperature (T); forumla (MR3.HQ.TBQ vs. MR3.HQ0.23 (2.9) 
 MR3.HQ.TAQ.RS forumla forumla, forumla,forumla, forumla 1 609.4 (REML) −773.7 (REML) [31] Across-stocks differences in the influence of T on alpha; forumla (MR3.HQ.TAQ.RS vs. MR3.HQ.TAQ0.5 (1.4) 
 MR3.HQ.AC forumla As in MR3.HQ 1 398.3 (ML) −651.2 (ML) [48] – – 
 MR3.HQ.TAQ.AC forumla As in MR3.HQ.TAQ 1 386.6 (ML) −643.3 (ML) [50] Quadratic dependence of αi on temperature (T) after first-order differencing of yit; forumla (MR3.HQ.TAQ.AC vs. MR3.HQ.AC) <0.01 (15.7) 
 MR3.HL.r As in MR3.HL, but fitted only to stocks with low autocorrelation in yit 665.3 (ML) −317.7 (ML) [15] – – 
 MR3.HL.TAL.r forumla As in MR3.HQ.TAQ 660.7 (ML) −314.4 (ML) [16] Linear dependence of αi on temperature (T); forumla (MR3.HL.TAL.r vs. MR3.HL.r) 0.01 (6.6) 

 
  DIC ΔDIC   

 
Bayesian models 
 BR3.HQ As in MR3.HQ 1 508    
 BR3.HQ.TAQ As in MR3.HQ.TAQ 1 489 −19   
 BR3.HQ.TAQ.RS As in MR3.HQ.TAQ.RS 1 480 −9   
 BR3.HQ.TAN.RS As in MR3.HQ.TAQ.RS, but using the inverse polynomial [Equation (9c)] to describe the dependence of alpha on T 1 500 +20   
Model code Model structure
 
AIC (fitting method) logLik (fitting method) [d.f.] Hypothesis tested (model comparison) LRT p-value (L. ratio) 
Fixed effects Random effects 
Mixed models 
 MR.1 forumla forumla 1 662.5 (REML) −825.2 (REML) [6] – – 
 MR.2 As in MR.1 forumla 1 661.6 (REML); −825.8 (REML) [5] Random effects are independent; ρ = 0 (MR.1 vs. MR.20.30 (1.1) 
 MR.3 forumla As in MR.2 1 595.5 (REML); 1 578.9 (ML) −772.7 (REML); −764.4 (ML) [25] Residual errors are stock-specific (MR.3 vs. MR.2) <0.01 (106.2) 
 MR3.HL forumla forumla, forumla 1 576.3 (ML) −762.2 (ML) [26] Linear dependence of βi on habitat size (H); fH1 = 0 (MR3.HL vs. MR.3) 0.03 (4.6) 
 MR3.HQ forumla As in MR3.HL 1 561.2 (ML) −753.6 (ML) [27] Quadratic dependence of βi on habitat size (H); fH2 = 0 (MR3.HQ vs. MR3.HL) <0.01 (17.2) 
 MR3.HQ.TAL forumla forumla, forumla1 562.4 (ML) −753.2(ML) [28] Linear dependence of αi on temperature (T); forumla (MR3.HQ.TAL vs. MR3.HQ.TAQ<0.01 (14.3) 
 MR3.HQ.TAQ forumla as in MR3.HL.TAL 1 550.2 (ML); 1 606.8 (REML) −746.1 (ML); −774.4 (REML) [29] Quadratic dependence of αi on temperature (T); forumla (MR3.HQ.TAQ vs. MR3.HQ) <0.01 (15.0) 
 MR3.HQ.TBL forumla forumla, forumla 1 562.6 (ML) −753.3 (ML) [28] Linear dependence of βi on temperature (T); forumla (MR3.HQ.TBL vs. MR3.HQ0.46 (0.6) 
 MR3.HQ.TBQ forumla As in MR3.HQ.TBL 1 562.2 (ML) −752.3 (ML) [29] Quadratic dependence of βi on temperature (T); forumla (MR3.HQ.TBQ vs. MR3.HQ0.23 (2.9) 
 MR3.HQ.TAQ.RS forumla forumla, forumla,forumla, forumla 1 609.4 (REML) −773.7 (REML) [31] Across-stocks differences in the influence of T on alpha; forumla (MR3.HQ.TAQ.RS vs. MR3.HQ.TAQ0.5 (1.4) 
 MR3.HQ.AC forumla As in MR3.HQ 1 398.3 (ML) −651.2 (ML) [48] – – 
 MR3.HQ.TAQ.AC forumla As in MR3.HQ.TAQ 1 386.6 (ML) −643.3 (ML) [50] Quadratic dependence of αi on temperature (T) after first-order differencing of yit; forumla (MR3.HQ.TAQ.AC vs. MR3.HQ.AC) <0.01 (15.7) 
 MR3.HL.r As in MR3.HL, but fitted only to stocks with low autocorrelation in yit 665.3 (ML) −317.7 (ML) [15] – – 
 MR3.HL.TAL.r forumla As in MR3.HQ.TAQ 660.7 (ML) −314.4 (ML) [16] Linear dependence of αi on temperature (T); forumla (MR3.HL.TAL.r vs. MR3.HL.r) 0.01 (6.6) 

 
  DIC ΔDIC   

 
Bayesian models 
 BR3.HQ As in MR3.HQ 1 508    
 BR3.HQ.TAQ As in MR3.HQ.TAQ 1 489 −19   
 BR3.HQ.TAQ.RS As in MR3.HQ.TAQ.RS 1 480 −9   
 BR3.HQ.TAN.RS As in MR3.HQ.TAQ.RS, but using the inverse polynomial [Equation (9c)] to describe the dependence of alpha on T 1 500 +20   

For every model, the AIC, the restricted log-likelihood (logLik; estimated using ML and/or REML), and the number of parameters (d.f.) are given. For the model comparisons, the p-values and the statistic (L. ratio) of the LRT are given. Models found to be superior in each pair of comparisons are emboldened.

The final MR model (MR3.HQ.TAQ) includes non-linear H and T effects on beta and alpha, respectively, common T-related slopes (forumla and forumla) across stocks, independent random effects, and heterogeneous (stock-specific) errors σyi (Table 2, Figure 2). The model was also implemented in the Bayesian framework (BR3.HQ.TAQ; Table 2), the only difference being that the H effects were introduced directly on beta (i.e. on Smax).

Figure 2.

Stock-specific residual standard errors (±s.e.) obtained from the Ricker-Bayesian model plotted vs. sample size.

Figure 2.

Stock-specific residual standard errors (±s.e.) obtained from the Ricker-Bayesian model plotted vs. sample size.

Although no significant differences were found among stocks for the T-related slopes (cT1i and cT2i), the terms were allowed to be stock-specific in a Bayesian model (BR3.HQ.TAQ.RS) by introducing a stock-level model on the corresponding parameters. This choice was also supported by the DIC model-selection criterion (1480 vs. 1489 of BR3.HQ.TAQ; Table 2). Therefore, the model was allowed to estimate the best possible parameters for each stock, while accounting for uncertainty. For comparison, the model omitting the T effect on alpha was also fitted. The DIC of that model was 1508, considerably higher than the DIC of BR3.HQ.TAQ.RS (1480; Table 2).

The Ricker model was also fitted using the inverse polynomial [Equation (9c)] to describe the dependence of alpha on temperature (BR3.HQ.TAN.RS). The Bayesian framework was used to implement this model, because the relationship is non-linear. The results showed that the quadratic dependence used in the BR3.HQ.TAQ.RS provided a superior fit (Table 2). The final Bayesian model (BR3.HQ.TAQRS) used to produce the results presented next becomes  

14
formula
where forumlaforumlaforumlaforumla and forumla The model provided satisfactory results when tested for distributional assumptions (independence of within-stock errors and normality of random effects).

Autocorrelation

Autocorrelation at lag 1 in log(Rt/St) was significant for certain stocks (Table 1). To identify whether autocorrelation was responsible for the improved goodness-of-fit of the MR3.HQ.TAQ model, including the T effect on the alpha parameter, two different approaches were employed. Initially, first-differencing for these stocks was used by introducing log(Rt−1/St−1) as an explanatory variable in the model, and allowing the corresponding coefficients to be stock-specific. The ML models with and without the T effect (MR3.HQ.TAQ.AC and MR2.HQ.AC, respectively) were fitted and compared with LRT. The test shows that the T-related terms remained significant (Table 2), although first-differencing decreases the statistical power by increasing the Type-II error rate (Pyper and Peterman, 1998). The second approach involved fitting MR3.HQ and MR3.HQ.TAQ only to the ten stocks exhibiting a low degree of autocorrelation. However, reducing the number of stocks resulted in non-significant, second-order terms related to the T and H effect, forumla and forumla respectively. Instead, therefore, the models MR3.HL.r and MR3.HL.TAL.r, including only the linear terms, were compared. The LRT also revealed that, in that case, including the T effect would improve the model fit, with p = 0.01 (Table 2).

Temperature effects on Ricker alpha

The identified species-level relationship between alpha and T can be estimated based on the means of the related terms, forumlaforumla and forumla (Figure 3a). The relationship shows a dome-shaped functional response to temperature; the curve is nearly flat between 4.5 and 6°C and peaks at ∼5°C, whereas alpha declines sharply at temperatures outside the 4.5–6°C range, roughly corresponding to the first and fourth quartiles of cod spawning-season T. The within-stocks alphaT relationships have a similar functional form, depending on the experienced T range (Figure 3b, Table 3). The positive and/or negative effects of temperature on cod recruitment dynamics were also evident for many stocks, based on the coherence between the recruitment survival patterns and the temperature-driven fluctuations in alpha (Figure S3).

Figure 3.

(a) The relationship between alpha and temperature (black line) on the cod species level, estimated by the Ricker-Bayesian model. The grey lines correspond to the 95% credibility intervals. (b) The stock-specific relationships between alpha and temperature estimated by the Ricker-Bayesian model for the individual stocks. Error bars (±s.e.) are also plotted for the observations at lowest and highest temperature. See Figure 2 for stock symbol codes.

Figure 3.

(a) The relationship between alpha and temperature (black line) on the cod species level, estimated by the Ricker-Bayesian model. The grey lines correspond to the 95% credibility intervals. (b) The stock-specific relationships between alpha and temperature estimated by the Ricker-Bayesian model for the individual stocks. Error bars (±s.e.) are also plotted for the observations at lowest and highest temperature. See Figure 2 for stock symbol codes.

Table 3.

Mean, maximum, and minimum stock-specific alpha values estimated by the Ricker multilevel model.

Stock Mean Minimum Maximum Change (%) 
cod-2224 2.84 2.59 2.91 −6.77 
cod-2532 1.87 1.50 2.00 −17.05 
cod-347d 4.52 4.25 4.70 −12.77 
cod-7ek 1.87 1.62 2.35 −23.72 
cod-arct 3.53 3.06 3.60 7.85 
cod-coas 0.64 0.39 0.95 −90.32 
cod-farp 2.08 1.65 2.33 −30.54 
cod-iceg 3.23 2.97 3.32 −9.01 
cod-kat 1.92 1.71 1.99 −11.56 
cod2j3kl 2.55 1.47 3.38 35.54 
cod3m 2.73 2.60 2.76 −2.99 
cod3no 1.16 0.29 1.42 36.63 
cod3pn4rs 1.05 0.53 1.70 53.11 
cod3ps 1.84 1.61 2.13 26.31 
cod4tvn 2.05 1.93 2.10 7.60 
cod4vsw 2.59 2.37 2.62 −0.54 
cod4x 1.59 1.55 1.60 5.15 
codgb 2.05 1.56 2.30 −22.36 
codgom 2.55 2.32 2.59 −6.43 
codvia 2.23 2.04 2.56 −29.59 
codviia 2.14 1.53 2.36 −20.75 
Stock Mean Minimum Maximum Change (%) 
cod-2224 2.84 2.59 2.91 −6.77 
cod-2532 1.87 1.50 2.00 −17.05 
cod-347d 4.52 4.25 4.70 −12.77 
cod-7ek 1.87 1.62 2.35 −23.72 
cod-arct 3.53 3.06 3.60 7.85 
cod-coas 0.64 0.39 0.95 −90.32 
cod-farp 2.08 1.65 2.33 −30.54 
cod-iceg 3.23 2.97 3.32 −9.01 
cod-kat 1.92 1.71 1.99 −11.56 
cod2j3kl 2.55 1.47 3.38 35.54 
cod3m 2.73 2.60 2.76 −2.99 
cod3no 1.16 0.29 1.42 36.63 
cod3pn4rs 1.05 0.53 1.70 53.11 
cod3ps 1.84 1.61 2.13 26.31 
cod4tvn 2.05 1.93 2.10 7.60 
cod4vsw 2.59 2.37 2.62 −0.54 
cod4x 1.59 1.55 1.60 5.15 
codgb 2.05 1.56 2.30 −22.36 
codgom 2.55 2.32 2.59 −6.43 
codvia 2.23 2.04 2.56 −29.59 
codviia 2.14 1.53 2.36 −20.75 

The column labelled “Change” refers to the change in mean alpha induced by an increase in current mean T by 3°C.

No significant differences were identified among the slopes cT1i and/or cT2i describing the non-linear dependence of alpha on T (Figure S4). Therefore, the critical temperature, the point after which the negative effects on alpha prevail, estimated as −cT1i/2cT2i, does not differ significantly across stocks and can be considered equal to the species estimate of ∼5°C (Figure 3a). However, stocks are expected to differ in the alpha rate of change [i.e. as given by differentiating Equation (9a)]: dαit/dT = cT1i + 2TpcT2i, where Tp is their current temperature. Therefore, the expected proportional change in αit, dαit resulting from, for example, a 3°C increase in the current average temperature of each stock can be estimated (Figure 4a, Table 3). Rates of change in alpha are generally positive for stocks with current mean spring temperature <4°C, but become increasingly negative above ∼5°C (Figure 4a, Table 3).

Figure 4.

(a) The ratios (±s.e.) between alpha estimated at mean stock-specific T and the alpha corresponding to a 3°C increase plotted against current mean temperature. The same result applies also to Keq. (b) The corresponding plot for Kmax.

Figure 4.

(a) The ratios (±s.e.) between alpha estimated at mean stock-specific T and the alpha corresponding to a 3°C increase plotted against current mean temperature. The same result applies also to Keq. (b) The corresponding plot for Kmax.

The stock-specific intercepts coi represent alpha at a common T (equal to the overall mean of the T observations) across stocks, hence represent the among-stocks differences in alpha left unexplained by the variability in T. There are significant differences between stocks: cod-coas, cod3no, and cod3pn4rs have the lowest, and cod-347d, cod-arct, cod2j3kl, and cod-iceg the highest estimates compared with the across-stocks mean, forumla (Figure 5). It is also clear that the deviations are greater, and negative usually, for stocks at intermediate or lower mean spring temperatures (Figure 5).

Figure 5.

Stock-specific coi estimates obtained from the Ricker (black error bars) and BH (grey error bars) hierarchical Bayesian models. The error bars correspond to 95% credibility intervals. The stocks are ordered by increasing mean T (right y-axis).

Figure 5.

Stock-specific coi estimates obtained from the Ricker (black error bars) and BH (grey error bars) hierarchical Bayesian models. The error bars correspond to 95% credibility intervals. The stocks are ordered by increasing mean T (right y-axis).

Evidence for temperature effects in no-pooling Ricker SR models

For comparison, individual-stock (no-pooling) Ricker models were also fitted, assuming either linear or quadratic dependence of alpha on T. In the former case, significant (p < 0.1) negative effects were found for three stocks (cod-347d, cod-farp, and codgb) and positive effects for cod3no and cod2j3kl, and in the latter (p < 0.1 for the quadratic term) for cod-arct, cod-coas, and codviia. It is notable that most of these stocks are located in the upper T range, whereas the effect was not significant for most stocks, with a considerable proportion of the observations in the mid-range 4.5–6°C. Also, pooling of the T-related slopes towards the mean in the multilevel Bayesian model is stronger for stocks with less data in the analysis [Figure 6a and b; see also Equation (8)].

Figure 6.

The ratio between the T-related slopes obtained from the no-pooling Ricker models assuming quadratic dependence of alpha on T and the corresponding slopes obtained from the Bayesian Ricker model [(a) cT1i, (b) cT2i] plotted against sample size. The pooling is less (ratio close to 1) for stocks with more data. See Figure 2 for stock symbol codes.

Figure 6.

The ratio between the T-related slopes obtained from the no-pooling Ricker models assuming quadratic dependence of alpha on T and the corresponding slopes obtained from the Bayesian Ricker model [(a) cT1i, (b) cT2i] plotted against sample size. The pooling is less (ratio close to 1) for stocks with more data. See Figure 2 for stock symbol codes.

The dependence of Ricker beta on habitat size

The spawner biomass producing maximum recruitment, according to the Ricker SR model (forumla), depends on H (log of habitat size, defined as the area between 40 and 300 m; Figure 7a). The parameter is relatively constant when the log of habitat size is below the across-stocks mean, but then increases exponentially (Figure 7a). The model was also explored in terms of the pooling introduced in the beta parameter by the stock-level model, forumla. Pooling of the individual estimates towards the stock-level model predictions is stronger in the cases where less information is available, i.e. for stocks with smaller sample size [Figure 7c, Equation (8)]. The pooling (or shrinkage) also gives plausible estimates of the parameters (positive values) for all stocks, even when the individual SR model gives results (cod-coas, cod3m, and cod3no) that are meaningless (i.e. negative estimates for forumla which has units of S).

Figure 7.

The fitted stock-level model of beta as a function of H (log habitat size) in (a) the Ricker (corresponding to forumla and (b) the BH (corresponding to forumla) multilevel Bayesian models. (c) The ratio of beta estimated from the Bayesian-Ricker model (HBR) and the corresponding estimate obtained from individual, stock-specific, Ricker models (SSR) plotted against sample size. The ratio is closer to 1 (representing less pooling) for stocks with more observations. (d) The beta (expressed as 2KBH for the BH and as forumla for the Ricker model; '000 t km−2) parameter estimates and 95% credibility intervals obtained from the Bayesian Ricker (black error bars) and BH (grey error bars) hierarchical models, ordered by increasing mean temperature (right y-axis). The vertical black and grey lines correspond to the across-stocks means of the Ricker and BH beta parameters, respectively. See Figure 2 for stock symbol codes.

Figure 7.

The fitted stock-level model of beta as a function of H (log habitat size) in (a) the Ricker (corresponding to forumla and (b) the BH (corresponding to forumla) multilevel Bayesian models. (c) The ratio of beta estimated from the Bayesian-Ricker model (HBR) and the corresponding estimate obtained from individual, stock-specific, Ricker models (SSR) plotted against sample size. The ratio is closer to 1 (representing less pooling) for stocks with more observations. (d) The beta (expressed as 2KBH for the BH and as forumla for the Ricker model; '000 t km−2) parameter estimates and 95% credibility intervals obtained from the Bayesian Ricker (black error bars) and BH (grey error bars) hierarchical models, ordered by increasing mean temperature (right y-axis). The vertical black and grey lines correspond to the across-stocks means of the Ricker and BH beta parameters, respectively. See Figure 2 for stock symbol codes.

The extent of across-stock variation in beta explained by differences in the habitat size can also be estimated using the r2-based statistic (Gelman and Hill, 2007):  

formula
where forumla and V is the finite-sample variance operator across stocks. The numerator represents the average variance in forumla, i.e. in the average variance of forumla left unexplained by H, and the denominator the average variance among the stock-specific beta values. Note that the expectations are averaging over the uncertainty of the model using posterior simulations, so leading to a lower estimate, comparable with the traditional adjusted r2 (Gelman and Hill, 2007). The intermediate value of r2 obtained (0.48) shows that habitat size can explain almost half the observed variation in beta among cod stocks. In particular, the beta values for stocks located in areas with low or intermediate average T tend to be higher than those predicted by the model (Figure 7a). This can also be shown by plotting the beta values estimated on a per-unit-area basis (Figure 7d). This comparison shows that stocks with higher estimates of beta per km2 tend to be located in waters with intermediate mean temperature.

Effects of temperature and habitat size on carrying capacity

The T effects on alpha also have implications for the K-related parameters Kmax and Keq, which depend on both alpha and beta Ricker parameters, as described above. Therefore, for a given stock, K is time-varying, following the fluctuations of alpha with T (Tables S1 and S2). The proportional change in the average Kmax induced by a 3°C increase in the mean stock-specific temperatures is given by exp(dαit) (Figure 4b, Table S2), whereas the change in Keq is equal to dαit (Figure 4a, Table 3). It follows that the across-stocks pattern in both K indices is similar to that for dαit/dT, with the impact becoming progressively more negative in areas currently having higher mean temperatures, whereas for stocks inhabiting colder water, it is expected that the change will be positive, other factors remaining equal (Figure 4a and b).

On a per-unit-area basis, it is evident that there are extensive differences (∼30-fold) in K across stocks (Figure S5, and Tables S1 and S2), similar to the respective pattern demonstrated by the beta values (Figure 7d). This was expected because the functional form of the relationship between K and habitat size is determined by the exponential-like, non-linear model of beta with H (Figure 7a), accounting for part of the variability across stocks. The variability in K is also driven by the across-stocks differences in the alpha values, represented by the pattern in the intercepts coi (Figure 5). Therefore, the stocks with higher K per unit area (e.g. cod in 3M, Iceland, North Sea, Kattegat, and the western Baltic; Figure S5) are among those with the higher estimates of beta per km2 (Figure 7d) and coi (Figure 5).

Comparison of BH and Ricker multilevel SR model results

The multilevel BH model, analogous to the final Ricker SR model given by Equation (14), is  

15
formula
where forumlaforumlaforumla and forumla

The two SR models, Ricker and BH, provided similar estimates for the alpha-related terms coi (Figure 5), cT1i, and cT2i (Figure S4). Consequently, the critical T and dαit, which depend on the previous terms, are not significantly different. In addition, the form of the stock-level model between forumla (beta) and log habitat size (H) is analogous to the corresponding Ricker submodel (Figure 7b). However, the BH submodel has a better fit, as quantified by the r2 statistic, explaining ∼70% of the across-stocks variability in the parameter. The model provided, in general, higher point estimates for beta. The difference, however, is considerable only in a few cases, namely for cod-iceg, cod4tvn, cod2j3kl, and cod-2532 (Figure 7d). Following the pattern in beta, there are substantial differences in the mean estimates of Keq provided by the two SR models (Table S1) and to a lesser extent in Kmax (Table S2). The relationship between K and log habitat is non-linear, as in the Ricker model, so per-unit-area comparisons of the parameters are not straightforward.

Reduction in parameter uncertainty

There was a substantial increase in precision for both alpha and beta parameters when estimated from the HBR compared with the SSR models, especially for stocks with shorter time-series (Figure 8a). The reduction in uncertainty was stronger for beta, for which the average reduction was 70%, and highest for northern Gulf of St Lawrence (cod3pn4rs), northern (cod2j3kl), and Celtic Sea (cod-7ek) cod. For alpha (estimated at mean stock-specific T), the decrease in uncertainty was 44%, and it was higher for Grand Bank (cod3no), Norwegian coastal (cod-coas), Celtic Sea (cod-7ek), and Flemish Cap (cod3m) cod stocks.

Figure 8.

(a) Stock-specific ratios of the estimation uncertainty indicators in the SSR (without temperature effects) to HBR models, for alpha (grey dots) and beta (black dots). Estimation uncertainty indicators for alpha and beta are defined as s.e./mean, i.e. expressing the standard error of a parameter in units of the mean. In the HBR, alpha estimates correspond to average stock-specific temperature. Estimation precision ratios for beta have been omitted in cases when the individual models gave implausible estimates (cod3m, cod3no, and cod-coas). The stocks are ordered by increasing sample size (from the bottom to the top). The uppermost larger symbols and the corresponding solid grey and black lines indicate the average of the alpha and beta ratios, respectively. The dashed line indicates unity, i.e. equal estimation precision between the models. (b) Ratios of the ESS (error sum of squares) obtained from the SSR to the HBR model against sample size. See Figure 2 for stock symbol codes. The dashed line indicates unity, i.e. equal ESS between the models.

Figure 8.

(a) Stock-specific ratios of the estimation uncertainty indicators in the SSR (without temperature effects) to HBR models, for alpha (grey dots) and beta (black dots). Estimation uncertainty indicators for alpha and beta are defined as s.e./mean, i.e. expressing the standard error of a parameter in units of the mean. In the HBR, alpha estimates correspond to average stock-specific temperature. Estimation precision ratios for beta have been omitted in cases when the individual models gave implausible estimates (cod3m, cod3no, and cod-coas). The stocks are ordered by increasing sample size (from the bottom to the top). The uppermost larger symbols and the corresponding solid grey and black lines indicate the average of the alpha and beta ratios, respectively. The dashed line indicates unity, i.e. equal estimation precision between the models. (b) Ratios of the ESS (error sum of squares) obtained from the SSR to the HBR model against sample size. See Figure 2 for stock symbol codes. The dashed line indicates unity, i.e. equal ESS between the models.

HBR model performance

The fitted SR curves obtained from the HBR model were better able to capture some of the trends in recruitment than an SSR model without temperature effects (Figure 9). The confidence intervals of the fit were also narrower, especially at higher spawner biomass, because of the greater estimation accuracy in the alpha and beta parameters (Figure 8a). For most stocks, the predicted curves seem to follow the patterns rather well (e.g. cod-7ek, cod-arct, cod-farp, cod-iceg, cod-kat, cod-2224, cod2j3kl, cod3pn4rs, cod3ps, codgb, codviia). The error sum of squares obtained from the HBR was reduced for many stocks compared with the SSR models (Figure 8b). This can be an effect of partial pooling which, by definition, assigns more weight to observations closer to the mean. Alternatively, for those stocks, T and S fluctuations could not explain effectively the variation in recruitment and that other factors were driving the patterns instead. It is also notable that, for most, there was substantial coherence between the recruitment survival and alpha (Figure S3), underlying the significance of thermal effects on cod productivity.

Figure 9.

Fitted Ricker SR curves (solid lines) and 95% confidence intervals (dashed lines). Grey lines correspond to the SSR model fits without the effect of temperature, and black lines to the Bayesian model with the effect of temperature (HBR). Symbols are observed data. Note that recruitment is expressed in replacement spawners.

Figure 9.

Fitted Ricker SR curves (solid lines) and 95% confidence intervals (dashed lines). Grey lines correspond to the SSR model fits without the effect of temperature, and black lines to the Bayesian model with the effect of temperature (HBR). Symbols are observed data. Note that recruitment is expressed in replacement spawners.

Discussion

The main objective of our study was to investigate the effects of ecosystem properties (e.g. spawning-season temperature and nursery-ground size) on cod alpha and beta SR parameters across the North Atlantic distribution of the species. We employed two complementary hierarchical methods, mixed modelling and Bayesian inference, to combine data on all cod stocks in an analysis of both Ricker and BH SR relationships. We found two main results.

First, there is a dome-shaped relationship between alpha, and hence also K, and temperature, identified using both the Ricker and BH SR models (Figure 3). Consequently, the influence of T is positive in colder waters down to ∼5°C, which is close to the mean of the cod spawning-season thermal range and the temperature at which cod egg survival is greatest in laboratory experiments (Pepin et al., 1997), and negative for stocks inhabiting warmer water. These findings imply that ocean warming will cause considerable impacts on both alpha and K, and hence on cod SR dynamics. More importantly, the models allow inference on both the sign and the extent of these effects at both species and individual stock level. Regarding the effect of habitat, beta (representing forumla in the Ricker, and forumla in the BH model), which is also a K component and also describes density-dependent regulation, depends on area size following an exponential-like relationship, partly causing the observed across-stocks differences in K per unit area. Additional biological interpretations of these results follow below.

Our second main result is that we quantified how much uncertainty in SR model parameters could be reduced by incorporating environmental information directly in the hierarchical Bayesian SR model framework. In some cases, it was possible to produce parameter estimates for stocks whose estimates in a single-stock (no-pooling) framework were otherwise ecologically meaningless or statistically insignificant. Hence, in cases where no functional relationship exists between spawner biomass and recruitment at the single-stock level, Bayesian-based approaches can be a means to estimate and predict the dynamics of these stocks under exploitation and environmental (temperature) scenarios (Hilborn and Liermann, 1998; Myers et al., 2001). Moreover, the uncertainty in the alpha and beta parameters of Ricker models can be reduced, on average, by 44 and 70%, respectively, if they are estimated in a hierarchical framework (Figure 8a). Given that uncertainty continues to be a major impediment to reliable forecasting of fish population dynamics, especially in changing environments and ecosystem structures (ICES, 2007a), our results could lead to more reliable forecasts (i.e. less uncertainty) and scenarios of stock dynamics (e.g. recoveries) under different exploitation and temperature scenarios.

The study was based on stock-and-recruit data extracted from published stock assessments. Therefore, and unavoidably, it suffers from shortcomings inherent to that framework and hence common to the SR modelling approaches (Hilborn and Walters, 1992; Needle, 2002). One important source of non-stationarity in SR models stems from the way stocks are defined for operational purposes, which may or may not agree with the “stock” concept from an ecological and/or evolutionary perspective (Hilborn and Walters, 1992; Waldman, 1999; Waples and Gaggiotti, 2006). For management, stocks are assumed to consist of discrete, homogeneous units, with negligible exchange with neighbouring stocks. However, such exchanges (e.g. migration, source/sink processes related to the drift of eggs and larvae) do occur (Metcalfe, 2006). Moreover, some stocks consist of several subpopulations within a management area and whose abundance and productivity might change over time (Robichaud and Rose, 2001; Metcalfe, 2006; Heath et al., 2008). Both these mechanisms can introduce uncertainty to SR models for single stocks, e.g. by increasing observation error in the estimation of spawning stock (Hilborn and Walters, 1992).

Also, the stock and recruit time-series are not direct observations, but estimates provided by VPA analyses, based on data subject to observation error. Therefore, the SR data are subject to uncertainty, introducing measurement error or error-in-variables bias in the SR models (Walters and Ludwig, 1981). In addition, spawner biomass is not always fully representative of the annual rate of egg production, because the age and length structure of a stock, also driven by exploitation, play an important role in reproductive output and success (Trippel, 1999; Marshall et al., 2003) and can also determine the sensitivity of recruitment to climate (Perry et al., 2009). In general, however, the measurement error in S does not affect seriously the alpha estimates in the Ricker model (Kehler et al., 2002). Second, the correlation of S levels with previous deviations of R, through system dynamics, introduces time-series bias in the analysis and can result in underestimation of optimum stock size (Walters, 1985). Nevertheless, it has been shown that time-series bias is usually of less importance when Ricker SR models are applied to cod stocks (Myers and Barrowman, 1995). A relatively new approach that can effectively address the above problems, at least for semelparous species, is the application of state–space methodology in a Bayesian framework for modelling SR relationships (Meyer and Millar, 2001). The application of this method to iteroparous species is not straightforward unless estimates of the observation error are available or can be assumed because the dependence of S and R on previous observations is more complex. Potential improvements in the present approach could involve changes in the estimation of the SPRF=0 parameters and of the habitat size, as discussed below.

Hierarchical modelling approaches in SR studies

A hierarchical multilevel approach offers a number of advantages which have been demonstrated to be particularly useful for the analyses of fisheries data (e.g. Hilborn and Liermann, 1998; Myers, 2002; Michielsens and McAllister, 2004; Su et al., 2004). The methods are based on the stock-level models describing the variation among stock-specific parameters across the species range. For beta, these models are extended to include habitat size as an individual stock predictor, which can partly explain the observed variation. In a Bayesian framework, the stock-level models can also be used to introduce existing knowledge into the model (Gelman and Hill, 2007). Here, however, an empirical Bayesian approach was used wherein priors are uninformative or, for beta, depend on H (log habitat size).

An alternative method to model across-stocks variation in parameters would have been to fit separate (no-pooling) models to each stock and then to model the stock-specific parameters. Multilevel methods, however, combine these two stages in a single model, so that inference is based on both within- and among-stocks variability, incorporating uncertainty in all parameters. Therefore, the stock-level models are used to convey information about the probability distribution of the parameter estimated by all-stocks data. As the inference about single-stock parameters is based on these priors, strength is “borrowed”, and the “loan” (pooling) is higher for “poorer” (limited or highly variable observations) stocks. Consequently, the empirical Bayesian inference is superior to that of no-pooling models, especially when the latter provide implausible estimates or lack the required power to demonstrate the significance of certain terms. Our analyses have provided several examples where the pooled parameter estimates or their uncertainty improved relative to the single-stock estimates.

Our hierarchical modelling approach also incorporated an alternative formulation of SR–environment models to evaluate specifically how the parameters of SR models vary with ecosystem properties. This approach disaggregates the parameter into a simple regression model, using the environmental effects as predictors [e.g. Equations (9a) or (9c) and (9b)]. Mathematically, the approach is similar to, for example, including additional environmental terms in the log-transformed version of the Ricker model. However, disaggregating the parameters explicitly and modelling them in response to ecosystem properties can be beneficial. For example, it was possible to quantify and visualize the functional form of the relationship between alpha and temperature and to use the Bayesian framework to quantify the reduction in uncertainty of these parameter estimates.

Effects of temperature on alpha and carrying capacity

Bringing together data, and particularly temperature, across the entire North Atlantic distribution of cod has allowed inference on the functional form of the relationship between alpha and T at a species level. It is worth noting that the T effect remains significant even after first-order differentiation of the R/S data, although the method is known to increase Type-II error probability, when employed for low-frequency environmental signals (Pyper and Peterman, 1998). Temperature has opposite effects at the upper (negative) and lower (positive) thermal extremes, roughly corresponding to waters with T above or below 5°C, respectively. Owing to the quadratic form, the strength of the effect is stronger for temperatures closer to the extremes, and weakest at the middle, “neutral” 4.5–6°C interval (Figure 3). Similar geographic patterns have been observed in the response of cod recruitment to temperature across certain stocks (Planque and Frédou, 1999; Brander, 2000). Consequently, the effect of a potential increase in mean temperature will be more pronounced for stocks inhabiting areas closer to the distribution limits. Therefore, alpha for cod in the northern Gulf of St Lawrence (3Pn4RS), with the lowest mean T, could be expected to increase by >50%, whereas in the Celtic Sea, where mean T is highest, the decrease would be ∼24% (Figure 4a; Table 3). These predictions are, of course, only preliminary and assume that other ecosystem properties will react in future to such a temperature change in the manner they have in the past. When they do not, the model predictions may not be valid. An example of such a case is cod in the southern Gulf of St Lawrence, which collapsed in the early 1990s and is expected to be extirpated within 40 years even if fishing mortality is eliminated. The reason for the failed recovery and expected extirpation is an increase in the natural mortality of adult cod (Swain and Chouinard, 2008).

Notwithstanding, we note that the effect of temperature is not necessarily (only) a direct physiological effect on a particular life-history stage, but rather an ecosystem-level thermal indicator averaged over modest temporal and spatial scales, i.e. a regional sea during spring. The temperature data, therefore, are indices that integrate a large number of processes which affect early life history in cod (including direct physiological effects on, for example, growth rates, but also effects on the foodweb that might also affect survival), and whose net effects within and among entire stocks is reflected in the results detected here. The mechanisms through which temperature affects cod early life history are many, complex, and non-linear (Planque and Frédou, 1999; Sundby, 2000; Drinkwater, 2005; Rijnsdorp et al., 2009). They include the reproductive potential of spawners, e.g. spawning time, weight-at-age, adult condition, and maturation, and the biological processes controlling early life stages, such as incubation time, size-at-hatch, growth rate, and stage duration. In addition, temperature is a proxy for other ecosystem-level processes, whose importance varies across stocks, but includes wind-induced turbulence, food availability, displacement of spawning areas, and predator abundance (Sundby and Fossum, 1990; deYoung and Rose, 1993; Beaugrand et al., 2003; Houde, 2008).

Our observation that alpha peaked at T ∼5°C is consistent with other macro-ecological investigations of cod ecology. For instance, the average temperature in cod spawning locations across the North Atlantic (Sundby, 2000), the narrower temperature range experienced by cod during their spawning season, as revealed by tagging studies (D. Righton, Cefas, per. comm.), and laboratory experiments on temperature effects on cod eggs and larvae (Rombough, 1997; Jordaan and Kling, 2003), suggest that the optimal thermal range for cod early life stages is around 5–7°C. Therefore, we believe that the temperature datasets used in this study are representative of the thermal conditions influencing cod recruitment dynamics by acting on processes relevant across its distribution. This suggestion is also supported by the fact that no significant differences were found among the stock-specific, T-related slopes in the alpha models (Figure S4), suggesting that the impact of temperature across stocks is an adequate description of the relationship. Consequently, the estimated species-level, dome-shaped relationship between reproductive rate and temperature can be used as a prior for the parameterization of cod SR models designed to incorporate the potential influence of ocean warming, hence incorporating more biological knowledge in stock assessments and management.

In the present context, where we use standardized recruitment time-series (multiplied by SPRF=0) and allow for T-effects, alpha can be interpreted as the maximum rate at which spawners produce replacement spawners, in the absence of fishing mortality, given the temperature conditions during the spawning season in a particular year. Therefore, “maximum” should not be interpreted in absolute terms because it has a temporal, temperature-dependent component. Therefore, it is possible to produce a set of SR curves, corresponding to the mean, minimum, and maximum alpha estimates, that show the most pronounced differences.

The SPRF=0 parameter depends on weight-at-age, which is affected by T (Brander, 1995, 2007; Dutil and Brander, 2003). Therefore, it can be inferred that the relationship between alpha and T also describes SPRF=0 fluctuations, assuming that that parameter is temperature-dependent and varies around a mean value corresponding to the mean estimates used for recruitment standardization. Hence, an alternative approach would be to use time-varying estimates of SPRF=0 (Shelton et al., 2006), either to standardize recruitment or as an explicit component of SR models.

It was also demonstrated that factors other than temperature influence alpha, causing the across-stocks variability in the intercepts of the alphaT relationship (Figure 5). The differences are more pronounced among stocks located in colder water, whereas in warmer areas, T seems to be the limiting factor. Indeed, there is evidence that environmental variables, such as the North Atlantic Oscillation (Stige et al., 2006), affect cod recruitment across the North Atlantic. In addition, fishing can impact the spawning potential of a stock by altering age- and size-at-maturity and/or growth rate (Ottersen et al., 2006; Perry et al., 2009; Planque et al., 2009). Biotic interactions with prey, predator, or competing species affect the productivity and survival rates of both early stage and adult cod (Lilly et al., 2008; Swain and Chouinard, 2008). The reasons for the among-stocks variability in alpha, whether of local- or species-level importance, bear further investigation. It would also be pertinent to employ similar methods to study cod response to other forcing factors, e.g. primary productivity, and to compare the magnitude and functional response to such factors with those documented here for temperature.

Finally, we note that our results for the parameter alpha agree with previous studies employing similar models to study SR dynamics (Myers et al., 1999, 2001), which have concluded that alpha is relatively constant across cod stocks. The average variation among the stock-specific alpha parameters in our study is approximately sevenfold, and the mean (corresponding to mean T) estimates are comparable with those obtained by Myers et al. (2001) applying a non-linear mixed-BH model to an older version of the cod dataset and using slightly different estimates of SPRF=0. Moreover, similar studies have suggested that the across-cod-stocks differences in alpha are related to mean bottom temperature in each region (Myers et al., 2002).

Carrying capacity, defined as the equilibrium S (Keq) or as the S producing the maximum replacement spawner biomass (Kmax), depends on alpha and is a dynamic variable varying temporally according to the fluctuations in alpha driven by T. Further, K is an exponential function of alpha, except for Keq obtained from the Ricker model, so temperature effects are more severe for K. For example, if mean T increased by 3°C, Kmax for northern cod (2j3kl) would be expected to increase approximately 2.5-fold, compared with a ∼35% increase in alpha (Figure 4). In the Celtic Sea, the reduction in K would be ∼44%, whereas the decrease in alpha was estimated at 24%.

Effects of habitat size on beta

Habitat size was defined as the area suitable for the settlement of juveniles and, hence, for density-dependent processes. Owing to the log-transformation of H, the model is quite robust to exact definitions of habitat size. The definition of H was able to represent some of the density-dependent processes in both BH and Ricker cases; this finding is, therefore, consistent with the results of earlier studies that showed that beta (=Smax), on a per-km2 basis, differed significantly among cod stocks (Myers et al., 2001), and that habitat size is a significant correlate of mean recruitment level in North Atlantic cod stocks (MacKenzie et al., 2003). The variance explained (r2) by the relationship is higher using the BH model (70 vs. 48% for the Ricker model), indicating that the employed habitat definition is mostly representative of the density-dependent mechanisms incorporated in that SR model. The dependence of beta on the log-transformed habitat size, under both SR models, was described by a quadratic submodel [Equation (12)], showing an exponential type of increase with increasing H for the most part (i.e. Figure 7a and b). Other types of submodel were also considered, but the posteriors of the parameters showed poor convergence.

Example case study interpretation: temperature and North Sea cod

Our results on the impacts of temperature on alpha and beta should be compared with other studies of cod ecology and might help interpret past dynamics of individual cod stocks. As an example, we briefly consider the North Sea cod stock, which has declined in recent decades as a consequence of both fishing and changing ecosystem conditions (especially temperature and zooplankton; Beaugrand et al., 2003; Kell et al., 2005a). We note initially that the long-term mean average condition in ten cod stocks in the North Atlantic is higher in warmer waters and that these stocks, including the North Sea stock, should therefore be able to withstand higher levels of exploitation (Rätz and Lloret, 2003). Moreover, alpha from Ricker SR models for these ten stocks was positively and linearly correlated with cod condition (Rätz and Lloret, 2003). Our analyses showed that alpha decreased at temperatures exceeding ca. 5°C. When considered together, our results and those of Rätz and Lloret (2003) suggest that, at warmer temperatures, the processes affecting different components of stock productivity, e.g. adult condition, growth, production of eggs and larvae, and survival of juveniles, may have different functional responses to temperature.

For example, in the North Sea, the negative temperature influence on recruitment documented here and by others (Planque and Frédou, 1999; Brander, 2000) is believed to be related to larval/pelagic 0-group feeding processes and zooplankton community dynamics (Beaugrand et al., 2003). However, the analysis of interactions between adult condition, alpha, and temperature (Rätz and Lloret, 2003) shows that cod condition should still be high in the North Sea even at warm temperatures. In addition, the scaling of the temperature variable in the different analyses was different: our study used temperature at spawning time, whereas the study of Rätz and Lloret (2003) employed annual means. The latter scaling may account for temperature effects during winter on cod condition, for example, which may yield a different functional response than the effect of spring temperature on cod recruitment and early juvenile survival.

Nevertheless, the current state of the North Sea cod stock suggests that exploitation has been so heavy that the beneficial effect of warm temperature on North Sea cod condition has been offset by greater fishing mortality of adults and juveniles and a reduced production of recruits associated with changes in the plankton community. Several modelling analyses suggest that less-intense exploitation would allow the stock to recover, despite the negative impact of warmer temperature on cod recruitment (Cook and Heath, 2005; Kell et al., 2005a). These modelling results are supported by archaeological studies which have shown that cod were common in coastal areas of the North Sea (and the neighbouring Kattegat) in the Atlantic Warm Period (ca. 7000–3900 BC) although temperatures were 2–3°C higher than late 20th century temperatures (Enghoff et al., 2007) and, therefore, similar to those expected near the end of the 21st century associated with global climate change. Moreover, recent tagging studies show that adult cod occasionally occupy warm water in the North Sea (Neat and Righton, 2007). Given these findings, we believe that reduced exploitation could help sustain North Sea cod populations well into the 21st century.

Perspectives and conclusions

Multilevel models are particularly useful for identifying and quantifying processes acting on a broad scale, such as in determining fish population dynamics across their distribution. We have described the response of cod recruitment and K to temperature, allowing for the effect of habitat size to approximate density-dependent mechanisms. It would also be interesting to develop similar approaches for other fish species. Applying hierarchical models to a study of environmental impacts on key forage or competing species of cod, either separately or combined in multispecies models, may also improve predictions of the implications of ocean warming for the structuring of local ecosystems and their trophic flows (Worm and Myers, 2003; Frank et al., 2007).

Developing hierarchical SR models under the Bayesian paradigm can also be useful for fisheries management, and especially for the implementation of the precautionary approach. A precautionary approach is based on defining limit and target reference points for the state of a stock and the exploitation rates, and requires formal consideration of scientific uncertainty in their estimation (FAO, 1995). For many ICES stocks, however, the management procedures ignore key sources of uncertainty, resulting in inconsistent or even inappropriate biomass and fishing mortality reference points (Kell et al., 2005b). Therefore, the development of operating models representing state-of-the-art knowledge of the underlying population dynamics and that are robust to uncertainty has been proposed to address this problem (Kell et al., 2005b). Also, it has been shown that only 17% of the ICES stocks actually had the necessary estimates of reference points needed to implement a precautionary harvest control rule and also that a large proportion of US stocks had uncertain stock status (Cadrin and Pastoors, 2008). A hierarchical Bayesian approach, such as that presented here, can therefore offer at least three advantages in terms of reference-point estimation: (i) integration over uncertainty in the parameters, including non-stationarity attributable to environmental effects, (ii) improved precision by borrowing strength from a multistock dataset and incorporating more biological knowledge in the form of priors, and (iii) derivation of consistent reference points also for data-poor stocks. In addition, the results of this study in terms of the dependence of cod maximum reproductive rate on temperature and of carrying capacity and density-dependence on habitat size can also be used as priors in similar studies, especially when investigating potential climate-change effects on cod recruitment.

Supplementary material

The following supplementary material is available at ICESJMS online: Hierarchical model estimation procedures: Bayesian inference and additional tables (S1–S2) and figures (S1–S5) as listed in the text.

Acknowledgements

This work is a contribution to the following EU projects: Marie Curie EST project METAOCEANS (MEST-CT-2005-019678), EurOceans Network of Excellence, UNCOVER (UNnderstanding the mechanisms of stock reCOVERy), and RECLAIM (REsolving CLimAtic IMpacts on fish stocks). We thank the Danish National Research Foundation (Danmarks Grundforskningsfond) for support to the Centre for Macroecology, Evolution and Climate (University of Copenhagen and DTU), Keith Brander (DTU-Aqua), Laurence Kell (ICCAT) and Bob Mohn (DFO Canada) for discussions and input, and especially Mette Bertelsen (ICES) and DFO scientists for their help in compiling the cod datasets. We are thankful to the anonymous reviewers for constructive comments and suggestions that contributed considerably to the improvement of the manuscript. We are also grateful to Else Juul Green (ICES) and Paul Dunphy (Maritimes) for assistance in extracting the temperature data and to Rasmus Borgstrøm and Andreas Espersen (both DTU-Aqua) for providing the estimates of habitat size. Funding to pay the Open Access publication charges for this article was provided by an EU Marie Curie Early Stage Training Network (Meta-oceans).

References

Beaugrand
G.
Brander
K. M.
Lindley
J. A.
Souissi
S.
Reid
P. C.
Plankton effect on cod recruitment in the North Sea
Nature
 , 
2003
, vol. 
426
 (pg. 
661
-
664
)
Berliner
L. M.
Hanson
K. M.
Silver
R. N.
Hierarchical Bayesian time series models
Maximum Entropy and Bayesian Methods
 , 
1996
Dordrecht
Kluwer Academic
(pg. 
15
-
22
480
Beverton
R. J. H.
Holt
S. J.
On the Dynamics of Exploited Fish Populations
1957
Fishery Investigations, London, Series II, 19
pg. 
533
 
Bishop
C. A.
Murphy
E. F.
Davis
M. B.
Baird
J. W.
Rose
G. A.
An assessment of the cod stock in NAFO Divisions 2J + 3KL
1993
NAFO Scientific Council Research Document, 93/86
Brander
K. M.
The effects of temperature on growth of Atlantic cod (Gadus morhua L.)
ICES Journal of Marine Science
 , 
1995
, vol. 
52
 (pg. 
1
-
10
)
Brander
K. M.
Effects of environmental variability on growth and recruitment in cod (Gadus morhua) using a comparative approach
Oceanologica Acta
 , 
2000
, vol. 
4
 (pg. 
485
-
496
)
Brander
K. M.
The role of growth changes in the decline and recovery of North Atlantic cod since 1970
ICES Journal of Marine Science
 , 
2007
, vol. 
64
 (pg. 
211
-
217
)
Brander
K. M.
Mohn
R. K.
Effect of the North Atlantic Oscillation on recruitment of Atlantic cod (Gadus morhua)
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2004
, vol. 
61
 (pg. 
1558
-
1564
)
Brander
K. M.
Mohn
R. K.
Erratum: effect of the North Atlantic Oscillation on recruitment of Atlantic cod (Gadus morhua)
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2004
, vol. 
61
 pg. 
2522
 
Brattey
J.
Cadigan
N. G.
Healey
B. P.
Lilly
G. R.
Murphy
E. F.
Shelton
P. A.
Mahé
C.
An assessment of the cod (Gadus morhua) stock in NAFO Subdivision 3Ps in October 2004
2004
DFO Canadian Science Advisory Secretariat Research Document, 2004/083
Cadrin
S. X.
Pastoors
M.
Precautionary harvest policies and the uncertainty paradox
Fisheries Research
 , 
2008
, vol. 
94
 (pg. 
367
-
372
)
Chouinard
G. A.
Currie
L.
Poirier
G. A.
Hurlbut
T.
Daigle
D.
Savoie
L.
Assessment of the southern Gulf of St Lawrence cod stock, February 2006
2006
Canadian Science Advisory Secretariat Research Document, 2006/006
Clark
D. S.
Gavaris
S.
Hinze
J. M.
Assessment of cod in Division 4X in 2002
2002
Canadian Science Advisory Secretariat Research Document, 2002/105
Clark
J. S.
Models for Ecological Data. An Introduction
2007
Princeton, NJ
Princeton University Press
pg. 
632
 
Cook
R. M.
Heath
M. R.
The implications of warming climate for the management of North Sea fisheries
ICES Journal of Marine Science
 , 
2005
, vol. 
62
 (pg. 
1322
-
1326
)
Demidenko
E.
Mixed Models: Theory and Applications
2004
Hoboken, NJ
Wiley
pg. 
736
 
deYoung
B.
Rose
G. A.
On recruitment and distribution of Atlantic cod (Gadus morhua) off Newfoundland
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1993
, vol. 
50
 (pg. 
2729
-
2741
)
Drinkwater
K. F.
The response of Atlantic cod (Gadus morhua) to future climate change
ICES Journal of Marine Science
 , 
2005
, vol. 
62
 (pg. 
1327
-
1337
)
Dutil
J-D.
Brander
K. M.
Comparing productivity of North Atlantic cod (Gadus morhua) stocks and limits to growth production
Fisheries Oceanography
 , 
2003
, vol. 
12
 (pg. 
502
-
512
)
Enghoff
I. B.
MacKenzie
B. R.
Nielsen
E. E.
The Danish fish fauna during the warm Atlantic period (ca. 7,000–3,900 BC): forerunner of future changes?
Fisheries Research
 , 
2007
, vol. 
87
 (pg. 
285
-
298
)
Fanning
L. P.
Mohn
R. K.
MacEachern
W. J.
Assessment of 4VsW cod to 2002
2003
Canadian Science Advisory Secretariat Research Document, 2003/027
pg. 
41
 
FAO.
Precautionary approaches to fisheries. 1. Guidelines on the precautionary approach to capture fisheries and species introductions
1995
FAO Fisheries Technical Paper, 350
pg. 
52
 
Frank
K. T.
Petrie
B.
Shackell
N. L.
The ups and downs of trophic control in continental shelf ecosystems
Trends in Ecology and Evolution
 , 
2007
, vol. 
22
 (pg. 
236
-
242
)
Fréchet
A.
Gauthier
J.
Schwab
P.
Pageau
L.
Savenkoff
C.
Castonguay
M.
Chabot
D.
, et al.  . 
The status of cod in the Northern Gulf of St Lawrence (3Pn, 4RS) in 2004
2005
Canadian Science Advisory Secretariat Research Document, 2005/060
pg. 
75
 
Gelman
A.
Carlin
J. B.
Stern
H. S.
Rubin
D. B.
Bayesian Data Analysis
2004
2nd edn.
London
CRC Press
Gelman
A.
Hill
J.
Data Analysis Using Regression and Multilevel/Hierarchical Models
2007
New York
Cambridge University Press
pg. 
625
 
Goodwin
N. B.
Grant
A.
Perry
A. L.
Dulvy
N. K.
Reynolds
J. D.
Life history correlates of density-dependent recruitment in marine fishes
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2006
, vol. 
63
 (pg. 
494
-
509
)
Gregory
D. N.
Climate: a database of temperature and salinity observations for the Northwest Atlantic
2004
Canadian Science Advisory Secretariat Research Document, 2004/075
Heath
M. R.
Kunzlik
P. A.
Gallego
A.
Holmes
S. J.
Wright
P. J.
A model of meta-population dynamics for North Sea and West of Scotland cod. The dynamic consequences of natal fidelity
Fisheries Research
 , 
2008
, vol. 
93
 (pg. 
92
-
116
)
Hilborn
R.
Liermann
M.
Standing on the shoulders of giants: learning from experience
Reviews in Fish Biology and Fisheries
 , 
1998
, vol. 
8
 (pg. 
273
-
283
)
Hilborn
R.
Walters
C. J.
Quantitative Fisheries Stock Assessment: Choice, Dynamics, and Uncertainty
1992
New York
Chapman and Hall
pg. 
570
 
Houde
E. D.
Emerging from Hjort's shadow
Journal of Northwest Atlantic Fishery Science
 , 
2008
, vol. 
41
 (pg. 
53
-
70
)
ICES.
Report of the Study Group on Precautionary Reference Points for Advice on Fishery Management
2003
ICES Document CM 2003/ACFM: 15
pg. 
338
 
ICES.
Spawning and life history information for North Atlantic cod stocks
2005
ICES Cooperative Research Report, 274
pg. 
154
 
ICES.
Report of the Working Group on the Assessment of Northern Shelf Demersal Stocks (WGNSDS), 9–18 May 2006
2006
ICES Document CM 2006/ACFM: 30
pg. 
870
 
ICES.
Report of the Workshop on the Integration of Environmental Information into Fisheries Management Strategies and Advice (WKEFA)
2007
ICES Document CM 2007/ACFM: 25
pg. 
182
 
ICES.
Report of the Working Group on the Assessment of Demersal Stocks in the North Sea and Skagerrak (WGNSSK), 5–14 September 2006
2007
ICES Document CM 2007/ACFM: 35
pg. 
1160
 
Jordaan
A.
Kling
L. J.
Browman
H. I.
Skiftesvik
A. B.
Determining the optimal temperature range for Atlantic cod (Gadus morhua) during early life
The Big Fish Bang. Proceedings of the 26th Annual Larval Fish Conference
 , 
2003
Bergen, Norway
Institute of Marine Research
(pg. 
45
-
62
)
Kehler
D. G.
Myers
R. A.
Field
C. A.
Measurement error and bias in the maximum reproductive rate for the Ricker model
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2002
, vol. 
59
 (pg. 
854
-
864
)
Kell
L. T.
Pilling
G. M.
Kirkwood
G. P.
Pastoors
M.
Mesnil
B.
Korsbrekke
K.
Abaunza
P.
, et al.  . 
An evaluation of the implicit management procedure used for some ICES roundfish stocks
ICES Journal of Marine Science
 , 
2005
, vol. 
62
 (pg. 
750
-
759
)
Kell
L. T.
Pilling
G. M.
O'Brien
C. M.
The implications of climate change for the management of North Sea cod (Gadus morhua)
ICES Journal of Marine Science
 , 
2005
, vol. 
62
 (pg. 
1483
-
1491
)
Kuparinen
A.
Merilä
J.
The role of fisheries-induced evolution
Science
 , 
2008
, vol. 
320
 (pg. 
47
-
48
)
Lilly
G. R.
Wieland
K.
Rothschild
B.
Sundby
S.
Drinkwater
K.
Brander
K.
Ottersen
G.
, et al.  . 
Kruse
G. H.
Drinkwater
K.
Ianelli
J. N.
Link
J. S.
Stram
D. L.
Wespestad
V.
Woodby
D.
Decline and recovery of Atlantic cod (Gadus morhua) stocks throughout the North Atlantic
Resiliency of Gadid Stocks to Fishing and Climate Change
 , 
2008
Fairbanks, Alaska
Alaska Sea Grant College Program
(pg. 
39
-
66
364
Lunn
D. J.
Thomas
A.
Best
N.
Spiegelhalter
D. J.
WinBUGS—a Bayesian modelling framework. Concepts, structure and extensibility
Statistics and Computing
 , 
2000
, vol. 
10
 (pg. 
325
-
337
)
Mace
P. M.
Relationships between common biological reference points used as thresholds and targets of fisheries management strategies
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1994
, vol. 
51
 (pg. 
110
-
122
)
MacKenzie
B. R.
Myers
R. A.
Bowen
K. G.
Spawner–recruit relationships and fish stock carrying capacity in aquatic ecosystems
Marine Ecology Progress Series
 , 
2003
, vol. 
248
 (pg. 
209
-
220
)
Marshall
C. T.
O'Brien
L.
Tomkiewicz
J.
Koster
F. W.
Kraus
G.
Marteinsdottir
G.
Morgan
M. J.
, et al.  . 
Developing alternative indices of reproductive potential for use in fisheries management: case studies for stocks spanning an information gradient
Journal of Northwest Atlantic Fishery Science
 , 
2003
, vol. 
33
 (pg. 
161
-
190
)
Martell
S. J. D.
Pine
W. E.
Walters
C. J.
Parameterizing age-structured models from a fisheries management perspective
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2008
, vol. 
65
 (pg. 
1586
-
1600
)
Mayo
R. K.
Col
L.
Mayo
R. K.
Terceiro
M.
Gulf of Maine cod
Assessment of 19 Northeast Groundfish Stocks through 2004. 2005 Groundfish Assessment Review Meeting (2005 GARM), Northeast Fisheries Science Center, Woods Hole, MA, 15–19 August 2005
 , 
2005
Northeast Fisheries Science Center Reference Document, 05-13
(pg. 
2.153
-
2.184
499
McCarthy
M. A.
Bayesian Methods for Ecology
2007
Cambridge, UK
Cambridge University Press
pg. 
310
 
Metcalfe
D. J.
Fish population structuring in the North Sea: understanding processes and mechanisms from studies of the movements of adults
Journal of Fish Biology
 , 
2006
, vol. 
69(Suppl. C)
 (pg. 
48
-
65
)
Meyer
R.
Millar
R. B.
George
E. I.
State-space models for stock–recruit time series
Bayesian Methods with Applications to Science, Policy, and Official Statistics; Selected Papers from ISBA 2000: the Sixth World Meeting of the International Society for Bayesian Analysis
 , 
2001
Monograph in Official Statistics, Eurostat
(pg. 
361
-
370
)
Michielsens
C. G. J.
McAllister
M. K.
A Bayesian hierarchical analysis of stock–recruit data: quantifying structural and parameter uncertainties
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2004
, vol. 
61
 (pg. 
1032
-
1047
)
Myers
R. A.
Stock and recruitment: generalizations about maximum reproductive rate, density dependence and variability
ICES Journal of Marine Science
 , 
2001
, vol. 
58
 (pg. 
937
-
951
)
Myers
R. A.
Hart
P.
Reynolds
J. D.
Recruitment: understanding density-dependence in fish populations
Handbook of Fish Biology and Fisheries, 1: Fish Biology
 , 
2002
Oxford
Blackwell
(pg. 
123
-
148
413
Myers
R. A.
Barrowman
N. J.
Time series bias in the estimation of density-dependent mortality in stock–recruitment models
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1995
, vol. 
52
 (pg. 
223
-
232
)
Myers
R. A.
Barrowman
N. J.
Hilborn
R.
Kehler
D. G.
Inferring Bayesian priors with limited direct data: applications to risk analysis
North American Journal of Fisheries Management
 , 
2002
, vol. 
22
 (pg. 
351
-
364
)
Myers
R. A.
Bowen
K. G.
Barrowman
N. J.
Maximum reproductive rate of fish at low population sizes
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1999
, vol. 
56
 (pg. 
2404
-
2419
)
Myers
R. A.
Cadigan
N. G.
Density-dependent juvenile mortality in marine demersal fish
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1993
, vol. 
50
 (pg. 
1576
-
1590
)
Myers
R. A.
MacKenzie
B. R.
Bowen
K. G.
Barrowman
N. J.
What is the carrying capacity of fish in the ocean? A meta-analysis of population dynamics of North Atlantic cod
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2001
, vol. 
58
 (pg. 
1464
-
1476
)
Myers
R. A.
Mertz
G.
The limits of exploitation: a precautionary approach
Ecological Applications
 , 
1998
, vol. 
8
 (pg. 
S165
-
S169
)
Myers
R. A.
Mertz
G.
Reducing uncertainty in the biological basis of fisheries management by meta-analysis of data from many populations; a synthesis
Fisheries Research
 , 
1998
, vol. 
37
 (pg. 
51
-
60
)
Neat
F.
Righton
D.
Warm water occupancy by North Sea cod
Proceedings of the Royal Society of London, Series B: Biological Sciences
 , 
2007
, vol. 
274
 (pg. 
789
-
798
)
Needle
G. L.
Recruitment models: diagnosis and prognosis
Reviews in Fish Biology and Fisheries
 , 
2002
, vol. 
11
 (pg. 
95
-
111
)
Nelder
J. A.
Inverse polynomials, a useful group of multi-factor response functions
Biometrics
 , 
1966
, vol. 
22
 (pg. 
128
-
141
)
Nissling
G.
Westin
L.
Egg mortality and hatching rate of Baltic cod (Gadus morhua) in different salinities
Marine Biology
 , 
1991
, vol. 
111
 (pg. 
29
-
32
)
Ntzoufras
I.
Bayesian Modeling using WinBUGS
2009
Hoboken, NJ
Wiley
pg. 
492
 
O'Brien
L.
Munroe
N. J.
Col
L.
Mayo
R. K.
Terceiro
M.
Georges Bank Atlantic cod
Assessment of 19 Northeast Groundfish Stocks through 2004. 2005 Groundfish Assessment Review Meeting (2005 GARM), Northeast Fisheries Science Center, Woods Hole, MA, 15–19 August 2005
 , 
2005
Northeast Fisheries Science Center Reference Document, 05-13
(pg. 
2.2
-
2.29
499
Ottersen
G.
Hjermann
D. Ø.
Stenseth
N. Ch.
Changes in spawning stock structure strengthen the link between climate and recruitment in a heavily fished cod (Gadus morhua) stock
Fisheries Oceanography
 , 
2006
, vol. 
15
 (pg. 
230
-
243
)
Ottersen
G.
Michalsen
K.
Kakken
O.
Ambient temperature and the distribution of Northeast Arctic cod
ICES Journal of Marine Science
 , 
1998
, vol. 
55
 (pg. 
67
-
85
)
Pepin
P.
Orr
D. C.
Anderson
J. T.
Time to hatch and larval size in relation to temperature and egg size in Atlantic cod (Gadus morhua)
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1997
, vol. 
54
 (pg. 
2
-
10
)
Perry
R. I.
Cury
P.
Brander
K.
Jennings
S.
Möllmann
C.
Planque
B.
Sensitivity of marine systems to climate and fishing: concepts, issues and management responses
Journal of Marine Systems
 , 
2009
, vol. 
79
 (pg. 
427
-
435
)
Pinheiro
J. C.
Bates
D. M.
Mixed-Effects Models in S and S-PLUS
2000
New York
Springer
pg. 
528
 
Planque
B.
Frédou
T.
Temperature and the recruitment of Atlantic cod (Gadus morhua)
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1999
, vol. 
56
 (pg. 
2069
-
2077
)
Planque
B.
Fromentin
J-M.
Cury
P.
Drinkwater
K. F.
Jennings
S.
Perry
R. I.
Kifani
S.
How does fishing alter marine populations and ecosystems sensitivity to climate?
Journal of Marine Systems
 , 
2009
, vol. 
79
 (pg. 
403
-
417
)
Power
D.
Healey
B. P.
Murphy
E. F.
Brattey
J.
Dwyer
K.
An assessment of the cod stock in NAFO Divisions 3NO
2005
NAFO Scientific Council Research Document, 05/67
Pyper
B. J.
Peterman
R. M.
Comparison of methods to account for autocorrelation in correlation analyses of fish data
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1998
, vol. 
55
 (pg. 
2127
-
2140
)
Quinn
T. J.
Deriso
R. B.
Quantitative Fish Dynamics
1999
New York
Oxford University Press
pg. 
542
 
Rätz
H. J.
Lloret
J.
Variation in fish condition between Atlantic cod (Gadus morhua) stocks, the effect on their productivity and management implications
Fisheries Research
 , 
2003
, vol. 
60
 (pg. 
369
-
380
)
Ricker
W. E.
Stock and recruitment
Journal of the Fisheries Research Board of Canada
 , 
1954
, vol. 
11
 (pg. 
559
-
623
)
Rijnsdorp
A. D.
Peck
M. A.
Engelhard
G. H.
Möllmann
C.
Pinnegar
J. K.
Resolving the effect of climate change on fish populations
ICES Journal of Marine Science
 , 
2009
, vol. 
66
 (pg. 
1570
-
1583
)
Robichaud
D.
Rose
G. A.
Multiyear homing of Atlantic cod to a spawning ground
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2001
, vol. 
58
 (pg. 
2325
-
2329
)
Robinson
G. K.
That BLUP is a good thing: the estimation of random effects
Statistical Science
 , 
1991
, vol. 
6
 (pg. 
15
-
51
)
Rombough
P. J.
Wood
C. M.
McDonald
D. G.
The effects of temperature on embryonic and larval development
Global Warming: Implications for Freshwater and Marine Fish
 , 
1997
Cambridge, UK
Cambridge University Press
(pg. 
177
-
223
)
Searle
S. R.
Casella
G.
McCulloch
C. E.
Variance Components
1992
New York
John Wiley
pg. 
501
 
Shelton
P. A.
Sinclair
A. F.
Chouinard
G. A.
Mohn
R.
Duplisea
D. E.
Fishing under low productivity conditions is further delaying recovery of Northwest Atlantic cod (Gadus morhua)
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2006
, vol. 
63
 (pg. 
235
-
238
)
Snijders
T. A. B.
Bosker
R. J.
Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling
1999
London
Sage
pg. 
266
 
Spiegelhalter
D.
Best
N.
Bradley
P.
van der Linde
A.
Bayesian measures of model complexity and fit
Journal of the Royal Statistical Society, Series B
 , 
2002
, vol. 
64
 (pg. 
583
-
639
)
Stige
L. C.
Ottersen
G.
Brander
K. M.
Chan
K-S.
Stenseth
N. Ch.
Cod and climate: effect of the North Atlantic Oscillation on recruitment in the North Atlantic
Marine Ecology Progress Series
 , 
2006
, vol. 
325
 (pg. 
227
-
241
)
Su
Z.
Peterman
R. M.
Haeseker
S. L.
Spatial hierarchical Bayesian models for stock–recruitment analysis of pink salmon
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2004
, vol. 
61
 (pg. 
2471
-
2486
)
Sundby
S.
Recruitment of Atlantic cod stocks in relation to temperature and advection of copepod populations
Sarsia
 , 
2000
, vol. 
85
 (pg. 
277
-
298
)
Sundby
S.
Fossum
P.
Feeding conditions of Arcto-Norwegian cod larvae compared with the Rothschild–Osborn theory on small-scale turbulence and plankton contact rates
Journal of Plankton Research
 , 
1990
, vol. 
12
 (pg. 
1153
-
1162
)
Swain
D. P.
Chouinard
G. A.
Predicted extirpation of the dominant demersal fish in a large marine ecosystem: Atlantic cod (Gadus morhua) in the southern Gulf of St Lawrence
Canadian Journal of Fisheries and Aquatic Sciences
 , 
2008
, vol. 
65
 (pg. 
2315
-
2319
)
Trippel
E. A.
Estimation of stock reproductive potential: history and challenges for Canadian Atlantic gadoid stock assessments
Journal of Northwest Atlantic Fishery Science
 , 
1999
, vol. 
25
 (pg. 
61
-
81
)
Vázquez
A.
Cerviño
S.
An assessment of the cod stock in NAFO Division 3M
2002
NAFO Scientific Council Research Document, 02/58
Waldman
J. R.
The importance of comparative studies in stock analysis
Fisheries Research
 , 
1999
, vol. 
43
 (pg. 
237
-
246
)
Walters
C. J.
Bias in the estimation of functional relationships from time series data
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1985
, vol. 
42
 (pg. 
147
-
149
)
Walters
C. J.
Nonstationarity of production relationships in exploited populations
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1987
, vol. 
44
 (pg. 
s156
-
s165
)
Walters
C. J.
Ludwig
D.
Effects of measurement errors on the assessment of stock–recruitment relationships
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1981
, vol. 
38
 (pg. 
704
-
710
)
Waples
R. S.
Gaggiotti
O. E.
What is a population? An empirical evaluation of some genetic methods for identifying the number of gene pools and their degree of connectivity
Molecular Ecology
 , 
2006
, vol. 
15
 (pg. 
1419
-
1439
)
West
B.
Welch
K.
Galecki
A.
Linear Mixed Models: a Practical Guide Using Statistical Software
2006
Boca Raton, FL
Chapman and Hall/CRC
Wikle
C. K.
Hierarchical models in environmental science
International Statistical Review
 , 
2003
, vol. 
71
 (pg. 
181
-
199
)
Worm
B.
Myers
R. A.
Meta-analysis of cod-shrimp interactions reveals top-down control in oceanic food webs
Ecology
 , 
2003
, vol. 
84
 (pg. 
162
-
173
)
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/2.5/uk/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Supplementary data