## Abstract

Understanding how temperature affects cod (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.

*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 (*F*_{MSY}). 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).

Stock | Area | SPR_{F}_{=0} (kg) | Habitat size (km^{2}) | Year range | Reference for stock data |
---|---|---|---|---|---|

codgb | Georges Bank (5Z) | 23.8^{a} | 85 306 | 1978–2004 | O'Brien et al. (2005) |

codgom | Gulf of Maine (5Y) | 27.9^{a} | 47 629 | 1982–2003 | Mayo and Col (2005) |

cod4x | Western Scotian Shelf (4X) | 14.7^{a} | 64 331 | 1983–2000 | Clark et al. (2002) |

cod4vsw | Eastern Scotian Shelf (4VsW) | 11.7^{a} | 91 656 | 1970–2001 | Fanning et al. (2003) |

cod4tvn* | Southern Gulf of St Lawrence (4TVn) | 7.0^{a} | 80 338 | 1950–2004 | Chouinard et al. (2006) |

cod3pn4rs* | Northern Gulf of St Lawrence (3Pn4RS) | 4.1^{a} | 89 734 | 1974–2003 | Fréchet et al. (2005) |

cod3ps* | Southern Newfoundland (3Ps) | 5.9^{a} | 55 509 | 1977–2002 | Brattey et al. (2004) |

cod3no* | Grand Banks (3NO) | 6.3^{a} | 114 967 | 1959–2004 | Power et al. (2005) |

cod3m | Flemish Cap (3M) | 9.8^{a} | 17 398 | 1972–1993 | Vázquez and Cerviño (2002) |

cod2j3kl* | Northern Newfoundland (2J3KL) | 7.4^{a} | 263 893 | 1962–1989 | Bishop et al. (1993) |

cod-iceg* | Iceland (Va) | 18.9^{b} | 63 158 | 1955–2003 | ICES DB 2006 |

cod-farp* | Faroe Plateau (Vb) | 11.6^{b} | 12 883 | 1961–2004 | ICES DB 2006 |

cod-arct* | Northeast Arctic (I, II) | 12.1^{b} | 858 979 | 1946–2003 | ICES DB 2006 |

cod-coas* | Norwegian Coastal (IIa) | 6.2^{c} | 234 615 | 1984–2004 | ICES DB 2006 |

cod-2224* | Western Baltic (IIId—west) | 5.3^{b} | 41 823 | 1970–2005 | ICES DB 2006 |

cod-2532* | Eastern Baltic (IIId—east) | 3.3^{c} | 131 019 | 1966–2003 | ICES DB 2006 |

cod-kat | Kattegat (IIIa—east) | 7.8^{b} | 19 926 | 1971–2004 | ICES DB 2006 |

cod-347d | North Sea (IIIa, IV, VIId) | 18.2^{b} | 418 821 | 1963–2004 | ICES (2007b) |

codviia | Irish Sea (VIIa) | 12.7^{b} | 23 940 | 1968–2005 | ICES (2006) |

cod-7ek | Celtic Sea (VIIe–k) | 19.9^{b} | 177 620 | 1971–2005 | ICES DB 2006 |

codvia | West of Scotland (VIa) | 12.9^{b} | 81 016 | 1978–2004 | ICES (2006) |

Stock | Area | SPR_{F}_{=0} (kg) | Habitat size (km^{2}) | Year range | Reference for stock data |
---|---|---|---|---|---|

codgb | Georges Bank (5Z) | 23.8^{a} | 85 306 | 1978–2004 | O'Brien et al. (2005) |

codgom | Gulf of Maine (5Y) | 27.9^{a} | 47 629 | 1982–2003 | Mayo and Col (2005) |

cod4x | Western Scotian Shelf (4X) | 14.7^{a} | 64 331 | 1983–2000 | Clark et al. (2002) |

cod4vsw | Eastern Scotian Shelf (4VsW) | 11.7^{a} | 91 656 | 1970–2001 | Fanning et al. (2003) |

cod4tvn* | Southern Gulf of St Lawrence (4TVn) | 7.0^{a} | 80 338 | 1950–2004 | Chouinard et al. (2006) |

cod3pn4rs* | Northern Gulf of St Lawrence (3Pn4RS) | 4.1^{a} | 89 734 | 1974–2003 | Fréchet et al. (2005) |

cod3ps* | Southern Newfoundland (3Ps) | 5.9^{a} | 55 509 | 1977–2002 | Brattey et al. (2004) |

cod3no* | Grand Banks (3NO) | 6.3^{a} | 114 967 | 1959–2004 | Power et al. (2005) |

cod3m | Flemish Cap (3M) | 9.8^{a} | 17 398 | 1972–1993 | Vázquez and Cerviño (2002) |

cod2j3kl* | Northern Newfoundland (2J3KL) | 7.4^{a} | 263 893 | 1962–1989 | Bishop et al. (1993) |

cod-iceg* | Iceland (Va) | 18.9^{b} | 63 158 | 1955–2003 | ICES DB 2006 |

cod-farp* | Faroe Plateau (Vb) | 11.6^{b} | 12 883 | 1961–2004 | ICES DB 2006 |

cod-arct* | Northeast Arctic (I, II) | 12.1^{b} | 858 979 | 1946–2003 | ICES DB 2006 |

cod-coas* | Norwegian Coastal (IIa) | 6.2^{c} | 234 615 | 1984–2004 | ICES DB 2006 |

cod-2224* | Western Baltic (IIId—west) | 5.3^{b} | 41 823 | 1970–2005 | ICES DB 2006 |

cod-2532* | Eastern Baltic (IIId—east) | 3.3^{c} | 131 019 | 1966–2003 | ICES DB 2006 |

cod-kat | Kattegat (IIIa—east) | 7.8^{b} | 19 926 | 1971–2004 | ICES DB 2006 |

cod-347d | North Sea (IIIa, IV, VIId) | 18.2^{b} | 418 821 | 1963–2004 | ICES (2007b) |

codviia | Irish Sea (VIIa) | 12.7^{b} | 23 940 | 1968–2005 | ICES (2006) |

cod-7ek | Celtic Sea (VIIe–k) | 19.9^{b} | 177 620 | 1971–2005 | ICES DB 2006 |

codvia | West of Scotland (VIa) | 12.9^{b} | 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. SPR_{F =0} refers to spawners produced per recruit (SPR) in the absence of fishing mortality (*F*). The areas are shown in Figure 1.

^{a}Shelton *et al*. (2006).

^{b}Goodwin *et al*. (2006).

^{c}Estimated 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 SPR_{F}_{=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 SPR_{F}_{=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):

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

*A*

^{RIC}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 SPR

_{F}

_{=0}. Therefore, as will be shown in detail below,

*A*

^{RIC}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 SPR

_{F}

_{=0}components.

Parameter *B* is related to carrying capacity because 1/*B*) equals spawner biomass when recruitment reaches the maximum (*S*_{max}, 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 *R*_{max} and *S*_{max} (Brander and Mohn, 2004a, b): . According to that formula, the parameters of the Ricker model can be expressed as *A*^{RIC}= exp(1)(*R*_{max}/*S*_{max}) and *beta* = 1/*B* = *S*_{max}. Therefore, *K*_{max}, representing the maximum number of recruits produced by the maximum number of spawners sustained by the ecosystem (i.e. *R*_{max}), is estimated as . *K*_{eq} 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: . The Ricker model can be linearized by natural log-transformation and assuming lognormal errors: log(*R*_{t}/*S*_{t}) = log(*A*^{RIC}) − *BS*_{t} + ε_{t}. To simplify notation, the Ricker model for stock *i* is written as

*y*

_{it}= log(

*R*

_{it}/

*S*

_{it}), , and

*i*denotes the stock. For simplicity, we denote as

*alpha*, understanding that

*alpha*= log(

*A*

^{RIC}). 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:

Ricker and BH SR models display similar behaviour at the limit of low*S*(i.e. compensation); therefore, parameter

*A*(denoted as

*A*

^{BH}for the BH model) has common interpretation and estimation in both models (Myers

*et al*., 1999). Parameter

*K*

^{BH}has the same dimensions as

*S*and can be interpreted as the “threshold biomass” resulting in half the maximum recruitment, which equals

*A*

^{BH}×

*K*

^{BH}.

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:

where*y*

_{it}= log(

*R*

_{it}/

*S*

_{it}),

*x*

_{it}=

*S*

_{it}, and

*i*denotes the stock. In this form, the model is linear to , which has the same interpretation as of the Ricker model in Equation (2) and will also be denoted as

*alpha*, and is comparable with the Ricker model , because it corresponds to

*S*, resulting in half the maximum recruitment given by . For simplicity, 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,

*K*

_{eq}, for the BH model as . 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 and , respectively. However, the same approaches apply to the BH model and its parameters and . 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 *y*_{it} = α + *β**x*_{it} + ε_{it}, with

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:

and where*μ*

_{α}and (or

*μ*

_{β}and ) 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 and 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: and where and 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: 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

It is assumed that the residuals associated with each stock are independent and that the residuals and random effects are independent of each other. The parameter 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 , corresponding to in Equation (2) and its mean across stocks (the fixed effect) *μ*_{α}:

*n*

_{i}the number of observations for stock

*i*(Gelman and Hill, 2007). Pooling is stronger when is small and for stocks with fewer observations and greater variability (). The stock-specific weight is defined as the reliability of the corresponding stock mean, and the ratio of the two weights is the ratio of the true variance to the error variance . As the true values are unknown, they are substituted by their estimates in Equation (8). The above estimate 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,

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:

Note that the temperature time-series in Equations (9a) and (9b) were centred to the overall (across stocks) mean (i.e. average of all *T*_{it} observations) to remove the correlation between the first- and the second-order term estimates. Therefore, the intercepts *c*_{oi} and *d*_{oi} 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):

*alpha*or

*beta*estimate and the across-stocks grand mean or , respectively. The coefficients describing the temperature effect on

*alpha*and

*beta*, (, ) and (, ), respectively, are assumed to be common across stocks. This assumption can be relaxed by imposing stock-level models on the coefficients. For instance, in Equation (9a) can be substituted by the stock-specific terms . Consequently, the Ricker SR model, updated to incorporate temperature effects on

*alpha*and

*beta*, is 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, *d*_{oi}, 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

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

with random effects and 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 *y*_{i} on *x*_{i} 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 *H*_{i}.

### 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 + 2*p*_{d}. 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 *p*_{d}, 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.

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 | 1 662.5 (REML) | −825.2 (REML) [6] | – | – | ||

MR.2 | As in MR.1 | 1 661.6 (REML); | −825.8 (REML) [5] | Random effects are independent; ρ = 0 (MR.1 vs. MR.2) | 0.30 (1.1) | |

MR.3 | 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 | , | 1 576.3 (ML) | −762.2 (ML) [26] | Linear dependence of β_{i} on habitat size (H); f_{H}_{1} = 0 (MR3.HL vs. MR.3) | 0.03 (4.6) | |

MR3.HQ | As in MR3.HL | 1 561.2 (ML) | −753.6 (ML) [27] | Quadratic dependence of β_{i} on habitat size (H); f_{H}_{2} = 0 (MR3.HQ vs. MR3.HL) | <0.01 (17.2) | |

MR3.HQ.TAL | , , | 1 562.4 (ML) | −753.2(ML) [28] | Linear dependence of α_{i} on temperature (T); (MR3.HQ.TAL vs. MR3.HQ.TAQ) | <0.01 (14.3) | |

MR3.HQ.TAQ | 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); (MR3.HQ.TAQ vs. MR3.HQ) | <0.01 (15.0) | |

MR3.HQ.TBL | , | 1 562.6 (ML) | −753.3 (ML) [28] | Linear dependence of β_{i} on temperature (T); (MR3.HQ.TBL vs. MR3.HQ) | 0.46 (0.6) | |

MR3.HQ.TBQ | As in MR3.HQ.TBL | 1 562.2 (ML) | −752.3 (ML) [29] | Quadratic dependence of β_{i} on temperature (T); (MR3.HQ.TBQ vs. MR3.HQ) | 0.23 (2.9) | |

MR3.HQ.TAQ.RS | , ,, | 1 609.4 (REML) | −773.7 (REML) [31] | Across-stocks differences in the influence of T on alpha; (MR3.HQ.TAQ.RS vs. MR3.HQ.TAQ) | 0.5 (1.4) | |

MR3.HQ.AC | As in MR3.HQ | 1 398.3 (ML) | −651.2 (ML) [48] | – | – | |

MR3.HQ.TAQ.AC | 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 y_{it}; (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 y_{it} | 665.3 (ML) | −317.7 (ML) [15] | – | – | |

MR3.HL.TAL.r | As in MR3.HQ.TAQ | 660.7 (ML) | −314.4 (ML) [16] | Linear dependence of α_{i} on temperature (T); (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 | 1 662.5 (REML) | −825.2 (REML) [6] | – | – | ||

MR.2 | As in MR.1 | 1 661.6 (REML); | −825.8 (REML) [5] | Random effects are independent; ρ = 0 (MR.1 vs. MR.2) | 0.30 (1.1) | |

MR.3 | 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 | , | 1 576.3 (ML) | −762.2 (ML) [26] | Linear dependence of β_{i} on habitat size (H); f_{H}_{1} = 0 (MR3.HL vs. MR.3) | 0.03 (4.6) | |

MR3.HQ | As in MR3.HL | 1 561.2 (ML) | −753.6 (ML) [27] | Quadratic dependence of β_{i} on habitat size (H); f_{H}_{2} = 0 (MR3.HQ vs. MR3.HL) | <0.01 (17.2) | |

MR3.HQ.TAL | , , | 1 562.4 (ML) | −753.2(ML) [28] | Linear dependence of α_{i} on temperature (T); (MR3.HQ.TAL vs. MR3.HQ.TAQ) | <0.01 (14.3) | |

MR3.HQ.TAQ | 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); (MR3.HQ.TAQ vs. MR3.HQ) | <0.01 (15.0) | |

MR3.HQ.TBL | , | 1 562.6 (ML) | −753.3 (ML) [28] | Linear dependence of β_{i} on temperature (T); (MR3.HQ.TBL vs. MR3.HQ) | 0.46 (0.6) | |

MR3.HQ.TBQ | As in MR3.HQ.TBL | 1 562.2 (ML) | −752.3 (ML) [29] | Quadratic dependence of β_{i} on temperature (T); (MR3.HQ.TBQ vs. MR3.HQ) | 0.23 (2.9) | |

MR3.HQ.TAQ.RS | , ,, | 1 609.4 (REML) | −773.7 (REML) [31] | Across-stocks differences in the influence of T on alpha; (MR3.HQ.TAQ.RS vs. MR3.HQ.TAQ) | 0.5 (1.4) | |

MR3.HQ.AC | As in MR3.HQ | 1 398.3 (ML) | −651.2 (ML) [48] | – | – | |

MR3.HQ.TAQ.AC | 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 y_{it}; (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 y_{it} | 665.3 (ML) | −317.7 (ML) [15] | – | – | |

MR3.HL.TAL.r | As in MR3.HQ.TAQ | 660.7 (ML) | −314.4 (ML) [16] | Linear dependence of α_{i} on temperature (T); (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 ( and ) 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 *S*_{max}).

Although no significant differences were found among stocks for the *T*-related slopes (*c*_{T}_{1}_{i} and *c*_{T}_{2}_{i}), 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

### Autocorrelation

Autocorrelation at lag 1 in log(*R*_{t}/*S*_{t}) 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(*R*_{t}_{−1}/*S*_{t}_{−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, and 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, and (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 *alpha*–*T* 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).

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 *c*_{T}_{1}_{i} and/or *c*_{T}_{2}_{i} 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 −*c*_{T}_{1}_{i}/2*c*_{T}_{2}_{i}, 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}/d*T* = *c*_{T}_{1}_{i} + 2*T*_{p}*c*_{T}_{2}_{i}, where *T*_{p} 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).

The stock-specific intercepts *c*_{oi} 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, (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).

### 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)].

### The dependence of Ricker *beta* on habitat size

The spawner biomass producing maximum recruitment, according to the Ricker SR model (), 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, . 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 which has units of *S*).

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

*V*is the finite-sample variance operator across stocks. The numerator represents the average variance in , i.e. in the average variance of 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

*r*

^{2}(Gelman and Hill, 2007). The intermediate value of

*r*

^{2}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 km

^{2}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 *K*_{max} and *K*_{eq}, 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 *K*_{max} 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 *K*_{eq} 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}/d*T*, 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 *c*_{oi} (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 km^{2} (Figure 7d) and *c*_{oi} (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

where andThe two SR models, Ricker and BH, provided similar estimates for the *alpha*-related terms *c*_{oi} (Figure 5), *c*_{T}_{1}_{i}, and *c*_{T}_{2}_{i} (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 (*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 *r*^{2} 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 *K*_{eq} provided by the two SR models (Table S1) and to a lesser extent in *K*_{max} (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.

### 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.

## 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 in the Ricker, and 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 SPR_{F}_{=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 SPR_{F}_{=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 SPR_{F}_{=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 SPR_{F}_{=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 SPR_{F}_{=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 *alpha*–*T* 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 SPR_{F}_{=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* (*K*_{eq}) or as the *S* producing the maximum replacement spawner biomass (*K*_{max}), 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 *K*_{eq} obtained from the Ricker model, so temperature effects are more severe for *K*. For example, if mean *T* increased by 3°C, *K*_{max} 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* (=*S*_{max}), on a per-km^{2} 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 (*r*^{2}) 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

## 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

*Gadus morhua*L.)

*Gadus morhua*) using a comparative approach

*Gadus morhua*)

*Gadus morhua*)

*Gadus morhua*) stock in NAFO Subdivision 3Ps in October 2004

*Gadus morhua*) off Newfoundland

*Gadus morhua*) to future climate change

*Gadus morhua*) stocks and limits to growth production

*Gadus morhua*) during early life

*Gadus morhua*)

*Gadus morhua*) stocks throughout the North Atlantic

*Gadus morhua*) in different salinities

*Gadus morhua*) stock

*Gadus morhua*)

*Gadus morhua*)

*Gadus morhua*) stocks, the effect on their productivity and management implications

*Gadus morhua*)

*Gadus morhua*) in the southern Gulf of St Lawrence