This investigation commenced by constructing principal maturation schedule equations as a function of fishing mortality (F), key biophysical factors and a term attributed to fisheries-induced adaptive change (FIAC). Following the onset of industrial trawl fishery on the model stock, Northeast Arctic cod (NEAC; Gadus morhua) (1934–2020), F on immature age group 5–8 years (F5-8) increased and mean age at 50% maturity (A50) decreased from ≈10 years in the late 1940s to ≈7 years today. Large annual fluctuations in total stock biomass (TSB), sea temperature (KolaT) and F5-8 were used to better understand A50 responses. In the model, the annual accumulation of F5-8 drives FIAC. The model includes the option that NEAC may sustain F5-8 up to a certain level (F_bal) before FIAC becomes statistically evident, with F_bal falling between 0.00 and 0.40 for A50. This dynamic range in F_bal indicates a sophisticated, underlying adaptive response. Independent of F_bal, our analysis clarifies that FIAC is necessary to explain the observed changes in A50.
La présente étude commence par l’établissement d’équations principales décrivant la chronologie de maturation en tant que fonction de la mortalité par pêche (F), de facteurs biophysiques clés et d’un terme attribué au changement adaptatif induit par la pêche (CAIP). Après le début de la pêche industrielle au chalut du stock utilisé comme modèle, la morue du nord-est de l’Atlantique (MNEA; Gadus morhua) (1934–2020), la F de groupes d’âge immatures de 5–8 ans (F5-8) a augmenté et l’âge moyen à 50 % de maturité (A50) est passé de ≈10 ans à la fin des années 1940 à ≈7 ans aujourd’hui. De grandes fluctuations annuelles de la biomasse totale du stock (BTS), de la température de la mer (KolaT) et de la F5-8 ont été utilisées pour mieux comprendre le comportement de l’A50. Dans le modèle, la F5-8 annuelle cumulative module le CAIP. Le modèle comprend l’option pour le stock de MNEA de supporter une F5-8 jusqu’à une certaine valeur (F_bal) avant que le CAIP ne devienne statistiquement évident, la F_bal chutant à des valeurs entre 0,00 et 0,40 pour l’A50. Cette fourchette dynamique de la F_bal indique une réaction adaptative sous-jacente complexe. Indépendamment de la F_bal, notre analyse fait clairement ressortir la nécessité du CAIP pour expliquer les variations de l’A50 observées. [Traduit par la Rédaction]
Generally, fisheries exploitation changes the age and body size structure of harvested populations (stocks; Ottersen et al. 2006; Shelton et al. 2015). The more intensively a stock is fished, the more likely it is made up of mostly younger and smaller fish (Law and Grey 1989; Trippel et al. 1997; Scott et al. 2006). Furthermore, there are associated trends towards an increased prevalence of early sexually mature individuals, as seen in, e.g., Canadian cod (Gadus morhua) stocks (Trippel 1995). Whereas these changes could reflect genetic responses to fishing mortality (F) (Law and Grey 1989; Heino 1998; Heino and Godø 2002), they may also possibly occur without any underlying changes in the genome as such, e.g., by epigenetic, transgenerational modifications (Heckwolf et al. 2020). In addition to these more fundamental aspects, phenotypic plasticity may relate to increased resource (prey) availability and faster growth of survivors when the stock biomass decreases (Olsen et al. 2005). Maturity ogives describe the proportion (or probability) of sexually mature individuals in a population as a function of age or size. The age at 50% proportion mature, traditionally called A50, is a standard indicator of age at sexual maturity (Trippel 1995). In other words, A50 depends on maturation disposition at the given age and size, implicitly also on survival and growth until that age and size. The factors introducing temporal changes in A50, detailed by Marty et al. (2014), can roughly be grouped as demographic, i.e., population age and size structure (across co-existing generations), phenotypic plasticity (within a given generation), and evolution (over succeeding generations). The demographic situation affects A50 because larger individuals within an age class are more likely to be mature and more vulnerable to fishing with size-selective fishing. This additive effect reduces the proportion mature within the age-class and increases A50 (Marty et al. 2014). In their study, Marty et al. (2014) established that A50 and L50 (length at 50% maturity) decreased from the 1970s to the 2000s for three North Sea gadid stocks, cod, haddock (Melanogrammus aeglefinus) and whiting (Merlangius merlangus), but not for Norway pout (Trisopterus esmarkii). These maturation trends were addressed by estimating probabilistic maturation reaction norms (Heino and Godø 2002). However, Marty et al. (2014) stated that the body of evidence would be more complete with empirical analysis of selection differential due to fishing. Our work aims to provide further insights into that research direction.
The Northeast Arctic cod (NEAC) stock inhabits the Barents Sea, an open shelf sea of around 1.4 million km2 in the Northeast Atlantic — north of Norway and northwestern Russia between around 70°N and 80°N (Bogstad et al. 2016). Mature specimens undertake extensive migrations, against the prevailing currents, during winter-spring to spawning grounds along the Norwegian coast (Yaragina et al. 2011). Spawning grounds show a northward displacement in warmer periods and a southward displacement in colder periods (Sundby and Nakken 2008). As opposed to many other Atlantic cod stocks, the NEAC stock is presently in a healthy state, i.e., classified by the International Council for the Exploration of the Sea (ICES) to have “full reproductive capacity” (ICES 2020). The total stock size (TSB; age 3+ years) was at a minimum in the 1980s, but in the last decade TSB has recovered to its post-WWII level of about 4 million tonnes before declining to the present (2020) level of between 2.5 and 3 million tonnes (ICES 2020). Increased stock productivity relates to an expanding habitat area in the Barents Sea at higher sea temperatures (Kjesbu et al. 2014).
Maturation-schedule phenotypes (hereinafter “phenotype”) of NEAC have changed in the post-war period. In the 1930–1950s, a small proportion of each year class matured at 6 and 7 years, the majority between 8 and 12 years, and the remainder at 13 and 14 years (Rollefsen 1954). The body growth rate increased from the 1930s to the 1970s, and, in the same period, the A50 shifted nearly 2 years downwards (Borisov 1978; Jørgensen 1990), while the interquartile range of the maturity curve (A75–A25) for the year classes 1923–1976 decreased (Jørgensen 1990). There is an apparent negative correlation between stock density and mean length at age, whereas there is a positive correlation between growth rate and feeding indices (Ponomarenko et al. 1980). A higher temperature positively affects the mean length of 2- to 5-year-olds but likely not the age of first-time spawners (Michalsen et al. 1998). Age at first-time spawners is sensitive to the level of exploitation, density dependency, food availability and year-class strength (Godø 2000). However, an analysis of the year classes 1923–1989 suggested that both fisheries-induced evolution (FIE) and the compensatory phenotypic response (caused by faster individual growth) resulted in earlier sexual maturation (Heino et al. 2002). A decrease in F will increase the proportion of slowly growing and late-maturing individuals that transmit hereditary traits to the offspring (Kovalev and Yaragina 2009). Regardless, there are apparent gaps in our understanding of phenotypic and genotypic adaptions to climatic fluctuations, especially when combined with resilience to fishing pressure (Yaragina et al. 2011). The NEAC stock is one of the most data-rich fish stocks in the world, with long time series available both on fisheries and biological data. Therefore, the NEAC stock is highly suitable for in-depth, longitudinal analyses, especially since it has been fished down to a very low level and then recovered to a healthy state, giving large variability in statutory population metrics within the mentioned data series (Kjesbu et al. 2014).
Here, we consider the phenotype A50 of NEAC as being a function of TSB, sea temperature, represented by the Kola hydrographic transect (KolaT; see below), and fishing mortality on immature age groups, age 5–8 years (F5-8; Fig. 1). As time series on food availability per se were comparably limited (Holt et al. 2019), this explanatory factor is presently left out. However, in line with the above elaborations, TSB is used as a proxy for intraspecific competition for food, whereas KolaT is set to relate to body growth rate and suitable feeding area extent. We adopted the key principle that high selective fishing pressure should induce changes in life-history traits: the higher the harvest rate and the longer the harvest period, the larger the expected shift in phenotypes, like A50 (Enberg et al. 2009; Dunlop et al. 2009; Marty et al. 2015; Eikeset et al. 2016; ICES 2018). Therefore, cumulative F on immature age groups as an explanatory variable is considered a reasonable approach to model associated changes of phenotypes between generations. However, as current, extensive research on the genome of NEAC (1907 vs. 2011–2014) as well on northern cod (1940 vs. 2013) by Pinsky et al. (2021) (and accompanying Commentary by Hutchings and Kuparinen 2021) provides no unequivocal empirical evidence of FIE. We, therefore, opted for fisheries-induced adaptive change (FIAC), described by Rijnsdorp et al. (2005), as the overall term in these specific respects. Multiple regression models may disentangle FIAC from phenotypic change (Heino et al. 2015). Furthermore, our scientific approach was encouraged by large temporal dynamics in explanatory factors (Fig. 2A), which should in principle allow disentangling variations in TSB and KolaT from any FIAC driven by F. Essentially, the present article focuses on the development of an intrinsic trait of a phenotype that principally defines the sexual maturity ogive: A50. We hypothesised that NEAC, and potentially other long-lived harvested fish populations, may sustain fishing mortality up to a certain level (F_bal) before FIAC is apparent statistically. To our knowledge, this proposition is new.
Materials and methods
All data originated from statutory monitoring programs, and all following analyses were done with the Add-ins packages “Analysis” and “Solver” in Microsoft EXCEL.
The very long population time series available for NEAC from complementary sources were consulted. Hylen (2002) did a virtual population analysis (VPA) assessment of the stock from 1913 to 1999. His analysis did not include cannibalism, and the natural mortality rate (M) was set equal to 0.20 for all age groups and years. However, annual data on cannibalism are available from 1946 (ICES 2015; Yaragina et al. 2018) and are included in the assessment done by ICES (2020). In the lack of historical data on cannibalism, M in 1930–1945 was set equal to the 1946–1983 average, i.e., 0.410 for age 3, 0.274 for age 4 and 0.219 for age 5 years. Numbers at age (N) since 1946 were taken from the most recent analysis (ICES 2020). For the period 1930–1945, terminal N values by year for the oldest true age group (14 years) were extracted from Hylen (2002), F for the plus group (15+ years) were set equal to F for the age group 14 years, and N for age group 3–13 years were calculated using cohort analysis (Pope 1972). All extracted, tabled body weight, maturity, and catch-at-age information was kept unchanged from Hylen (2002). Compared to the assessment in Hylen (2002); an introductory study clarified that the inclusion of cannibalism on age group 3–5 years had minimal effects on total stock biomass (TSB) and F on age group 5+ years. However, it changed the estimated recruitment at age 3 years (N3). As N3 played the role as background information (see below), this effect was not further pursued. TSB estimates refer to the beginning of the year (Fig. 2A).
We organised the maturity data for the year classes 1923–1976 in Jørgensen (1990, see his table 3) by year and estimated annual values of A50. This analysis was limited to 1934–1983, i.e., years when the available data included one or more age groups below and above 50% maturity. The available A50 series are then from Jørgensen (1990) for 1934–1983 and 1946–2020 from ICES (2020). The procedure used by Jørgensen (1990) differs from the ICES standardised procedure (ICES 2001). For the overlapping years, 1946–1983, these two A50 series are close, especially when the values are in the range of 9–11 years (Fig. 2B). Therefore, the A50 series from Jørgensen for 1934–1945 was appended (unadjusted) to the ICES series (1946–2020). This extended series indicates that the major decline in A50 took place when fishing mortalities on immature fish (F5-8) were high (Figs. 2A, 2C). Note that F on immature fish markedly dropped during WWII (Figs. 2A, 2D) as there was no trawl fishery in the Barents Sea. In 1945–1946, the trawlers quickly resumed their pre-war fishing in the Barents Sea (Fig. 2D).
The physical environment of NEAC in the Barents Sea was, as commonly done (Ingvaldsen et al. 2003), reflected by the Kola Transect temperature (KolaT; Boitsov et al. 2012). KolaT is the annual average from seven hydrographic stations along the 33°30′E meridian from 70°30′N to 72°30′N and vertically between 0 and 200 m depth (Bochkov 1982; Tereshchenko 1996; recent data provided by VNIRO Polar Branch, Murmansk). Monthly data were available for 1921–2019, except for the last four months in 2018. As body growth of NEAC mainly occurs between May and November (Jørgensen 1990); annual sea temperature averages over these seven months were presently estimated, starting from 1930 (Fig. 2A). Due to the missing 2018 observations mentioned, KolaT in that particular year is the average from May to August.
We used a multivariable linear model with a minimum number of coefficients. The model incorporates the expectations that sexual maturity at the beginning of the year depends on the intraspecific competition for food (with TSB as a proxy) and environmental temperature (KolaT) in the previous year. The third explanatory variable considered, FIAC, was assumed to have long-term effects across generations. More specifically, FIAC changes slowly due to the accumulated sum of the annual difference between average fishing mortalities on immatures and a noticeable “threshold”, represented by F_bal, in the preceding years. This approach presumes that F_bal remains constant over the years and is independent of TSB and KolaT. However, data on fishing mortalities on immature and mature fish split by age are unavailable; most immature and mature fish catches are jointly taken by trawl in the Barents Sea. Therefore, the trend in F5-8 is expected to represent the harvest rate on both immature and mature fish across age groups. In the pilot study, we tried out different alternatives to F5-8. In one, the lower age was fixed at 5 years, and the upper age varied with A50; different fixed ranges of age groups were tested in other alternatives, but they all made little difference to the results. The apparent explanation should be that fishing mortalities on the different age groups within a year are highly correlated. We decided to use F5-8, the average for age group 5 to 8 years.
The annual impact of TSB and KolaT on A50 was set relative to the base year 1933. Consequently, the mathematical formulation (eq. 1) reads as follows:
where A50(y) is the recorded value of A50 in the beginning of year y; TSB is the given total stock biomass of age group 3+ years (in 1000 t); KolaT is average sea temperature measured in the Kola Section in May–November (in °C) (May–August in 2018, see above); F5-8 is average annual fishing mortality on age group 5–8 years; F_bal is the average fishing mortality on age group 5–8 years the stock can sustain without a long-term, noticeable change in A50 related to FIAC; α is an estimate of A50(1933), β is the rate of change in A50 related to TSB; θ is the rate of change in A50 related to KolaT; γ is the rate of change in A50 related to accumulated fishing pressure in the preceding years (FIAC); and the error term shows a normal distribution: εy ∼ N(0,σ2). It is assumed that there is no correlation between consecutive residuals. The accumulation of F5-8 starts with i = 1933, one year before the first year in the A50 series used in the regression.
The following equation modelled sexual maturity (Mat) as a function of year and age:
where A50(y) is age at 50% sexual maturity in year y, and S50(y) is the slope of the sexual maturity curve in year y at Age = A50(y). Minimisation of gave estimates of A50 and S50 for each year. Note that for the years 1934–1937 the summations were truncated before age 15 years as maturity data for the year classes 1922 and older were unavailable.
The sensitivity analyses also included the possibility that the sexually mature stock consists of two parts, an early-maturing component with a fixed value A50_E and a late-maturing component with a fixed value A50_L. The following combined model for age at maturity (A50) would then be
In these two latter equations, δ(y) is the proportion of the total spawning stock belonging to the early-maturing component each year, and γδ is the FIAC factor for δ(y). When eqs. 3 and 4 are merged, they are reduced to eq. 1, where β and are unchanged, while α = δ(1933)·(A50_E – A50_L) + A50_L and γ = γδ·(A50_E – A50_L). This rearrangement means that eq. 1 gives the necessary parameter values to model the changes in the spawning stock (eq. 4) for different values of A50_E and A50_L, i.e.:
The regression of eq. 1 to the A50 data series is independent of the choice of A50_E and A50_L. Therefore, eqs. 3 and 4 cannot be used to estimate the most likely values of A50_E and A50_L. The first choice of A50_E and A50_L (6.2, 10.8) represents the lowest and highest values of A50 on the record (1934–2020) for the total spawning stock. To test the sensitivity of the results from eq. 4 to A50_E and A50_L, we chose somewhat arbitrarily two additional scenarios. The span between A50 for the early and late maturing components was set to increase (5.7, 12.0). In the third alternative, the span decreased (6.7, 10.0). This split of the spawning stock into two components with fixed A50 values (eq. 3) simplified the probable continuous spectrum of A50. In this two-component model, seven parameters in eqs. 3 and 4 had to be parameterized in advance or estimated in the regression (see below).
Principally, F5-8 < F_bal reverses the trend in A50(y) in the overarching model, eq. 1. The following rationale formed a basis for calculating possible recovery times of A50 to its 1933 level when TSB and KolaT are set equal to their 1933 values. In eq. 1, the FIAC part of the A50 change depends only on the accumulated fishing mortalities in the previous years. Let 2021 be the first year in the recovery phase, the average fishing mortality on the 5- to 8-year-olds be constant (F_const) in the following years, and Yr is the year A50 is restored to its pre-WWII level when F_const < F_bal. Then:
A50(Yr) = A50(1933) gives
So, eq. 8 provides estimates of the recovery time, depending on F_bal and F_const. This model assumes a linear relationship. Hence, it implicitly precludes any magnifying or dampening effect of fishing above F_bal consistently over many years, such as genetic variability of traits. It also assumes that F_bal is constant and independent of environmental change.
Choice of model and F_bal
F_bal is a fixed input parameter in eq. 1, and the trajectory of depends on the value of F_bal (Fig. 2C). The changing trends of A50 allow estimating the likely range of F_bal. We ran a regression analysis (eq. 1) for different values of F_bal in the range of 0–0.50 (Fig. 3A). The R2 curve for this extended A50 series was almost flat when F_bal is in the range 0.0–0.40 (Fig. 3A), followed by a rapid drop. The regressions’ residuals are similar (Figs. 3B–3D), particularly for F_bal = 0.00 and F_bal = 0.20, and our choice of a linear model (see above) seems reasonable. Given the wide range of possible values of F_bal, we chose to do the subsequent analysis using the following three values of F_bal: 0.00, 0.20 and 0.40.
Multilinear regression of time series implicitly assumes that consecutive residuals are not correlated. However, autocorrelation in a time series reduces the effective sample size (Bretherton 2015). Consequently, we explored how replacement of actual sampling size (N) with effective sampling size (N*) affected the calculation of P for the studied impact factors, TSB, KolaT and FIAC at F_bal = 0.00, 0.20 and 0.40 (see above).
Impacts of TSB, KolaT and FIAC
Studying firstly the outputs from separate regressions, FIAC showed the highest explanatory power (Table 1). The further inclusion of TSB and KolaT, when FIAC was already included, contributed little to R2 (Table 1; F_bal = 0.20). The same type of regression analyses with F_bal = 0.00 and 0.40 supported this finding.
Alternative ways of mathematically expressing TSB(y – 1) and KolaT(y – 1) in eq. 1 were also tested. The increase in the goodness of fit (R2) was minimal with 2–3 years running averages, the apparent explanation being that the autocorrelation of TSB in 1934–2020 was as high as 0.93 (see further elaborations below). Therefore, running averages for TSB over several years had minimal effect on R2. The autocorrelation for KolaT was, however, lower (0.47), but, as mentioned below, KolaT generally had little impact on A50, and consequently, several years running average did not improve R2.
For all the three values of F_bal considered (Fig. 2C), outputs from eq. 1 mimicked the observed declining trend in A50 (Figs. 4A, 4C, 4E). Furthermore, there was no significant effect of KolaT on A50, but so for TSB (Table 2), where the more extensive, absolute influence appeared in the middle of the study period (Figs. 4B, 4D, 4F). The long-term influence of FIAC was evidently in place (Figs. 4B, 4D, 4F; Table 2). Analyses of the A50 series from the shorter time intervals (1934–1983 from Jørgensen (1990) and 1946–2020 from ICES (2020)) supported patterns seen in the presently adopted, combined 1934–2020 time series results, but provided no further, relevant insights.
Autocorrelations (lag 1) between A50 residuals were easily detectable (Figs. 3B–3D). As expected, the act of replacing the actual with the lower effective sampling size increased P values (as demonstrated for KolaT), but the high significance level for both TSB and FIAC remained as before for all three values of F_bal (Table 2 vs. Table 3). Thus, the observed autocorrelations did not invalidate our analysis.
In our model, spawning stock demography appeared sensitive to the selected values of A50_E and A50_L (Fig. 5A). However, for the reasons described above, the regressions did not indicate a more likely set of A50_E and A50_L. In this simulation (Fig. 5A), we used the values of α, β and γ for F_bal = 0.20 (Table 2). The respective alternatives of A50_E and A50_L provided estimates of δ(1933) and γδ (eqs. 5 and 6) to be used in eq. 4. Values of α appeared rather insensitive to F_bal, particularly for F_bal = 0.20 compared to F_bal = 0.00 (Table 2), which added an explanation to similar demographic changes for these two values of F_bal (Fig. 5B). However, for F_bal = 0.40, the trajectory of the SSB composition exhibited more inter-annual variability (Fig. 5B). To add, these model investigations should be considered in view of the fact that the mean age of the spawners was reduced by 3–4 years from the post-war level to the 1980s, while in recent years, 2016–2020, the relative age composition has approached the post-war level (Fig. 5C).
Projected recovery time
The higher the value of F_bal, the shorter became the recovery time, ranging from a few decades to more than a thousand years (Fig. 5D). This simulation, using eq. 8, started with a virtual moratorium on age 5–8 years in 2021, i.e., F_const = 0 in eq. 8. For the year 2020, F5-8 was set to 0.261, i.e., equal to F5-8 for 2019 (ICES 2020), implying ∑(F5-8) = 39.55 for 1933–2020.
Background information on recruitment
The number of recruits to the NEAC stock at age 3 years (N3) showed large temporal fluctuations since the 1930s but with a reasonably stabilized level of recruitment since the 1990s (Fig. 6A). Relative recruitment, i.e., N3 divided by parent SSB, showed a declining trend since the 1950s (Fig. 6B).
First, we are not aware of any ecogenetic models that explicitly provide values of fishing mortality a stock may sustain (F_bal) without changing life-history traits noticeably. In these models, FIE is made possible by stochastic processes, which may include variations both in the offspring genetic trait value for a given pair of parents, and in the expression of the phenotype (Enberg et al. 2009; Dunlop et al. 2009; Marty et al. 2015; Eikeset et al. 2016). These stochastic processes also make the recovery of life-history traits possible. Enberg et al. (2009) modelled a stock resembling Atlantic cod in the northern part of its range, defining 100 years of harvest followed by a moratorium. In their model, genetic traits took thousands of years to evolve back to preharvest levels, indicating that natural selection driving recovery of these traits is weaker than fisheries-induced selection. In a similar study, Dunlop et al. (2009) concluded that the recovery of life-history traits from FIE during a moratorium was slow to nonexistent. Given the long recovery times in these ecogenetic studies, they have an “effective” F_bal close to 0 (Fig. 5D).
One would expect that F_bal > 0, otherwise all in-built, natural phenotypic changes would have taken place well before the intense trawl fishery started on immature fish. Logically, there must also be an upper limit for the possible range of F_bal. In the last decade, TSB has recovered to its post-war level, whereas A50 has not done so. In the six decades, from 1950 to 2010, when the decline in A50 and the impact of FIAC was at its highest (Figs. 4A–4F), F5-8 was high (Fig. 2A), on average 0.58 (ICES 2020). One would, therefore, expect F_bal to be well below 0.58. In simulations made by Borisov (1978), the spawning stock of NEAC was divided into three subgroups, the early mature fish (6–7 years), medium (8–10 years) and late-maturing (11+ years). With an annual survival rate of 65% on immature age groups, he found that the ratio of reproduction output from the three subgroups remained nearly constant over time. A survival rate of 65% corresponds to F = 0.23 (when M = 0.20). With the present fishing pattern, the average fishing mortality on the age group 5–10 years that maximises sustainable yield (MSY) is 0.40 (ICES 2020). This FMSY was 61% of the average F5-10 in 1950–2010, whereas the average F5-8 was 0.58 in the same period. A prorate calculation (0.61 × 0.58) gives FMSY on age group 5–8 years to be 0.35. This estimate of FMSY is close to the upper limit of the likely range of F_bal (0.40), while F_bal derived from the ecogenetic studies referred to above is near the lower limit (0.00).
Recent development in A50
As mentioned, in contrast to TSB, A50 has not recovered to its level in 1930–1950. From the early 1980s, when the NEAC stock was at a historical minimum, to the present, the maturity of NEAC has shown only minor fluctuations, with an A50 of around 7 years. The high A50 values around 1987 were related to adverse feeding conditions and low body condition factor (Mehl and Sunnanå 1991). However, the fishing mortality on age group 5–8 years was high (0.52–0.94) from 1980 to 2006 (except for 1990–1992), and low from 2007 to 2019 (0.18–0.39), on the average 0.25. Since the 1980s, the impact of TSB and FIAC on A50 have generally worked in opposite directions (Figs. 4B, 4D, 4F). The net result is that the steep downward trend of A50 halted during the last four decades (Figs. 4A, 4C, 4E). The FIAC is the probable reason why A50, in contrast to TSB, has not recovered to its level in the 1940s.
Composition of the spawning stock
The intense post-WWII trawl fishery altered the spawning stock’s demography (Ottersen 2008). The size-specific selectivity has also changed over time. Before the 1920s, the exploitation was mainly on older fish, using passive gears at the spawning sites. With the development of trawl fisheries, the pattern of exploitation expanded to include also immature fish (Godø 2003). From the 1970s onwards, the exploitation moved towards larger and older fish again. This pattern was related to increased minimum mesh size and mandatory sorting grids in trawls and limits on bycatch of juveniles (Yaragina et al. 2011). The percentage of fish of age 10+ years in SSB decreased from 70% (1946–1955) to 11% in 1980–1989 and then increased to 50% in 2016–2020 (ICES 2020; Fig. 5C). The latter pattern is primarily the result of reduced fishing mortality on younger age groups since 2007 (Fig. 2A). The regressions did not signal the most likely value of (A50_E, A50_L; see above). However, our analysis supports the conclusion that the size-selective fishery since 1946 has increased the share of fish in the (spawning) stock that become sexually mature at an early age. This depiction explains why A50, contrary to TSB, has not recovered to its level 70 years ago. This inference does not depend on the value of F_bal.
The recovery time of A50 depends very much on model settings, represented by F_bal. The estimated recovery time is about 300 years if F_bal = 0.10, decreases to about 100 years if F_bal = 0.20, and is close to 10 years if F_bal = 0.40. As discussed above, ecogenetic models estimate much longer recovery times for age at maturity, i.e., in the order of several thousand years (Enberg et al. 2009). Such ecogenetic models also give different recovery times for different life-history traits. In a laboratory study of Atlantic silverside (Menidia menidia), populations were exposed to size-selective fishing (Conover et al. 2009). After five generations, this experimental design generated a twofold difference among populations in mean body weight. Then followed five generations without size-selective harvest. Conover et al. (2009) concluded that if enough genetic variations remained, populations have an intrinsic capacity to recover from disadvantageous evolutionary changes caused by fishing. More specifically, they predicted evolutionary recovery in about 12 generations. Depending on the value of A50 (6–11 years), the time for 12 generations for NEAC is 100 ± 30 years. A recovery time of this length corresponds to F_bal ≈ 0.20. However, in contrast to cod, Atlantic silverside is short-lived with maximum longevity of 2 years. Experimentally-based extrapolation to a long-lived fish like NEAC may therefore not be valid (see below). The extrapolation of eq. 1 hundreds of years into the future is only an indication that the recovery of A50 may be a slow process. Another factor that may influence the recovery is the sex difference in the onset of maturity, as NEAC males, on average, sexually mature one year earlier and at a smaller length than NEAC females (Ajiad et al. 1999). To correct for this factor, the sex difference in fishing mortality needs to be considered.
Possible causes for FIAC
Instead of the more general term FIAC, FIE refers to evolutionary changes in life-history traits in exploited stocks (Heino et al. 2015). That FIE is a possibility is indicated by aquaculture studies of NEAC, which suggest high variance in body growth under identical environmental conditions (Bangera et al. 2015; Gjerde et al. 2004). Therkildsen et al. (2013) did an exploratory screening of allele frequency shifts at certain loci in four overfished populations of Atlantic cod off Canada, stating that the findings were consistent with temperature- and fisheries-induced selection pressures acting over the 80 years study period. Furthermore, a theoretical study concluded that the expected rates of FIE are ≈0.1%–0.6% per year (Andersen and Brander 2009). This range compares well with the average impact of FIAC on A50 of about 0.60% (Table 2). This estimate is for the period 1950 to 2010 when F5-8 was at its highest and nearly identical for all three values of F_bal used in the sensitivity analysis (Fig. 2A; Table 2). Eikeset et al. (2016) used a multitrait ecoevolutionary model to explore age and length at maturity of NEAC since the 1930s. The noticed decreasing trends were explained by a possible combination of several processes, truncation of the population's age structure, phenotypic plasticity in maturity arising through density-dependent growth, and FIE favouring faster-growing or earlier-maturing fish. Their model also predicted that the role of evolution diminishes when density-dependent growth is present.
Nevertheless, despite the evidence that fishing induces changes in phenotypes, unequivocal empirical evidence of underlying genetic change resulting from fishing is generally lacking (Pandolfi 2009; Heino et al. 2013, 2015). However, in the most recent genome-wide analysis by Pinsky et al. (2021) on NEAC and cod off Newfoundland (see above), phenotypic changes in these populations are apparently not constrained by any irreversible loss in genetic variation, implying that former traits could be re-established with demographic recovery. Therefore, we firmly stuck to FIAC rather than FIE throughout this article.
The cost of early maturation
Earlier maturity (decreased A50) and the accompanying increased survival costs of reproduction tend to result in higher natural mortality and a shorter life-span (Jørgensen and Fiksen 2010; Swain 2011). Beverton et al. (1994) analysed NEAC otoliths for the year classes 1920–1935, collected on the main spawning grounds in Lofoten. They stated that cohorts maturing at 8 years and upwards had nearly the same body growth and adult life-spans as far as it could be judged. The highest natural mortality rate (M) for these older maturing cohorts consistent with the data are 0.15; on that basis, M for cohorts maturing at 6 and 7 years would have been about 0.25 and 0.17, respectively (Beverton et al. 1994). There has been a declining proportion of repeat spawners in the NEAC stock from 1913 to 2004 (Ottersen 2008). The facts accumulated from several experimental and field studies suggest that repeat (and larger) spawners will have larger eggs with less mortality than recruit (first-time) spawners as well as higher fecundity (higher gonad weight/body weight ratio; Kjesbu et al. 1992, 1996; Solemdal et al. 1995; Marshall et al. 1998). Further, older and large NEAC spawn over a wider area than the younger and smaller ones (Stige et al. 2017). It follows that the population's reproductive potential could be enhanced by a heterogeneous stock composition through size- or age-dependent differences in timing, duration, or spawning location, ensuring that a sufficient number of eggs and larvae encounter favourable environmental conditions. Spawning over a large area or a long season should reduce the risk of recruitment failure. Also, there is statistical substantiation that reduced age and or length of spawners will increase climate change sensitivity (Ottersen et al. 2006). However, Ottersen (2008) could not find any clear connection between recruitment and age structure in the population. Still, we indirectly find that earlier maturity, not examined by Ottersen (2008), causes the spawning stock’s total fecundity to decrease and potentially negatively affects recruitment to the fishery, as indicated by a reduction in the average N3/SSB ratio during the last seven decades (Fig. 6B). It is also worth noting that the five largest year classes on record since 1946 were hatched before 1971 (Fig. 6A), that is, in a period when the proportion of late-maturing cod in the spawning stock exceeded 50% (Figs. 5A–5B).
As stated in the Introduction, the present status of the NEAC stock is considered to have “full reproductive capacity” (ICES 2020). In 1946–1954, SSB was, on average, 17% of TSB. In the last ten years, the SSB has been the highest on record, on average 55% of TSB. This is partly due to year classes 2004 and 2005, which at age 3 years were twice the size of the average in 1946–2019 (Fig. 6A). Two other factors to count in are low fishing mortalities facilitated by implementing harvest control rules (Kjesbu et al. 2014) and a reduced age at 50% maturity. Furthermore, the increased influx of warm Atlantic water may increase juveniles’ survival rate (Sætersdal and Loeng 1987; Ottersen and Stenseth 2001) and have expanded the suitable feeding area in the Barents Sea (Kjesbu et al. 2014). Since the turn of the century, the NEAC stock has generally benefited from favourable climate conditions (Gullestad et al. 2020).
Possible consequences for management
According to our hypothesis, the present fishing mortality on age group 5–8 years (0.25 in 2007–2019) is at a level that has reduced the annual change in A50 and possibly halted it or turned into a very early phase of reversing this long-term trend (Pinsky et al. 2021). However, further reducing F towards zero to recover the maturation schedule to historical levels could have large consequences, both biologically and socioeconomically. Biologically, reduced F would certainly increase stock size, but likely cause decreased body growth rates and increased cannibalism during the rebuilding period. One possibility could be to shift the exploitation pattern towards larger and older fish, which likely increases the stock yield (see, e.g., Kvamme and Bogstad 2007). This maneuver could be done by increasing the current 44 cm minimum landing size, increasing the current 130 mm trawl mesh size, and shifting the fisheries towards other gears such as gillnets (ICES 2020). At present, about 70% of total catches (and almost all catches by other countries than Norway) are taken by trawl (ICES 2020), so such changes will inevitably have large political ramifications. It should also be considered that other species such as haddock, which do not grow to such large body lengths as cod do, are fished together with cod (ICES 2020). Hence, increasing mesh size may, in such cases, reduce the yield from other stocks, although there are no model studies on this topic. NEAC is the most important species for Norwegian fisheries and the Russian fleet fishing in the Barents Sea (Gullestad et al. 2020). Therefore, reducing the yield from this stock to rebuild the stock’s maturity structure, which likely will take many decades, would be a tough sell to managers and the industry. Fortunately, the stock’s age and size structure are now largely rebuilt and resemble the structure in the late 1940s when the stock had just experienced a period of low fishing mortality due to WWII (Hylen 2002; ICES 2020). Hutchings and Kuparinen (2020) concluded that since many FIE-implicated populations have recovered rapidly to management-based targets following cessation of overfishing, FIE is generally of minor importance to recovery compared to overfishing, the magnitude of depletion and natural mortality. So, avoiding overfishing and reducing fishing pressure when population abundance becomes low remains a key management strategy, particularly on immature fish.
Our analysis of the maturation traits of NEAC had the advantage of data from more than eight decades. The large variations in A50, the explanatory variables TSB, KolaT and the size-selective fishery (F5-8) certainly helped to detangle the impacts of the various variables. Regardless, there is an implicit assumption in linear regression analyses on time series data that there is no correlation between consecutive residuals. This assumption is often not fully met in biological time series. However, the available time series was long enough to overcome these obstacles in the present analysis of NEAC. The impact of fisheries-induced adaptive change (FIAC) is statistically significant at the P < 0.001 level. FIAC is our primary explanation for why the present value of age at 50% maturity (A50) has not recovered to its pre-war level. The annual accumulation of fishing mortalities appears potentially to be a strong indicator of the long-term impact of FIAC. Phenotypic plasticity is in part covered in our model set-up by the inclusion of TSB and KolaT as factors, but only TSB showed a clear impact on the maturity curve. The model includes the option that the average fishing mortality on the 5- to 8-year-olds (F5-8) needs to be above F_bal before FIAC is detectable in model runs. However, the possible interval of F_bal for NEAC A50 is broad, ranging from 0.00 to 0.40, whereas ecogenetic models indicate a F_bal value near 0. However, FMSY on ages 5–8 years appears close to the upper end of the F_bal range. The recovery time of A50 depends critical on F_bal. Pending on the intensity of fishing mortalities on immatures, the development of A50 in the next decades may help narrow down the possible range in F_bal. Given the uncertainty of the value of F_bal, we do not propose this to be a management target, but it is a concept that may have practical consequences and is worth pursuing.
We thank the Associate Editor and two anonymous reviewers for comments, which greatly improved the content of our manuscript. The Russian Federal Research Institute of Fisheries and Oceanography (VNIRO), Polar branch in Murmansk, Russia, kindly provided us with the Kola Transect sea temperature data. Terje Jørgensen and Knut Korsbrekke, both at Institute of Marine Research, were most helpful during the initial preparation phase of this study. The writing-up of this article was supported by the Norwegian Research Council (NFR) projects Scaling Climate Effects from Individual Physiology to Population Responses (ScaleClim; No. 268336; OSK) and Cod: DIet and food web dyNAmics (CoDINA; No. 255460; GO). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Ajiad A., Jakobsen T., and Nakken O. 1999. Sexual difference in the maturation of Northeast Arctic cod. J. Northwest Atl. Fish. Sci. 25: 1–15.
Bangera R., Drangsholt T.M.K., Nielsen H.M., Sae-Lim P., Ødegard J., Puvanendran V., et al. 2015. Genotype by environment interaction for growth in Atlantic cod (Gadus morhua L.) in four farms of Norway. J. Mar. Sci. Eng. 3: 412–427.
Bogstad B., Yaragina N.A., and Nash R.D.M. 2016. The early life-history dynamics of Northeast Arctic cod: levels of natural mortality and abundance during the first 3 years of life. Can. J. Fish. Aquat. Sci. 73(2): 246–256.
Eikeset A.M., Dunlop E.S., Heino M., Storvik G., Stenseth N.C., and Dieckmann U. 2016. Roles of density-dependent growth and life history evolution in accounting for fisheries-induced trait changes. Proc. Natl. Acad. Sci. U.S.A. 113: 15030–15035.
Gullestad P., Sundby S., and Kjesbu O.S. 2020. Management of transboundary and straddling fish stocks in the Northeast Atlantic in view of climate-induced shifts in spatial distribution. Fish Fish. 1: 1008–1026.
Heckwolf M.J., Meyer B.S., Hasler R., Hoppner M.P., Eizaguirre C., and Reusch T.B.H. 2020. Two different epigenetic information channels in wild three-spined sticklebacks are involved in salinity adaptation. Sci. Adv. 6: eaaz1138.
Kjesbu O.S., Kryvi H., Sundby S., and Solemdal P. 1992. Buoyancy variations in eggs of Atlantic cod (Gadus morhua L.) in relation to chorion thickness and egg size: Theory and observations. J. Fish Biol. 41: 581–599.
Kjesbu O.S., Bogstad B., Devine J.A., Gjøsæter H., Howell D., Ingvaldsen R.B., et al. 2014. Synergies between climate and management for Atlantic cod fisheries at high latitudes. Proc. Natl. Acad. Sci. U.S.A. 111: 3478–3483.
Marshall C.T., Kjesbu O.S., Yaragina N.A., Solemdal P., and Ulltang Ø. 1998. Is spawner biomass a sensitive measure of the reproductive and recruitment potential of Northeast Arctic cod? Can. J. Fish. Aquat. Sci. 55(7): 1766–1783.
Olsen E., Lilly G.R., Heino M., Morgan M.J., Brattey J., and Dieckmann U. 2005. Assessing changes in age and size at maturation in collapsing populations of Atlantic cod (Gadus morhua). Can. J. Fish. Aquat. Sci. 62(4): 811–823.
Ottersen G., Hjermann D.Ø., and Stenseth N.C. 2006. Changes in spawning stock structure strengthen the link between climate and recruitment in a heavily fished cod (Gadus morhua) stock. Fisheries Oceanogr. 15: 230–243.
Pinsky M.L., Eikeset A.M., Helmerson C., Bradbury I.R., Bentzen P., Morris C., et al. 2021. Genomic stability through time despite decades of exploitation in cod on both sides of the Atlantic. Proc. Natl. Acad. Sci. U.S.A. 118(15): e2025453118.
Scott B.E., Marteinsdottir G., Begg G.A., Wright P.J., and Kjesbu O.S. 2006. Effects of population size/age structure, condition and temporal dynamics of spawning on reproductive output in Atlantic cod (Gadus morhua). Ecol. Model. 191: 383–415.
Shelton A.O., Hutchings J.A., Waples R.S., Keith D.M., Akcakaya H.R., and Dulvy N.K. 2015. Maternal age effects on Atlantic cod recruitment and implications for future population trajectories. ICES J. Mar. Sci. 72: 1769–1778.
Stige L.C., Yaragina N.A., Langangen O., Bogstad B., Stenseth N.C., and Ottersen G. 2017. Effect of a fish stock’s demographic structure on offspring survival and sensitivity to climate. Proc. Natl. Acad. Sci. U.S.A. 114: 1347–1352.
Therkildsen N.O., Hemmer-Hansen J., Als T.D., Swain D.P., Morgan M.J., Trippel E.A., et al. 2013. Microevolution in time and space: SNP analysis of historical DNA reveals dynamic signatures of selection in Atlantic cod. Mol. Ecol. 22: 2424–2440.
Trippel, E.A., Kjesbu, O.S., and Solemdal, P. 1997. Effects of adult age and size structure on reproductive output in marine fishes. In Early life history and recruitment in fish populations. Edited by R.C. Chambers and E.A. Trippel. Chapman and Hall, London. pp. 31–62.
Yaragina, N.A., Aglen, A., and Sokolov, K.M. 2011. Cod. In The Barents Sea, ecosystem, resources, management. Half a century of Russian–Norwegian cooperation. Edited by T. Jakobsen and V.K. Ozhigin. Tapir Academic Press, Trondheim, Norway. pp. 225–270.
Carl Jakob Rørvik, Bjarte Bogstad, Geir Ottersen, and Olav Sigurd Kjesbu. 2021. Long-term interplay between harvest regimes and biophysical conditions may lead to persistent changes in age at sexual maturity of Northeast Arctic cod (Gadus morhua). Canadian Journal of Fisheries and Aquatic Sciences.
79(4): 576-586. https://doi.org/10.1139/cjfas-2021-0068
If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.
1. Growth portfolios buffer climate‐linked environmental change in marine systems