ICES Journal of Marine Science

Mantzouni, Irene; Sørensen, Helle; O'Hara, Robert B.; MacKenzie, Brian R.

doi: 10.1093/icesjms/fsp291pmid: N/A

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. Open in new tabDownload slide Open in new tabDownload slide 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 (km2) . Year 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 (km2) . Year 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. aShelton et al. (2006). bGoodwin et al. (2006). cEstimated by the authors. Open in new tab 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 (km2) . Year 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 (km2) . Year 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. aShelton et al. (2006). bGoodwin et al. (2006). cEstimated by the authors. Open in new tab 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 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): ⁠. 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 ⁠. 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: ⁠. 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 where yit = log(Rit/Sit), ⁠, and i denotes the stock. For simplicity, we denote 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 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 where yit = log(Rit/Sit), xit = Sit, 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, Keq, 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 yit = α + βxit + ε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: 5a and 5b 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: 6a and 6b 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 7 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) μα: 8 where is the stock-specific model variance and ni 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, 9a and 9b 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 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 and 10b These models describe the (constant in time) divergence between the stock-specific 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 11 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 Combining Equations (9a), (9b), (10a), (11), and (12), the full model is 13 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 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 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); fH1 = 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); fH2 = 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 yit; (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 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); fH1 = 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); fH2 = 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 yit; (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 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. Open in new tab 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 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); fH1 = 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); fH2 = 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 yit; (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 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); fH1 = 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); fH2 = 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 yit; (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 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. Open in new tab 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 Smax). Figure 2. Open in new tabDownload slide 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 where and 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, 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). Figure 3. Open in new tabDownload slide (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. Open in new tab 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. Open in new tab 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. Open in new tabDownload slide (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, (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. Open in new tabDownload slide 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. Open in new tabDownload slide 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 (⁠⁠), 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). Figure 7. Open in new tabDownload slide The fitted stock-level model of beta as a function of H (log habitat size) in (a) the Ricker (corresponding to and (b) the BH (corresponding to ⁠) 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 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): where and 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 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 where and 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 (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. Open in new tabDownload slide (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. Open in new tabDownload slide 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 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 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 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 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 ) Google Scholar Crossref Search ADS PubMed WorldCat 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 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Beverton R. J. H. , Holt S. J. . On the Dynamics of Exploited Fish Populations , 1957 Fishery Investigations, London, Series II, 19 pg. 533 Google Scholar 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 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 Google Scholar Crossref Search ADS WorldCat 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 Google Scholar Cadrin S. X. , Pastoors M. . Precautionary harvest policies and the uncertainty paradox , Fisheries Research , 2008 , vol. 94 (pg. 367 - 372 ) Google Scholar Crossref Search ADS WorldCat 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 Google Scholar 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 Google Scholar Clark J. S. . Models for Ecological Data. An Introduction , 2007 Princeton, NJ Princeton University Press pg. 632 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat Demidenko E. . Mixed Models: Theory and Applications , 2004 Hoboken, NJ Wiley pg. 736 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 Google Scholar 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 Google Scholar 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 ) Google Scholar Crossref Search ADS PubMed WorldCat 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 Google Scholar Gelman A. , Carlin J. B. , Stern H. S. , Rubin D. B. . Bayesian Data Analysis , 2004 2nd edn. London CRC Press Google Scholar Gelman A. , Hill J. . Data Analysis Using Regression and Multilevel/Hierarchical Models , 2007 New York Cambridge University Press pg. 625 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat Gregory D. N. . Climate: a database of temperature and salinity observations for the Northwest Atlantic , 2004 Canadian Science Advisory Secretariat Research Document, 2004/075 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Hilborn R. , Walters C. J. . Quantitative Fisheries Stock Assessment: Choice, Dynamics, and Uncertainty , 1992 New York Chapman and Hall pg. 570 Google Scholar Houde E. D. . Emerging from Hjort's shadow , Journal of Northwest Atlantic Fishery Science , 2008 , vol. 41 (pg. 53 - 70 ) Google Scholar Crossref Search ADS WorldCat ICES. Report of the Study Group on Precautionary Reference Points for Advice on Fishery Management , 2003 ICES Document CM 2003/ACFM: 15 pg. 338 Google Scholar ICES. Spawning and life history information for North Atlantic cod stocks , 2005 ICES Cooperative Research Report, 274 pg. 154 Google Scholar 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 Google Scholar 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 Google Scholar 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 Google Scholar 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 ) Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Kuparinen A. , Merilä J. . The role of fisheries-induced evolution , Science , 2008 , vol. 320 (pg. 47 - 48 ) Google Scholar PubMed OpenURL Placeholder Text WorldCat 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 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC McCarthy M. A. . Bayesian Methods for Ecology , 2007 Cambridge, UK Cambridge University Press pg. 310 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar OpenURL Placeholder Text WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Myers R. A. , Mertz G. . The limits of exploitation: a precautionary approach , Ecological Applications , 1998 , vol. 8 (pg. S165 - S169 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Needle G. L. . Recruitment models: diagnosis and prognosis , Reviews in Fish Biology and Fisheries , 2002 , vol. 11 (pg. 95 - 111 ) Google Scholar Crossref Search ADS WorldCat Nelder J. A. . Inverse polynomials, a useful group of multi-factor response functions , Biometrics , 1966 , vol. 22 (pg. 128 - 141 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Ntzoufras I. . Bayesian Modeling using WinBUGS , 2009 Hoboken, NJ Wiley pg. 492 Google Scholar 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 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Pinheiro J. C. , Bates D. M. . Mixed-Effects Models in S and S-PLUS , 2000 New York Springer pg. 528 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat Quinn T. J. , Deriso R. B. . Quantitative Fish Dynamics , 1999 New York Oxford University Press pg. 542 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat Ricker W. E. . Stock and recruitment , Journal of the Fisheries Research Board of Canada , 1954 , vol. 11 (pg. 559 - 623 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Robinson G. K. . That BLUP is a good thing: the estimation of random effects , Statistical Science , 1991 , vol. 6 (pg. 15 - 51 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Searle S. R. , Casella G. , McCulloch C. E. . Variance Components , 1992 New York John Wiley pg. 501 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat Snijders T. A. B. , Bosker R. J. . Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling , 1999 London Sage pg. 266 Google Scholar 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Sundby S. . Recruitment of Atlantic cod stocks in relation to temperature and advection of copepod populations , Sarsia , 2000 , vol. 85 (pg. 277 - 298 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Vázquez A. , Cerviño S. . An assessment of the cod stock in NAFO Division 3M , 2002 NAFO Scientific Council Research Document, 02/58 Google Scholar Waldman J. R. . The importance of comparative studies in stock analysis , Fisheries Research , 1999 , vol. 43 (pg. 237 - 246 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat Walters C. J. . Nonstationarity of production relationships in exploited populations , Canadian Journal of Fisheries and Aquatic Sciences , 1987 , vol. 44 (pg. s156 - s165 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS PubMed WorldCat West B. , Welch K. , Galecki A. . Linear Mixed Models: a Practical Guide Using Statistical Software , 2006 Boca Raton, FL Chapman and Hall/CRC Google Scholar Wikle C. K. . Hierarchical models in environmental science , International Statistical Review , 2003 , vol. 71 (pg. 181 - 199 ) Google Scholar Crossref Search ADS WorldCat 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 ) Google Scholar Crossref Search ADS WorldCat © 2010 The Author(s) 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. © 2010 The Author(s)

Garrison, Lance P.; Link, Jason S.; Kilduff, D. Patrick; Cieri, Matthew D.; Muffley, Brandon; Vaughan, Douglas S.; Sharov, Alexei; Mahmoudi, Behzad; Latour, Robert J.

doi: 10.1093/icesjms/fsq005pmid: N/A

Garrison, L. P., Link, J. S., Kilduff, D. P., Cieri, M. D., Muffley, B., Vaughan, D. S., Sharov, A., Mahmoudi, B., and Latour, R. J. 2010. An expansion of the MSVPA approach for quantifying predatorprey interactions in exploited fish communities. ICES Journal of Marine Science, 67: 856870.Ecosystem-based fisheries management requires tools to place fish-stock dynamics in the broader context of fishery, predator, and competitive removals. Multispecies virtual population analysis (MSVPA) is an approach to quantifying predatorprey interactions and estimating the rates of predation mortality for exploited fish populations. Here, an extended MSVPA (MSVPA-X) is presented as an alternative to existing MSVPA approaches. Notably, MSVPA-X uses index-tuned VPA methods, applies a more flexible feeding model, and includes an alternative functional feeding response. The MSVPA-X model is applied to a western Atlantic fish community, focusing on Atlantic menhaden and its major fish predators, and a sensitivity analysis of major model parameters is presented. The sensitivity analysis highlights the need for adequate diet sampling. The MSVPA-X represents an improvement over previous approaches by increasing the flexibility to model seasonal and interannual dynamics in the strength of predatorprey interactions. Model results demonstrate that, for menhaden in particular, and forage fish in general, quantifying predation mortality is an important part of effective assessments of forage fish, their predators, and the fisheries of both.

Di Franco, Antonio; Ferruzza, Gaetano; Baiata, Pasquale; Chemello, Renato; Milazzo, Marco

doi: 10.1093/icesjms/fsq007pmid: N/A

Di Franco, A., Ferruzza, G., Baiata, P., Chemello, R., and Milazzo, M. 2010. Can recreational scuba divers alter natural gross sedimentation rate? A case study from a Mediterranean deep cave. ICES Journal of Marine Science, 67: 871874.Submarine caves are environments with features distinguishing them from other littoral habitats but, despite their ecological importance, their response to anthropogenic disturbance has been seldom verified. One potential threat affecting natural communities within caves is represented by recreational scuba diving. Divers' disturbance within marine caves is mainly related to physical contacts and increased sediment resuspension potentially affecting sessile organisms. The aim was to assess the potential effect of recreational divers' frequentation on the natural gross sedimentation rate (GSR) in a Mediterranean deep-water cave. To achieve this, sediment traps were deployed along a scuba trail before, during, and after the peak season for tourist-related diving. No effects of divers' frequentation were evident in terms of alteration of natural GSR, but findings will need to be validated for other caves to asses whether potential disturbance depends on frequentation levels and/or cave characteristics, i.e. sediment grain size or cave depth.

Hufnagl, Marc; Temming, Axel; Siegel, Volker; Tulp, Ingrid; Bolle, Loes

doi: 10.1093/icesjms/fsq003pmid: N/A

Hufnagl, M., Temming, A., Siegel, V., Tulp, I., and Bolle, L. 2010. Estimating total mortality and asymptotic length of Crangon crangon between 1955 and 2006. ICES Journal of Marine Science, 67: 875884.Total mortality (Z, year1) of southern North Sea brown shrimp (Crangon crangon) was determined as Z K, based on the von Bertalanffy lengthgrowth constant (K, year1) and derived from length-based methods. Mortality estimates were based on length frequency distributions obtained from four long-term dataseries (19552006): German Demersal Young Fish Survey, Dutch Demersal Fish Survey, and two German Bycatch series (Bsum and East Frisia). Four methods to estimate and L were evaluated. Highest total mortality (Z 8 year1) was estimated for the early 1990s, and the lowest (Z 4 year1) for the 1960s. Accounting for these differences, a median Z rather than mean values was calculated for the whole series, and the value ranged from 5.74 (Ssentongo and Larkin method), through 5.65 (Beverton and Holt method) and 5.64 (Jones and Zalinge method), to 5.35 (length-converted catch curves). Over the whole period, an increase in and a decrease in the proportion of shrimps >60 mm in the catch was observed, whereas asymptotic length L remained constant (at 79.3 mm total length).

Bierman, Stijn M.; Dickey-Collas, Mark; van Damme, Cindy J. G.; van Overzee, Harrit M. J.; Pennock-Vos, M. G.; Tribuhl, Silja V.; Clausen, Lotte A. W.

doi: 10.1093/icesjms/fsp300pmid: N/A

Bierman, S. M., Dickey-Collas, M., van Damme, C. J. G., van Overzee, H. M. J., Pennock-Vos, M. G., Tribuhl, S. V., and Clausen, L. A. W. 2010. Between-year variability in the mixing of North Sea herring spawning components leads to pronounced variation in the composition of the catch. ICES Journal of Marine Science, 67: 885896.North Sea herring (Clupea harengus) are managed as a single stock, but maintaining a diversity of spawning components is considered important. However, the total catch from each of these components cannot be estimated easily because the components mix during the summer feeding season. The spawning origin of herring is determined from patterns in the microstructure of the otolith core, from samples taken in the central and northern North Sea during summer of 20042007. The annual catch composition of Dutch vessels is determined within a statistical framework that takes account of the spatial patterns in mixing of spawning components and the classification success of the method. Mixing of components varied between years, with steep latitudinal gradients in compositions in some years, resulting in pronounced between-year differences in estimated catch compositions. Differences in lengths-at-age between spawning components, in particular of the 2000 year class, may have caused the observed between-year changes in mixing of components. Our results indicate that estimates of compositions change when assumptions of perfect spatial mixing and perfect classification are relaxed, and can be uncertain in particular as a result of misclassifications, and that it may not be appropriate to assume that ratios between components are constant through time.

Morote, Elvira; Olivar, Mara Pilar; Villate, Fernando; Uriarte, Ibon

doi: 10.1093/icesjms/fsp302pmid: N/A

Morote, E., Olivar, M. P., Villate, F., and Uriarte, I. 2010. A comparison of anchovy (Engraulis encrasicolus) and sardine (Sardina pilchardus) larvae feeding in the Northwest Mediterranean: influence of prey availability and ontogeny. ICES Journal of Marine Science, 67: 897908.The feeding ecology of the larvae of the two most important small pelagic species in the western Mediterranean, anchovy (Engraulis encrasicolus) and sardine (Sardina pilchardus), is investigated. Samples were collected in early summer (anchovy) and autumn (sardine). Composition of the diet and prey selectivity were compared over the course of ontogeny in anchovy (2.714.5 mm standard length, SL) and sardine (5.515.8 mm SL) larvae. Anchovy larvae begin feeding on prey items >150 m at smaller sizes than sardine larvae, and the diets of both species are based mainly on copepod nauplii and postnauplii. Seasonal differences in the composition of the plankton influenced the contributions of prey types to the diets of the two species, e.g. the cladoceran Evadne spp. in anchovy and the tintinnid Codonellopsis sp. in sardine. Although copepod eggs are generally present in the diets of larval clupeoids, they were not major components of the diets of the species considered here. Despite morphological similarity, selection patterns were different between the species and changed through ontogeny.

Spitz, Jrme; Mourocq, Emeline; Schoen, Valrie; Ridoux, Vincent

doi: 10.1093/icesjms/fsq008pmid: N/A

Spitz, J., Mourocq, E., Schoen, V., and Ridoux, V. 2010. Proximate composition and energy content of forage species from the Bay of Biscay: high- or low-quality food? ICES Journal of Marine Science, 67: 909915.Collapses of high-energy dense concentrations of prey species induce negative effects on populations of top predators. Knowledge of prey quality appears to be crucial in ecosystem modelling and management. The aim here was to provide baseline data of forage species quality in the Bay of Biscay. Proximate composition (water, ash, protein, and lipid) and energy content have been determined to assess the quality of 78 species, including jellyfish, crustaceans, cephalopods, cartilaginous fish, and bony fish. Results show broad variations between species, with energy densities ranging from 2 to 10 kJ g1. Lipids are the most structuring component and largely determined prey quality, and prey species are not necessarily interchangeable for the fulfilment of a predator's energy and food requirements. In ecosystem models, therefore, multispecies compartments of forage organisms would ideally be constituted using prey species of equivalent quality and hence of equivalent benefit to top predators.

De Oliveira, Jos A. A.; Darby, Chris D.; Roel, Beatriz A.

doi: 10.1093/icesjms/fsp290pmid: N/A

De Oliveira, J. A. A., Darby, C. D., and Roel, B. A. 2010. A linked separableADAPT VPA assessment model for western horse mackerel (Trachurus trachurus), accounting for realized fecundity as a function of fish weight. ICES Journal of Marine Science, 67: 916930.The western horse mackerel stock covers a large area of the Northeast Atlantic and is characterized by the sporadic occurrence of strong year classes. The only fishery-independent data available for its assessments are the estimates of total annual egg production from triennial egg surveys. Horse mackerel are indeterminate spawners, so a direct conversion of total annual egg production to spawning-stock biomass using prespawning-season estimates of fecundity is not viable. There is also evidence that potential fecundity per kg spawning female increases with fish weight. A linked separableADAPT VPA (SAD) model was developed that combines data on total catch, catch-at-age, and total annual egg production with data on potential and realized fecundity to provide population estimates for western horse mackerel. The model accounts for potential fecundity as a function of fish weight and for the development of a targeted fishery on the strong 1982 year class. Simulations confirm that the SAD model is able to reproduce population estimates without bias under a range of scenarios, except where there is a trend in realized fecundity. This underscores the need for improved information on realized fecundity, or alternatively the need to develop management plans that are robust to this source of uncertainty.

Farrell, Edward D.; Mariani, Stefano; Clarke, Maurice W.

doi: 10.1093/icesjms/fsp295pmid: N/A

Farrell, E. D., Mariani, S., and Clarke, M. W. 2010. Age and growth estimates for the starry smoothhound (Mustelus asterias) in the Northeast Atlantic Ocean. ICES Journal of Marine Science, 67: 931939.Age, growth, and longevity were estimated for the starry smoothhound (Mustelus asterias), based on the interpretation of sectioned vertebrae of 106 males (35104 cm stretch total length, STL) and 114 females (44112 cm STL). Growth curves were fitted to length-at-age data using von Bertalanffy and Gompertz models. The three-parameter von Bertalanffy growth function (VBGF) provided the best statistical fit to the male data (L 103.7 cm STL, L0 38.1 cm STL, k 0.195 year1). The Gompertz growth function provided the best statistical fit to the female data, although the estimated parameters were biologically unreasonable. Therefore, the three-parameter VBGF was also accepted for females (L 123.5 cm STL, L0 34.9 cm STL, k 0.146 year1). Longevity was estimated to be 13 and 18.3 years for males and females, respectively. The lengthweight relationship is also presented for 304 male and 424 female M. asterias. The von Bertalanffy model was fitted to the weight-at-age data. These estimates can form the basis of future work on the assessment and management of this species.

Zavalaga, Carlos B.; Halls, Joanne; Dell'Omo, Giacomo

doi: 10.1093/icesjms/fsp298pmid: N/A

Zavalaga, C. B., Halls, J., and Dell'Omo, G. 2010. Marine habitat use of Peruvian boobies: a geographic and oceanographic comparison between inshore and offshore islands. ICES Journal of Marine Science, 67: 940951.Foraging areas of Peruvian boobies (Sula variegata) from Isla Lobos de Tierra (LT, inshore) and Isla Lobos de Afuera (LA, offshore) were overlaid with concurrent data on wind patterns and remotely sensed hydrographic features to identify the main abiotic factors that affect booby distribution and to compare habitat use between birds from inshore and offshore islands. Birds used across winds and across headwinds when commuting to their feeding grounds and across tailwinds when returning to their colonies. This strategy increased flight speeds by 3853 in relation to headwinds, probably resulting in a reduction in energy costs when birds returned to their nests with heavy loads of food. Sea surface temperature of the foraged and the available areas were similar around LT and LA. Boobies from LT fed over the continental shelf and alongshore in areas with chlorophyll a (Chl a) concentration higher than expected from a random distribution. Conversely, birds from LA never came inshore, feeding over the shelf break and in waters with concentrations of Chl a similar to those in the available area. Remote-sensing measurements of primary productivity predicted the distribution of Peruvian boobies only for LT. Plumes of enriched water alongshore are formed during the austral summer, favouring the aggregation of anchoveta (Engraulis ringens), the main prey of Peruvian boobies. Boobies from LA may have potentially foraged inshore in an area also used by birds from LT, but it is likely that competition for food and energy constraints to raise three-chick broods restricted foraging range to waters far from the mainland.