Reproductive investment in Atlantic cod populations off Newfoundland: Contrasting trends between males and females

Three cod populations were considered in this study (Fig. 1): the northern cod stock (NAFO Divisions 2J, 3K, and 3L), for which data from only Division 3L were available; the southern Grand Bank cod stock (NAFO Divisions 3N and 3O); and the St. Pierre Bank cod stock (NAFO Subdivision 3Ps). The pattern in fishing mortality experienced by each of these populations since the late 1950s or early 1960s is available from cohort analyses (Bishop et al. 1993; Brattey et al. 2001; Rideout et al. 2013). Because cohort analyses have not been accepted as being representative of the population dynamics of the northern cod stock since 1993 and of the St. Pierre Bank stock since 2005, trends in total mortality of these stocks have also been deduced from catch rates during bottom-trawl surveys. Mortality values derived from these two approaches (models of population dynamics and catch-rate analysis) do not always coincide and can display a high interannual variability.

Nevertheless, periods with different mortality regimes can be identified. Fishing mortality increased during the 1960s and early 1970s and then decreased considerably in the late 1970s, as limitations on international fishing fleets resulting from the establishment of the Canadian 200 nautical mile exclusive economic zone came into effect. Fishing mortality increased again through the latter half of the 1980s and into the early 1990s, prompting the introduction of moratoria between 1992 and 1994 on directed fishing for each of the three cod stocks. For the northern cod stock, the moratorium remains in effect for the offshore, but relatively small directed fisheries have been permitted in the inshore between 1998 and 2002 and from 2006 to the present. The southern Grand Bank stock remains under a moratorium, but directed fishing on the St. Pierre Bank stock recommenced in 1997. For the northern cod stock, the cessation of directed fishing did not immediately result in a reduction of mortality, as total mortality computed from the research surveys continued to be very high for at least a decade (average total instantaneous mortality rate Z = 1.05 year⁻¹ over the period from 1995 to 2005, for individuals between ages 4 and 6 years), presumably because of heightened natural mortality (Lilly et al. 2013). Total mortality appears to have declined during the most recent decade (DFO 2014). For the southern Grand Bank stock, fishing mortality declined to very low levels after the introduction of the moratorium (average Z = 0.40 year⁻¹ between 1995 and 2000 for ages 4–6 years) but increased to moderate levels for about a decade because of high bycatch (average Z = 0.53 year⁻¹ between 2001 and 2010 for ages 4–6 years; Shelton and Morgan 2005; Rideout et al. 2013). Fishing mortality has been low since 2008 (Rideout et al. 2013). For the St. Pierre Bank stock, fishing mortality declined to a very low level for the three full years of the moratorium (average F = 0.005 year⁻¹ for ages 4–6 years; Brattey et al. 2001) and then increased to a moderate level when fishing recommenced (average F = 0.076 year⁻¹ between 1997 and 2001 for ages 4–6 years; Brattey et al. 2001). Total mortality, as estimated from research vessel surveys, has been moderately high for the past decade (DFO 2015).

All samples were collected during stratified-random bottom-trawl surveys conducted in the spring by Fisheries and Oceans Canada (DFO) from 1978 to 2013. The surveys took place during the early part of the spawning season (with the middle of the spawning season usually occurring in May or June; Hutchings and Myers 1994). For each trawl haul, the date, position, and bottom depth were recorded. The fish caught were subsampled by length classes, sexed, and aged by otolith reading; their length was measured to the nearest centimeter; and their gutted weight, as well as the weights of their livers and gonads, was recorded. The maturity stage of the gonads was defined based on morphological criteria, according to a classification scheme adapted from Templeman et al. (1978).

Only females with ovaries exclusively containing opaque eggs visible to the naked eye, corresponding to stage “Mat A-P” of Templeman et al. (1978), were selected for analyses. This selection was made to prevent the weight of gonads being affected by the hydration of the oocytes, which would lead to an increase in gonad weight that did not reflect changes in the actual energetic reproductive investment, and to reduce the risk of the inclusion of fish that had already begun to spawn. Correspondingly, only males with thick but non-running testes, corresponding to maturity stage “Mat P” of Templeman et al. (1978), were selected for analyses. However, a few isolated peaks were noticed in the distributions of GSI values. To remove these suspicious observations, individuals with a gonadosomatic index (GSI, defined here as gonad weight divided by gutted weight W_gutted) below 2% were removed from the analyses (10.2% of individuals on average). We suspect that some females had already started to spawn at the time of sampling and no longer contained hydrated oocytes; these individuals would, thus, have been erroneously classified as pre-spawning with an abnormally low gonad weight. An analogous filtering was applied to males, because a bimodal distribution of their GSI values, with a peak at very low GSI values (< 1%), suggested that some of them may also have been misclassified (4.5% of individuals on average). Furthermore, to reduce the influence of declining numbers of old cod in the later part of the time series, individuals older than 10 years were excluded from the analysis. The final numbers of individuals used for the analyses are given in Table 1.

All experiments proceeded according to the Canadian Council on Animal Care guidelines (CCAC 2005).

To explain the variability in gonad weight, we explored a number of variables describing the sampled individuals, sampling times and locations, and overall environmental conditions.

To account for individual body condition, the relative liver index was used as an explanatory variable, which is not confounded by the allometry between liver weight and body weight. This index is defined as W _liver/W _liver(predicted), where the predicted liver weight is obtained from the allometric relationship W _liver(predicted) = aW^b _gutted, with parameters a and b estimated separately for each sex and stock.

To account for seasonal effects, the sampling time during the year was used as an explanatory variable, expressed as the day of the year on which the sampling occurred (from 1 to 365 or 366).

To account for spatial effects, sampling position was used as an explanatory variable, described by its latitude (°N), longitude (°W), and bottom depth.

To account for density-dependent effects, stock biomass was used as an explanatory variable. Indices of stock biomass originate from the same surveys as individual biological data. The first year for which stock indices are available is 1985 for Division 3L, 1984 for Divisions 3NO, and 1983 for Subdivision 3Ps. No suitable biomass indices are available from earlier surveys because no method yet exists to make data collected with earlier and present sampling protocols comparable. This limitation restricts the time span over which stock biomass indices can be included as an additional explanatory variable in the statistical analyses below. The most recent year of the series considered is 2013.

To account for temperature effects, ocean temperature was used as an explanatory variable, considering two time series reflecting ocean temperature (Fig. 2). The first time series is the mean annual temperature across all depths (0–175 m) measured at the hydrographic monitoring Station 27 (47°32.8′N, 52°35.2′W), located about 7 km off St. John’s (Newfoundland, Canada; Fig. 1). The deviation of this temperature from its long-term mean (1981–2010) (i.e., the temperature anomaly) was used to represent changes in the average temperature conditions for all stocks. Although Station 27 is located in the inshore branch of the Labrador Current, temperatures measured there are representative of ocean temperatures over the Newfoundland Shelf up to a depth of 175 m (Petrie et al. 1991). The second time series is the area of the vertical cross section of the cold intermediate layer (CIL) along the Bonavista transect (Fig. 1) measured during summer. The CIL is a layer of cold water (< 0 °C) between the seasonally heated surface layer and the warmer bottom layer originating from the continental slope (DFO 2008). In comparison with water at Station 27, the CIL represents a deeper water column, has a broader cross-shelf extent, and is less influenced by near-surface conditions. The larger the CIL area, the colder the oceanographic conditions during the previous winter and spring. As measurements are not available for the years 1978 and 1982 (water temperature anomaly) or 1980 and 1982 (CIL area), the missing values were estimated through linear interpolation. For both of these time series, the extent to which they represent conditions at a biological sampling site likely declines with the distance from that site (Petrie et al. 1991).

Following Roff (1983) and Tomkins and Simmons (2002), we assume that gonad weight is allometrically related to body weight, W _gonad = aW^b _body, where a is a scaling parameter and b is the allometric exponent. This formulation is preferable to using GSI as the response variable because it allows the estimation of the allometric relationship, rather than being confounded by it. As a special case, isometry (b = 1) applies when gonad weight is proportional to body weight, whereas b > 1 corresponds to positive allometry (relative gonad weight increases with body weight). By logarithmic transformation of weight measurements, the above equation can be fitted to data using a generalized linear model. The simplest model, without additional explanatory variables, thus, takes the form:We estimated this model as a generalized linear model with a logarithmic link function and gamma error distribution; the latter was chosen to account for the observed increase of variance in gonad weight with body weight. Body weight was always measured as gutted weight (rather than as total body weight) to remove the confounding effect of the weight of the contents of the stomach and intestine.

We treated stocks and sexes independently and estimated separate models for the entire time span covered by the biological data (1978–2013) and for the shorter time span also covered by biomass data (1983/1984/1985–2013, depending on the stock). For both time spans, we considered the following as potential explanatory variables: gutted weight (log-transformed), age, body condition (relative liver index), sampling day, the geographical coordinates and depth of the sampling station, NAFO division (for the southern Grand Bank stock), annual water temperature anomaly, and CIL area. For all tested variables, we applied a significance threshold of 0.01, to reduce false positives (type-I errors).

The model selection in our analyses was performed as follows. We started from a model including all explanatory variables apart from sampling year, as well as three interactions judged a priori as being potentially important, as follows: the interaction between age and gutted weight, the interaction between relative liver index and sampling day (to account for a potential decrease of the relative liver index at the expense of gonads during a season), and the interaction between temperature and sampling day (to account for the potential influence of temperature-induced changes in spawning time). Regarding this last interaction, both CIL area and water temperature anomaly were tested to represent temperature variability. We then identified the best model using AIC-based stepwise model selection, i.e., stepwise deletion or addition of a single term that results in the lowest AIC during each step. After that, the significance of the linear and quadratic effects of sampling year, added as the last explanatory variables, was evaluated using a likelihood-ratio test. The quadratic effect was included to allow for nonlinear time trends; in particular, for the leveling off and reversal of time trends following reductions in fishing pressure. Because the latter effects are expected to be gradual, a quadratic term is more suitable than a change-point approach to account for them. To be able to test and interpret linear and quadratic effects independently (Schielzeth 2010), the quadratic term for sampling year was centered such that the middle of each time series corresponded to the value zero. When considering the shorter time spans for which biomass data are available (1983/1984/1985–2013, depending on the stock), the best model was selected in the same way as for the full time series, but with the stock biomass index included as a potential additional explanatory variable. This resulted in 12 final models, one for each combination of sex (2 levels), stock (3 levels), and time span (2 levels).

We illustrate the results of our statistical analyses by presenting model-predicted GSI values for a standardized, “typical” individual with a gutted weight of 900 g. Whenever relevant, the following reference values are used: an age of 5 years; a relative liver index of 1; a sampling day of 1 May; a sampling depth of 100 m; the average longitude and latitude of samples in the considered NAFO area(s); a water temperature anomaly of 0 °C; a CIL area of 26 km² (corresponding to the average over the years 1978–2013); and stock biomass indices of 9.17 10⁴ metric tons for Division 3L, 1.32 10⁵ metric tons for Divisions 3NO, and 6.58 10⁴ metric tons for Subdivision 3Ps (corresponding to the averages for each stock during the study period).

The single most important variable explaining variability in gonad weight is, as expected, body weight: larger fish have larger gonads. According to our findings, this relationship is isometric only for females in one of the three studied populations. For females in the other two populations and for males in all three populations, the relationship shows positive allometry, so that GSI increases with body weight. We are not aware of earlier studies documenting such positive allometry in male cod, although it has been demonstrated in other species (Stoltz et al. 2005).

Because the primary drivers of variation in gonad weight, beyond the effect of body weight, are poorly known in the cod populations studied here, we have tested an array of candidate explanatory variables, using standard model selection techniques to find the “best” model for each combination of sex, stock, and time span. When reporting on the resulting 12 models, there is an increased risk of a significant effect appearing by chance alone. We have reduced this risk of spurious effects by testing only variables and interactions that we considered biologically reasonable a priori, and when in doubt, we have checked that our main conclusions do not depend on including potentially spurious effects.

Among individual-level explanatory variables, the best models often include the relative liver index, which shows a positive relationship with gonad weight. This result is biologically expected and in agreement with earlier studies by Yoneda and Wright (2004), who reported a positive correlation between the relative liver index and the potential fecundity of female cod for various stocks off Scotland, and by Skjæraasen et al. (2006), who found a weak positive relationship between potential fecundity and lipid content at the time of spawning in cod (in cod, lipids are mainly stored in the liver; Kjesbu et al. 1991; Lambert and Dutil 1997). The influence of the relative liver index turns out to be relatively weak in our study as well. Similarly, McIntyre and Hutchings (2003) concluded that, after accounting for the effect of body weight, accounting for the effect of the body condition index did not significantly improve the description of the potential fecundity of cod.

Among other significant explanatory variables, those accounting for sampling day (day of the year) and sampling location (longitude and latitude) frequently explain some variability in gonadal investment (with sampling day and longitude accounting for more variability than latitude). As sampling date and sampling location are only weakly correlated, we consider it likely that these statistically estimated effects are real biological effects. The increase in gonad weight throughout the season, accounted for in all our 12 models via the inclusion of the sampling day, illustrates the variability in gonad weight within the considered maturity stages. A significant geographical gradient could reflect a movement of individuals in the most advanced maturity stages toward spawning sites. This is in line with earlier reports of the geographic distribution of cod individuals depending on their maturity stage (Morgan and Trippel 1996; Lawson and Rose 2000; Robichaud and Rose 2003). These studies, together with direct tank observations (Hutchings et al. 1999), suggest the existence of a lek-like mating system in the studied cod populations, with male cod in spawning condition found on the continental shelf together with spawning females, whereas females with premature ovaries and spent females are more abundant on the continental slope.

Interpreting the relationships between gonad weight and water temperature anomaly (showing a positive correlation when present) and between gonad weight and the area of the vertical cross section of the CIL (showing positive or negative correlations when present) requires some care. The CIL area and water temperature anomaly are negatively correlated (Pearson’s correlation coefficient r = −0.70 for the years 1978–2013), and they never occur together in any of the 12 best models. Although the measurements underlying water temperature anomaly and CIL area are partly taken in Division 3L and are, thus, thought to reflect conditions in this division, they are also thought to be representative of the prevailing environmental conditions in Divisions 3NO, but less so those in Subdivision 3Ps. The estimated effect on gonad weight could reflect direct effects of temperature on physiological rates, or effects of oceanographic conditions on cod distribution and food availability. A positive effect of water temperature anomaly and a negative effect of CIL area on reproductive investment seem reasonable for populations living near the cold end of the species’ range (Kjesbu et al. 1998; Yoneda and Wright 2005; Wright et al. 2011). Moreover, Pinhorn (1984) proposed the influence of the water temperature as a potential cause of differences in fecundity observed among the stocks considered here. However, in Subdivision 3Ps, CIL area has a positive effect on the gonadal investment of females, which suggests a different mechanism linking oceanographic conditions and reproductive investment. For example, Hutchings and Myers (1994) found a positive effect of cold water on spawning time and attributed this to cold on-shelf conditions in the area forcing the fish to stay in deeper warmer waters. Water temperature anomaly and CIL area do not have a large influence on our results. Importantly, disregarding these variables does not affect the significance of the trends reported in this study except for females in Subdivision 3Ps, where a significant negative trend in gonad weight appears for the full time span. In this sense, the consideration and potential inclusion of temperature effects in our models is conservative.

Biomass does not seem to be a key factor influencing gonad weight in the studied cod populations. The biomass index is retained in models for both sexes in Divisions 3NO and for females only in Subdivision 3Ps; when present, its influence on gonad weight is consistently negative. A negative effect of population density on fecundity has been shown in many fishes (e.g., Kipling and Frost 1969; Healey 1978; Koslow et al. 1995), including the cod stocks studied here (Stares et al. 2007). As a robustness check, we have confirmed that the choice between considering the stock biomass either in the sampling year or in the previous year (when cod start to divert energy to the maturation of gametes to be released during the next spawning event) does not modify our conclusions about time trends, except for males in Divisions 3NO, where the latter choice suppresses the linear trend.

Of course, a stock biomass index is only a proxy for the per capita availability of resources and does not capture changes originating from the overall availability of resources. For example, the population size and availability of capelin are known to be important for the growth and reproduction of northern cod (Rose and O’Driscoll 2002), but quantifying this interaction over long time periods would be challenging. Nevertheless, because the influence of stock biomass on fecundity—or, in our study, gonad weight—is mediated through the nutritional status of individuals, we may expect this influence to be stronger when the energetic condition of individuals (e.g., through the relative liver index) is not accounted for. In agreement with this expectation, the effect of biomass on gonad weight becomes significant in one case (females in Division 3L) when the relative liver index is excluded from the model. In some other cases, however, the biomass effects on gonad weight remain significant even when the relative liver index is included, which is indicative of a more complex link between stock biomass and gonad weight.

Even though some interannual variability has been reported in the spawning time of cod (Hutchings and Myers 1994; Morgan et al. 2013), this variability seems to be captured by the various individual characteristics (e.g., relative liver index) or sampling characteristics (e.g., longitude) we have considered. As a consequence, the interaction between sampling day and CIL area or water temperature anomaly appears not to bring any additional information, and hence, has never been retained in any of the 12 best models.

Our results reveal positive residual long-term trends in the gonadal investment of male cod in each of the three populations after the effects of individual, sampling-related, and environmental explanatory variables have been accounted for. Although we cannot rule out the possibility of unaccounted variables acting on phenotypic plasticity, the positive trends observed in males at least during the first halves of the time series are in accordance with the hypothesis that the period of heavy fishing has driven the evolution of higher reproductive investment.

A traditional assessment of the reproductive biology of cod suggests that males should be less likely to show evolutionary responses in gonadal investment than females. In a fish species like cod, featuring female mate choice (Hutchings et al. 1999; Nordeide and Folstad 2000; Rowe et al. 2007), male reproductive success should depend less strongly on the number of produced gametes than female reproductive success. In other words, the importance of gonad weight in the determination of fitness is expected to be lower for males than for females. In support of this rationale, Rakitin et al. (2001) did not find a significant relationship between reproductive success and the weight loss of male cod during spawning under experimental conditions. However, GSI values are similar between the two sexes (Fig. 4; Burton et al. 1997), suggesting that an investment in large gonads is also important for male cod. Indeed, the production of higher numbers of spermatozoa can confer an advantage, as male cod mate with several females over the course of a spawning season (Hutchings et al. 1999; Bekkevold et al. 2002). Given that cod sperm remains viable for long periods in water, this may allow male cod to fertilize unfertilized eggs from other mating pairs (Trippel and Morgan 1994).

Our analyses have provided no support for the hypothesis that increased mortality rates have caused the evolution of increased gonadal investment in female cod. However, there are several reasons why an actual evolutionary change might not be detected by our analyses, i.e., why they might be prone to false negatives (type-II errors). First, our methodology is conservative, in that we tested for the significance of the effect of sampling year only after all other significant covariates had been accounted for; specifically, we used a likelihood-ratio test to evaluate whether the inclusion of sampling year could improve a model in which several other significant explanatory variables were already present.

A second explanation for the absence of evidence for an evolutionary increase in female gonadal investment stems from the massive plasticity associated with this metric in cod, as already shown by McIntyre and Hutchings (2003). Naturally, the high variability in gonad weight leads to a reduced statistical power of the analyses. We have tried to minimize this problem by focusing on a single maturity stage for females. Oocytes go through significant changes that influence egg diameter during the late maturity stages (Kjesbu et al. 1990; Kjesbu 1991), and the incorporation of proteins into the oocytes during these stages leads to an increase in gonad weight (Eliassen and Vahl 1982). GSI is not static in males either, but even though male GSI displays a strong increase in the course of spermatogenesis (Dahle et al. 2003), the production of spermatozoa takes place several months prior to spawning (Burton et al. 1997; Rideout and Burton 2000), leading to the relative stability of male GSI values during the maturity stage considered here (Almeida et al. 2008). Indeed, our models for females explain a systematically lower portion of deviance and have higher standard errors for estimates of the effect of sampling year than our models for males. Nevertheless, the point estimates for linear time trends are consistently negative for females (Table 1), suggesting that the absence of positive trends is genuine, rather than representing false negatives (type-II errors).

The prediction that fishing favors increased reproductive investment is qualitatively well established. Only recently, however, have quantitative predictions from eco-genetic models, parameterized for female life histories in cod-like species, shown that fisheries-induced selection not only favors increased GSI, but also that the resultant GSI evolution is slow relative to maturation evolution (Dunlop et al. 2009; Enberg et al. 2009; but see Andersen and Brander 2009, who reached divergent conclusions). It is, therefore, possible that the time period during which two indispensable conditions hold jointly in the three studied populations—that we have adequate quantitative observations and that fisheries-induced selection was strong—has been too short (up to 20 years, and less for analyses including stock biomass) to cause detectable evolutionary change. Moreover, since the exploitation of the studied cod populations had already intensified in the 1950s (Lilly et al. 2008), it is possible that adaptations to high mortality rates occurred before the time period covered in this study. Yet, we argue that adaptation to fishing is unlikely to have halted just a few decades after the industrialization of the Newfoundland cod fishery.

Our results show that changes in gonadal investment in males have leveled off or reversed in all studied stocks toward the end of the study period. This is in agreement with the prediction that a reduction in total mortality should reverse the evolution of reproductive investment, albeit at a slower pace (Dunlop et al. 2009; Enberg et al. 2009). Such a reversal in gonadal investment should become apparent only with a delay of several years, because it becomes expressed first when the progeny of selected parents reach maturity. The reported negative trends in the gonadal investment of females could also represent such an evolutionary reversal. However, these potential reversals should be interpreted with caution, as their degree is difficult to predict for the stocks studied here, for at least two reasons. First, even after the directed cod fisheries were closed, natural mortality may have remained high (Shelton et al. 2006; Lilly et al. 2013), and selection for increased reproductive investment may, thus, have persisted, albeit at reduced intensity. Second, some cod fishing has taken place during the moratoria in all stocks, and the moratorium for St. Pierre Bank stock was lifted in 1997 (DFO 2015).

In summary, the results of the analyses performed here show evidence for an evolutionary change in gonadal investment for cod populations off Newfoundland, but only for male cod. This suggests that gonad weight has greater importance for the reproductive success of male cod than has previously been acknowledged. The lack of a similar response in female cod could be explained by the relatively short time span with strong fisheries selection, in conjunction with the higher phenotypic variability in gonadal investment in females compared with males; but this does not explain why the available evidence points to negative trends in females, if any. Our findings highlight that investigating the ecological and evolutionary determinants of reproductive investment in wild fish is not straightforward, and underscore how little we know about the reproduction of wild cod in particular, thus calling for more studies unraveling the reproductive biology of this important species.