Running Title: Thermal Effects on Cod Recruitment

Productivity responses of a wide-spread marine piscivore, Gadus morhua, to oceanic thermal extremes and trends

Running title: Thermal effects on cod recruitment

Irene Mantzouni1* and Brian R. MacKenzie1, 2, 3

Supplementary Material

1National Institute for Aquatic Resources

TechnicalUniversity of Denmark (DTU-Aqua)

Jægersborg Allé 1

CharlottenlundCastle

DK-2920 Charlottenlund

Denmark

2Department of Marine Ecology, University of Aarhus,

c/o DTU-Aqua,

Jægersborg Allé 1

CharlottenlundCastle

DK-2920 Charlottenlund

Denmark

3Center for Macroecology, Evolution and Climate

Department of Biology

University of Copenhagen

Universitetsparken 15

DK-2100 Copenhagen

Denmark

*Corresponding author (IM): Tel: +45 33963422;

Datasets

Cod spawner stock and recruitment data

We compiled a database including population time-series for the 21 major cod stocks in the N. Atlantic (Supplementary Table S1, Supplementary Figure S1). Population data include time-series of spawner stock biomass (SSB; in 000’s tons) and recruitment (REC; in 000’s individuals). These numbers are estimated from sequential population analysis (SPA) standardized, in most cases, with fisheries independent (such as research trawl survey) data and were extracted from published stock assessment reports (Supplementary Table S1).

Temperature

Since recruitment success in many stocks and years is determined during the egg-larval-pelagic juvenile phase (Brander 2003), spring temperature time-series averaged over the area occupied by each stock at the surface (0-100m) layer were used. Temperature was estimated across the entire fisheries statistical areas of the stocks, rather than within the spawning locations, except for the cases mentioned below. The main reason for this choice is that we are aiming at studying thermal effects on the survival of the early stages, rather than only eggs, operating through a broader range of mechanisms. The exact distribution of these stages is only partly known, and can vary inter-annually, depending on oceanographic processes rather than active habitat selection. Also, apart from the ambient conditions altering the physiological rates, temperature can affect cod also indirectly, e.g., by influencing the trophic interactions with prey and predators (Houde 2008; Rijnsdorp et al. 2009). Moreover, substantial spatial correlation has been found among sub-areas comprising large statistical regions (Planque & Frédou 1999; MacKenzie & Schiedek 2007), and thus no considerable bias is expected due to this choice in most cases.

Some special considerations apply to certain NE and NW areas. For eastern Baltic cod stock (ICES Subdivisions 25-32) we used temperature estimates in Subdivisions 25-29, since Subdivisions 30-32 are unfavorable for cod reproduction due to low salinity (Nissling & Westin 1991). For Barents Sea cod, given that stock distribution can be limited by cold waters (Ottersen et al. 1998), we estimated temperature in the area south of 78oN. Low water temperature can also limit the distribution of Icelandic cod to the southern part (ICES 2005). Therefore we used temperature estimates applying to the region south of 62 oN in ICES subdivision Va. We also applied spatial restrictions in the 4 NW Atlantic areas (Flemish Cap-3M, Grand Bank-3NO, W and E Scotian Shelf- 4VsW and 4X) mostly affected by the Gulf Stream, by excluding temperature observations southern that 42o. For the 3M area, only temperature within the Flemish Cap was used.

The temperature time-series for the NE Atlantic stocks were provided by the ICES (International Council for the Exploration of the Sea) data centre ( For the NW Atlantic, the datasets were extracted from the DFO (Fisheries and Oceans Canada) Hydrographic Climate database (Gregory 2004). The temperature datasets were explored in terms of consistency in spatial coverage throughout years and results were satisfactory for most areas.

Methods

Effect Sizes

1. Risk Ratio (RR)

The Risk Ratio was employed to test the null hypothesis that the probabilities (risks) of “successful” year-classes were equal during extremely warm (“Exposed” or “Hot” group; TT75th%ile) and cold (“Control” or “Cool” group; TT25th%ile) seasons. “Successful” year-classes are considered those exceeding the corresponding time-series mean, and vice versa for the “failed” ones (for Ricker residuals, successful and failed events correspond to the positive and negative values, respectively).

In other words, the aim is to investigate whether the chances of higher recruitment survival [log(REC/SSB) or Ricker model residuals] differ between colder and warmer seasons. Therefore, for each stock, we count the number of observations falling in each cell of the following fourfold table:

From these frequencies we estimate the RRi for each stock i:

where and . Thus, if RR is above 1, there is higher probability (risk) of “successful” events during ”Hot” seasons, and vice versa if RR1.

For the statistical analysis the natural logarithm of RRi is used, because of its better statistical properties (Cooper & Hedges 1994: 248). The associated sampling variance is:

2. Hedges’g (HG)

HG belongs to a broader category of t-test related metrics, used to quantify the difference between the means of two groups, subjected to different treatments (Hedges & Olkin 1985). Hence, HG is used to compare average recruitment survival between “Cool”and ”Hot”seasons, testing the null hypothesis that there is no substantial difference. The stock- specific HGi estimator is:

where () and () is the mean (sample size) in the warmer and colder of seasons, respectively. is the correction for small sample size and approaches 1 as the number of observations increases (Hedges & Olkin 1985). The is the pooled standard deviation between the two groups of observations:

where and is the standard deviation of the observations during the warmer and colder seasons, respectively. The variance of HGi is given by:

A negative estimate of HG implies that mean recruitment survival is lower during relatively warm seasons, and vice versa for HG0.The HG can be transformed into Cohen’s d (CD), an effect-size based also on the standardized difference between the means (Cohen 1988, Cooper & Hedges 1994: 239). The CD has useful interpretations, such as the percent of non-overlap between the distributions of the recruitment survival observations during the warmer and colder season, respectively (Cohen 1988).

3. Fisher’s z correlation coefficient (FZ)

FZ is simply the Fisher’s z-transform of the correlation (ri) between recruitment survival and temperature. T his transformation is usually preferred, since it is nearly normally distributed (Cooper & Hedges 1994:240):

The effective sample size , corrected for autocorrelation at lag 1 within the stock and temperature time-series, was estimated using the “modified Chelton” method (Pyper & Peterman 1998):

,

where ni is the number of observations and wi,P, wi,Tthe autocorrelation at lag 1 in the stock and temperature time-series, respectively. The following formula was used to estimate the variance, in order to avoid inflation due to low sample size (Stuart & Ord 1987: 533):

Random-Effects Meta-Analyses

The above metrics (Risk Ratio-RR, Hedges’ g-HG, Fisher’s z correlation coefficient- FZ) were analyzed separately for the Cold and Warm stocks using random-effects meta-analyses. The method uses the stock-specific estimates to produce a weighted, average metric, representing the across-stocks effect within each group (Cold or Warm). The stock-specific proportional contribution (weight) to the estimation of the mean effect-size depends on both the sampling error (estimated using the above presented formulas) and the random-effects variance. The latter represents true differences among the individual metrics, due to variability in the underlying conditions (Cooper & Hedges 1994: 316). In other words, there is a (normal) distribution of metrics across the stocks, with each stock-specific estimate being drawn from that distribution. This approach is of hierarchical nature, since it considers two levels of variability, both within (sampling variance) and across (random-effects variance) stocks (Cooper & Hedges 1994: 366). Thus, this method provides more conservative estimates of significance compared to a fixed-effects model, considering only the former source of variability (Cooper & Hedges 1994: 275). The choice of random-effects meta-analysis is more justified in our study where we are combining estimates obtained for various stocks, distributed across N. Atlantic and thus occupying a number of ecosystems that differ in many biotic and abiotic aspects (Osenberg et al. 1999; Worm & Myers 2003; Lilly et al. 2008). The meta-analyses were conducted using the MS Excel add-in tool, MIX 1.7 (Bax et al. 2008).

Auto-correlation (AC)

An important aspect that should be considered for the first two sets of meta-analyses is the presence of auto-correlation (AC) in the stock time-series. In order to identify this effect, we plotted the empirical auto-correlation function referring to the entire time-series of each stock. Positive AC at lag 1 was found significant in a number of cases as illustrated in the Supplementary Table S1. The issue is readily dealt with for Fisher’s z, by using the appropriate formulas for the estimation of the effective sample size, as presented above.

In order to reduce AC effects within the stock-specific “Cool”and “Hot” groups for the Risk Ratio estimation, when present, a possible approach is to delete the appropriate observations so that consecutive events are eliminated from the groups (e.g., if the recruitment observations in the “Exposed” group of a given stock correspond to years 1990, 1991 and 1992, the 1992 observation was omitted, so that there is at least a two-year interval between the events). RRi is based on the number of successful (Gi or Ci) and failed (gi or ci) events, and hence the value of the ratio is not sensitive to which observations are excluded.

Regarding Hedges’ g, this metric is sensitive to the values of the observations included in the groups under comparison [“Hot” (or “Exposed”) and “Cool” (or“Control”) seasons]. Therefore, the method used in order to eliminate possible AC within the groups in RR estimation, could not be applied to this analysis. The Wilcoxon test described below is an alternative approach used in order to accommodate possible AC effects between the observation groups, in this case.

Wilcoxon signed rank test

This nonparametric test is an alternative to the Student’s paired t-test, applicable to cases when means are compared between two correlated groups (Zar 1999). In the present case, it is used as an alternative to the Hedges’ g, presented above, in order to consider also the auto-correlation in the recruitment survival observations (discussed later), resulting in correlated groups of observations. The test is based on a similar metric:

where and is the mean recruitment survival in the ”Hot and “Cool”seasons, respectively. The null hypothesis tested is that the median of the estimates distribution is 0, i.e., that there is no difference in recruitment survival under the two temperature regimes. The alternative, one-sided hypotheses were that the median is either negative or positive. For these tests, the exact probability was used, since the number of stocks within either stock group is 25 (Zar 1999).

References

Bax, L., Yu, L.M., Ikeda, N., Tsuruta, H., Moons, K.G.M. 2008 MIX: comprehensive free software for meta-analysis of causal research data Version 17

Bishop, C. A., Murphy, E. F., Davis, M. B., Baird, J. W.Rose, G. A. 1993 An assessment of the cod stock in NAFO Divisions 2J+3KL. NAFO Scientific Council Research Document, 93/86.

Brander, K.M. 2003 What kinds of fish stock predictions do we need and what kinds of information will help us to make better predictions? Scientia Marina,67, 21–33.

Brattey, J., Cadigan, N. G., Healey, B. P., Lilly, G. R., Murphy, E. F., Shelton, P. A.Mahé, J-C. 2004 An assessment of the cod (Gadus morhua) stock in NAFO Subdivision 3Ps in October 2004. DFO Canadian Science Advisory Secretariat Reseach Document, 2004/083.

Chouinard, G. A., Currie, L., Poirier, G. A., Hurlbut, T., Daigle, D. & Savoie, L. 2006 Assessment of the southern Gulf of St Lawrence cod stock, February 2006. Canadian Science Advisory Secretariat Research Document, 2006/006.

Clark, D. S., Gavaris, S. Hinze, J. M. 2002 Assessment of cod in Division 4X in 2002. Canadian Science Advisory Secretariat Research Document, 2002/105.

Cohen, J. 1988. Statistical power analysis for the behavioral sciences (2nd ed). HillsdaleNJ: Lawrence Earlbaum Associates

Cooper, H. Hedges, L.V. (Eds.) 1994 The handbook of research synthesis. Russell Sage Foundation, New York, 573 pp.

Fanning, L. P., Mohn, R. K. MacEachern, W. J. 2003 Assessment of 4VsW cod to 2002. Canadian Science Advisory Secretariat Research Document, 2003/027. 41 pp.

Fréchet, A., Gauthier, J., Schwab, P., Pageau, L., Savenkoff, C., Castonguay, M., Chabot, D., et al. 2005 The status of cod in the Northern Gulf of St Lawrence (3Pn, 4RS) in 2004. Canadian Science Advisory Secretariat Research Document, 2005/060. 75 pp.

Gregory, D.N. 2004 Climate: A Database of Temperature and Salinity Observations for the Northwest Atlantic. Canadian Science Advisory SecretariatResearch Document- 2004/075

Hedges, L.V. Olkin, I.1985 Statistical methods for meta-analysis. Academic Press, New York, 369 pp.

Houde, E. D. 2008 Emerging from Hjort’s Shadow. J. Northw. Atl. Fish. Sci., 41, 53–70.

ICES. 2005. Spawning and life history information for North Atlanticcod stocks. ICES Cooperative Research Report, 274. 154 pp.

ICES. 2006. Report of the Working Group on the Assessment of Northern Shelf Demersal Stocks (WGNSDS), 9–18 May 2006. ICES Document CM 2006/ACFM: 30. 870 pp.

ICES. 2007 Report of the Working Group on the Assessment of Demersal Stocks in the North Sea and Skagerrak (WGNSSK), 5–14 September 2006. ICES Document CM 2007/ACFM: 35. 1160 pp.

Lilly, G.R., Wieland, K., Rotschild, B., Sundby, S., Drinkwater, K., Brander, K., Ottersen, G., Carscadden, J., Stenson, G., Chouinard, G., Swain, D., Daan, N., Enberg, K., Hammill, M., Rosing-Asvid, A., Svedäng, H. Vázquez, A. 2008 Decline and recovery of Atlantic cod (Gadus morhua) stocks throughout the North Atlantic. In: Resiliency of gadid stocks to fishing and climate change (eds Kruse GH, Drinkwater K, Ianelli JN, Link JS, Stram DL, Wespestad V, Woodby D) pp. 39-66. AlaskaSeaGrantCollege Program, Faribanks, Alaska

MacKenzie, B.R. & Schiedek, D. 2007 Daily ocean monitoring since the 1860s shows record warming of northern European seas. Global Change Bioogyl, 13(7), 1335-1347.

Mayo, R. K. Col, L. 2005 Gulf of Maine cod. In 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, pp. 2.153–2.184. Ed. by R. K. Mayo, and M. Terceiro. NortheastFisheriesScienceCenter Reference Document, 05–13. 499 pp.

Nissling, G., Westin, L. 1991 Egg mortality and hatching rate of Baltic cod (Gadus morhua) in different salinities. Marine Biology,111, 29–32.

O’Brien, L., Munroe, N. J.Col, L. 2005 Georges Bank Atlantic cod. In Assessment of 19 Northeast groundfish stocks through 2004. 2005 Groundfish Assessment Review Meeting (2005 GARM), Northeast Fisheries Science Center, Woods Hole, Massachusetts, 15–19 August 2005, pp. 2.2–2.29. Ed. by R. K. Mayo, and M. Terceiro, NortheastFisheriesScienceCenter Reference Document, 05–13. 499 pp.

Osenberg, C.W., Sarnelle, O., Cooper S.D. Holt D.H. 1999 Resolving ecological questions through meta-analysis: Goals, metrics, and models. Ecology, 80, 1105-1117.

Ottersen, G., Michalsen, K., Kakken, O. 1998 Ambient temperature and the distribution of Northeast Arctic cod. ICES Journal of Marine Science, 55, 67–85.

Planque, B. Frédou, T. 1999 Temperature and the recruitment of Atlantic cod (Gadus morhua). Canadian Journal of Fisheries and Aquatic Sciences, 56, 2069-2077.

Power, D., Healey, B. P., Murphy, E. F., Brattey, J., Dwyer, K. 2005 An assessment of the cod stock in NAFO Divisions 3NO. NAFO Scientific Council Research Document, 05/67.

Pyper, B.J. & Peterman, R.M. 1998 Comparison of methods to account for autocorrelation in correlation analyses of fish data. Canadian Journal of Fisheries and Aquatic Sciences, 55, 2127-2140, plus the erratum printed in CJFAS,55, 2710.

Rijnsdorp, A., Peck, M.A., Engelhard,G.H., Möllmann,C., & Pinnegar,J.K. 2009 Resolving the effect of climate change on fish populations. ICES Journal of Marine Science66, 1570–1583.

Stuart, A., & Ord, J.K. 1987 Kendall’s advanced theory of statistics, Volume 1 Distribution theory. OxfordUniversityPress, New York, USA, 704 pp.

Vázquez, A. Cerviño, S. 2002 An assessment of the cod stock in NAFO Division 3M. NAFO Scientific Council Research Document, 02/58.

Worm, B., Myers, R.A. 2003 Meta-analysis of cod-shrimp interactions reveals top-down control in oceanic food webs. Ecology, 84, 162–173.

Zar, J.H., 1999 Biostatistical Analysis, 4th edition. Prentice-Hall, Inc., Upper Saddle River,NJ, 931 pp.

1

Supplementary Figures

Figure S1. Fisheries statistical areas in the western (a) and eastern (b) N Atlantic. The locations of the cod stocks are listed in Supplementary Table S1.The maps are reproduced with the permission of FAO.

1

Figure S2. Stock-specific (squares) Risk Ratios (RR; see Table 1 for interpretation) and the overall (diamond) RR with 95% confidence intervals (CI’s) estimated with random effects meta-analysis across the Warm cod stocks for (a) the recruitment survival [log(REC/SSB)]and (b) the Ricker model residuals. The width of the diamond represent the 95% CI’s of the overall RR. The vertical black line corresponds to the mean overall RR and isplotted for comparison with the individual estimates. The size of the squares is proportional to the weight of the individual RR in the meta- analysis. The number of successful events in the “Cool” (C) and“Hot” (G) groupsof seasons(see Table 1 and Supplementary Methods) are also given for each stock. The total number of observations within groups is denoted as G+g and C+c, respectively.

Figure S3. Cumulative meta-analyses of temperature effects on recruitment survival of the Warm cod stocks using (a) Risk Ratio (RR), (b) Hedges’ g (HG) and on the Ricker model residuals using Fisher’s z (FZ) for the (c) Warm and (d) Cold stocks. For this purpose, the meta-analysis is repeated by adding each stock progressively in the analysis. Thus, every result in each plot is the meta-analysis outcome based on the stock denoted on the left column and all the above. The point meta-analytic estimates with confidence intervals (95% for (a)-(c) and 90% for (d)) are plotted with grey symbols and grey horizontal lines. The vertical black lines correspond to the final meta-analytic result estimated based on all the stocks of each group, and are plotted for comparison. See Supplementary Table 1 for stock codes.

Supplementary Tables

Table S1. Summary of cod stocks and temperature data used in the study. The codes, geographic locations (shown in Supplementary Figure S1), time-series length, mean spring temperature (T) and data sources for recruitment and SSB are given for each stock. The %C or %W (in italics) is the index of “coldness”[proportion of T observations in the lower (T4oC)] or “warmness”[the upper range (T6.5oC)], respectively; only those stocks having > 25% of their temperature observations in either interval were included in analyses.The presence of auto-correlation (AC) at lag 1 (p 0.1) in (a) SSB, (b) recruitment, (c) recruitment survival and (d) residuals of the Ricker model is also indicated.

The data for the NE Atlantic stocks were extracted from the ICES Stock Assessment Summary Database (ICES DB 2006; unless otherwise stated.

Stock Code / Areas / AC / Years / spring T / Reference for stock data
mean / %C or %W
cod-2224 / W Baltic (IIId -west) / a, b, c, d / 1970 / - / 2005 / 4.62 / 0.25 / ICES DB 2006
cod-2532 / E Baltic (IIId -east) / a, b, c, d / 1966 / - / 2003 / 4.47 / 0.32 / ICES DB 2006
cod-arct / NE Arctic (I, II) / a, b, c, d / 1953 / - / 2003 / 4.37 / 0.34 / ICES DB 2006
cod-3no / Grand bank (3NO) / a, b, c, d / 1959 / - / 2004 / 3.53 / 0.63 / Power et al.(2005)
cod-3m / Flemish Cap (3M) / a / 1972 / - / 1993 / 2.71 / 0.91 / Vázquez & Cerviño (2002)
cod-2j3kl / N Newfoundland (2J3KL) / a, b, c, d / 1962 / - / 1989 / 1.22 / 0.96 / Bishop et al.(1993)
cod-3pn4rs / N Gulf of St. Lawrence (3Pn4RS) / a, b, c, d / 1974 / - / 2003 / 0.54 / 0.97 / Fréchet et al.(2005)
cod-3ps / S Newfoundland (3Ps) / a, b, c, d / 1977 / - / 2002 / 1.24 / 1 / Brattey et al.(2004)
cod-4tvn / S Gulf of St. Lawrence (4TVn) / a, b, c, d / 1953 / - / 2004 / 0.83 / 1 / Chouinard et al.(2006)
cod-347d / North Sea (IIIa-IV-VIId) / a / 1963 / - / 2005 / 6.45 / 0.45 / ICES (2007)
cod-gb / Georges Bank (5Z) / a / 1978 / - / 2004 / 7.07 / 0.67 / O'Brien et al.(2005)
cod-iceg / Iceland (Va) / a, c / 1956 / - / 2003 / 7.16 / 0.82 / ICES DB 2006
cod-viia / Irish Sea (VIIa) / a / 1968 / - / 2005 / 7.55 / 0.92 / ICES (2006)
cod-7e-k / CelticSea (VIIe-k) / a / 1971 / - / 2005 / 10.63 / 1 / ICES DB 2006
cod-farp / Faroe Plateau (Vb) / a, b, c, d / 1961 / - / 2004 / 7.74 / 1 / ICES DB 2006
cod-via / W Scotland (VIa) / a / 1978 / - / 2004 / 8.79 / 1 / ICES (2006)
cod-gom / Gulf of Maine (5Y) / a / 1982 / - / 2004 / 4.86 / Mayo & Col (2005)
cod-4x / W Scotian Shelf (4X) / a / 1983 / - / 2000 / 4.27 / Clark et al.(2002)
cod-4vsw / E Scotian Shelf (4VsW) / a / 1970 / - / 2001 / 5.05 / Fanning et al.(2003)
cod-coas / Norwegian Coastal (IIa) / a, b, c, d / 1984 / - / 2004 / 5.93 / ICES DB 2006
cod-kat / Kattegat (IIIa -east) / a, b / 1971 / - / 2004 / 5.68 / ICES DB 2006

Table S2. The number of observations (,), the mean (,) and the standard deviation (,)for the “Cool” (or Control) and “Hot” (orExposed) groups of seasons (denoted with C and E subscripts, respectively) used to estimate Hedges’ g (HG) for the Warm stocks. See Supplementary Table 1 for stock codes.