Demographics and recovery potential of exploited marine teleosts

TL;DR

This study examines the demographics of 77 marine fish stocks and finds that recovery potential varies more by region than by species traits.

Contribution

The study applies conservation biology metrics to marine fish stocks, revealing the importance of regional factors over life history traits in predicting recovery.

Findings

Median population growth rates (λ) were close to 1, with only slight increases when fishing was excluded.

Population growth for many stocks relied on a few years of high growth rates.

Recovery potential was strongly influenced by regional factors rather than species-specific life history traits.

Abstract

Equilibrium concepts and the expectation of compensatory density dependence remain fundamental to fisheries science, but stock collapses and an increasing appreciation of environmental factors have raised questions about their real-world applicability. To explore the demographic variability of harvested marine fishes, we have calculated metrics commonly used in conservation biology to describe the demographics for 77 assessed stocks from the North Atlantic and Northeast Pacific Oceans using life-tables. We found that median annual population growth rates (λ) were centered around 1, and surprisingly, they were only slightly higher when the effect of fishing was excluded. For most stocks, as abundance declined, λ tended to increase and become more variable as would be expected from compensatory dynamics. The population growth of several stocks was sustained by a limited number of years with exceptionally high rates. However, the ability of a stock to increase from low abundance appeared largely independent of life history characteristics and exhibited stronger geographical differences among stocks of the same species (notably Atlantic cod). Life history characteristics alone were poor predictors of annual population growth or future recovery potential, whereas regional factors appeared to be more influential. Overall, recovery potential remained relatively high, with simulations indicating that 62 of the stocks would be highly likely to double in size within 20 years in the absence of fishing. Low recovery potential was exclusively observed in stocks with a low median λ and low variability in λ. These results suggest that understanding stock-specific (rather than species-specific) demographic parameters is necessary to promote sustainable management or develop rebuilding plans for collapsed stocks.

Figures

Locations of the stocks used in this analysis.Colors indicate the number of stocks within each subarea. In the Atlantic the subareas were divided using the ICES and NAFO subareas. The ICES region was grouped into seven subareas which broadly aligned with the ecoregions definitions developed by ICES. The NAFO subareas were divided into two regions, the first combining NAFO Regions 2J, 3, and 4 with the second combining NAFO regions 5 and 6(A-C). In the Pacific, the subareas were divided into a Gulf of Alaska subarea and a subarea covering the Bering Sea and Aleutian Islands. The divisions were made using the Alaska Department of Fish and Game Commercial Fisheries Groundfish Registration Areas.

Boxplot of theλ estimates for each stock of G. morhua, M. aeglefinus, P. virens, C. harengus, P. platessa, and S. solea. Blue are the λdem estimates and grey are the λreal estimates.The red stars indicate the median and the dashed vertical line represents a λ of one. The x axis is on the log scale.

Boxplot of the λ estimates for the remainder of the stocks in this analysis, plotted by the location of the stock.The abbreviated species name is provided on the y axis. Blue are the λdem estimates and grey are the λreal estimates. The red stars indicate the median and the dashed vertical line represents a λ of one. The x axis is on the log scale.

The average standard deviation of λdem for Atlantic herring (C. harengus), Atlantic cod (G. morhua), haddock (M. aeglefinus), and pollock (*P. virens) stocks in the western Atlantic (NAFO regions 2-5), Northern Atlantic (Iceland, Faroe’s, and the northeast Arctic), and Eastern Atlantic (Irish, Celtic, North, and Baltic Seas). The error bars represent one standard deviation.

Density dependence of λ.Panel (a) compares the GAM smooths for λdem and λreal and the GAM smooth and the GAM estimated 95% confidence interval. Panel (b) shows λdem against the proportion of maximum abundance. The points represent the estimate for each year and stock along with the GAM smooth and the GAM estimated 95% confidence interval. Panel (c) shows λreal against the proportion of maximum abundance. Panel (d) shows the data and GAM smooths for two Atlantic cod (G. morhua) stocks in the Northeast Atlantic and two stocks in the Northwest Atlantic.

𝐑0 Estimates for each stock of G. morhua, M. aeglefinus, P. virens, C. harengus, P. platessa, and S. solea.The blue values are calculated using natural mortality and grey includes both natural mortality and fishing mortality. The red stars indicate the median and the dashed vertical line represents a lifetime reproductive output of one. The x axis is on the log scale.

Generation time estimates for each stock of G. morhua, M. aeglefinus, P. virens, C. harengus, P. platessa, and S. solea.The blue values are calculated using natural mortality and grey includes both natural mortality and fishing mortality. The red stars indicate the median.

Proportion of the simulations in which a given stock has doubled in size when there is no fishing mortality.Stocks are grouped by taxonomic order, with 39 Gadiform stocks, 2 Scombriform stocks, 1 Scorpaeniform stock, 20 Pleuronectiform stocks, 11 Clupeiform stocks, and 4 Perciform stocks. The x-axis represents the year of the simulation, with the first 20 years shown.

The doubling time recovery potential, the colors represent low (grey), moderate (white), or high(blue) for each stock after 20 years compared to the medianλdem and the standard deviation of λdem. The red dashed line on the left panel is a λdem value of one.

Equations

Full text

Introduction

Marine fishes were once believed to be inexhaustible and immune to the effects of fishing [1–2]. However, striking collapses of once-prolific fisheries and their delayed or failed recovery [3–5] have contradicted this perception. The assumption of resilience originated from early theoretical population dynamics modeling (e.g., [6]) that culminated in the development of modern fisheries science, largely built around the concept of maximum sustainable yield (MSY, [7–9]). Despite long-standing concerns regarding the use of MSY [10,11], its theoretical underpinning in equilibrium concepts and the associated ideas that annual per-capita reproductive output and/or juvenile survival rates increase as abundance declines have persisted in fisheries science [8].

Indeed, many modern assessments incorporate a stock-recruit function in which per-capita recruitment increases as biomass declines [7,9,12,13]. Such fisheries assessments implicitly assume that lifetime reproductive rates in the absence of fishing should be greater than one [12,14] and that even heavily-exploited fish stocks retain an innate capacity for recovery [15]. Although MSY and equilibrium concepts are relatively straightforward, intuitive, and provide a means of developing tractable mechanistic models, questions remain regarding their applicability [10,11]. There are concerns about whether they accurately reflect the system being modeled [12,16,17], particularly for species with markedly different life-history characteristics [18].

Life-history characteristics (e.g., longevity, maximum size, age at maturity) are often considered indicative of a species’ natural mortality rate, reproductive capacity and overall productivity [19–21], and are increasingly used within data-limited assessment approaches to provide management advice [22,23]. Ecological theory suggests that ‘r-selected’ species characterized by small body size, numerous offspring, rapid growth, early maturity and short lifespans, should have greater recovery potential and be better able to withstand fishing pressure [24,25]. Conversely, “K-selected” species with life-histories characterized by longer lifespans and later maturity, or those with lower per-capita reproductive output, are considered less resilient to fishing pressure and are expected to have a lower recovery potential [21]. The life-history traits of ‘r-selected’ species should manifest as a high innate capacity for population growth in the absence of fishing, resulting in rapid recovery from low abundance [26]. Furthermore, owing to a shared evolutionary history, the life-history traits of closely related species are anticipated to be more similar than those of distantly related species [27]. This implies that there should be less heterogeneity in the productivity among stocks of a single species than between stocks of different species with similar life-history traits. However, it remains unclear if correlated and compounding ecological and environmental processes (e.g., climate change effects, predator-prey dynamics) can obscure the importance of life-history on realized productivity. This raises concerns that evolutionary history and life-history characteristics may serve as a weak foundation for estimating demographic rates of species for which this information is not available.

Stock assessment data from commercially-harvested marine fishes provide a unique opportunity to explore demographic variability and realized compensatory responses. Focusing on species with a long assessment history means that their life-history characteristics are likely to be well-known and that it is possible to quantify realized variability in survival and reproductive output while accounting for the impact of fishery removals. We based our analyses on 77 marine fish stocks included in a database of age-specific data [28]. By constructing life-tables (i.e., age-specific estimates of survival and reproduction rates) from abundance-at-age and catch-at-age data, we calculated standard demographic metrics, including annual population growth rates, doubling time, lifetime reproductive success, and generation time [29]. Our objectives are to: 1) quantify the potential for population growth ( $[eqn]$ ) and the strength of density dependence for each stock, 2) quantify recovery potential of each stock, and 3) explore the variability of these demographic rates both regionally and taxonomically. Our results have implications for the sustainable management of these commercially-harvested marine fishes and could be used to inform future management actions for stocks that have collapsed.

Methods

Data

These analyses were conducted using stock assessment data from a repository which included sufficient data for 77 stocks from the Northeast Atlantic, Northwest Atlantic, and the North Pacific; including stocks managed by International Council for the Exploration of the Sea (ICES; 44 stocks), NOAA (NOAA;25 stocks), and Fisheries and Oceans Canada [DFO, 8 stocks, Fig 1; this repository has been updated with the data used in this analysis [28,30]. Data were available between 1947 and 2018, with an average time series length of 38.6 years (range 13–68 years). The assessed stocks represented 6 taxonomic orders, with the majority of stocks representing Gadiformes (39 stocks), Pleuronectiformes (20 stocks) and Clupeiformes (11 stocks). The other three orders had relatively few stocks included in the analysis, with 4 Perciform stocks, 2 Scombriform stocks, and just 1 Scorpaeniform stock. The majority of the data were for stocks from six species: Atlantic cod (Gadus morhua; 14 stocks), herring (Clupea harengus; 9 stocks), haddock (Melanogrammus aeglefinus; 8 stocks), pollock (Pollachius virens; 5 stocks), common sole (Solea solea; 5 stocks), and European plaice (Pleuronectes platessa; 4 stocks).

Locations of the stocks used in this analysis.Colors indicate the number of stocks within each subarea. In the Atlantic the subareas were divided using the ICES and NAFO subareas. The ICES region was grouped into seven subareas which broadly aligned with the ecoregions definitions developed by ICES. The NAFO subareas were divided into two regions, the first combining NAFO Regions 2J, 3, and 4 with the second combining NAFO regions 5 and 6(A-C). In the Pacific, the subareas were divided into a Gulf of Alaska subarea and a subarea covering the Bering Sea and Aleutian Islands. The divisions were made using the Alaska Department of Fish and Game Commercial Fisheries Groundfish Registration Areas.

The repository included age-specific estimates of abundance, weight (kg), biomass (tonnes), maturity, natural mortality, and fishery removals (tonnes) for each stock (Table A1). While uncertainty is inherent in assessment processes, each assessment is the basis for management advice. As such, the underlying data represents the best available understanding of each stock’s population dynamics at the time of the assessment.

Population growth rates

We used life-table reconstructions to calculate survival-at-age and fecundity-at-age over time from the available data series. We applied the Euler-Lotka equation [31] to estimate the corresponding population growth rate ( $[eqn]$ ) in a given year. This allowed us to use information on recruitment, mortality, and maturity to describe each stock’s demographics and how they changed with age and time using a tractable and consistent model.

The discrete time Euler-Lotka equation was first proposed by Euler [later published by 32] and is a special case of the continuous time demographic model subsequently developed by Lotka [31,33,34]:

$$[eqn]$$

Where $[eqn]$ is the age, $[eqn]$ is the age at first reproduction, $[eqn]$ is the maximum age, $[eqn]$ is the discrete rate of population growth ( $[eqn]$ , where $[eqn]$ is the intrinsic rate of increase), $[eqn]$ is the survivorship (proportional) to age $[eqn]$ , and $[eqn]$ is the fecundity of individuals of age $[eqn]$ . This equation is solved for each year ( $[eqn]$ ) with available data. The discrete-time Euler-Lotka equation is the characteristic polynomial of the Leslie Matrix, where $[eqn]$ is the fecundity of the first age class with mature individuals [29,35]. We calculated survivorship-at-age and fecundity-at-age using the available stock assessment data, which enabled $[eqn]$ to be estimated. This assumes that each stock represents a discrete population with independent dynamics.

When calculating survival and fecundity inputs, we treated the last observed age class as the terminal age, beyond which there was no contribution to productivity. The first age class was age 1 for the majority of stocks (N = 42), where only 21 of the stocks had a first age of 2 (N = 10) or greater (N = 11). For simplicity, we assumed survivorship to the first age class represented by the data to be 1, which can result in a slight overestimation of productivity from the Euler-Lotka equation [36]. However, we tested the influence of this assumption by back-calculating the abundance for the earliest ages using the available time series of age-specific natural mortality rates, and this resulted in negligible differences to the annual productivity estimates (Figs S10 and S11).

Fecundity-at-age.

We parameterized fecundity as the number of offspring (recruits) produced by each female in each age class rather than as the number of eggs, because egg production by size, fertilization rates, and egg-to-age-0 survival rates are unknown. We first calculated an overall recruitment rate as the number of recruits produced per kilogram of spawning biomass $[eqn]$ :

$$[eqn]$$

Where $[eqn]$ is the number of individuals in the youngest age class available for a specific stock (hereafter referred to as recruits), and $[eqn]$ is the total spawning stock biomass in year $[eqn]$ , offset by the recruit age ( $[eqn]$ ). For example, if the first age class in the abundance time series is 2 years, we assumed that the recruits in 2010 were produced by the spawning stock biomass in 2008 ( $[eqn]$ ).

Next, we calculated the total number of recruits produced by each age class ( $[eqn]$ ) by multiplying the spawner biomass in each age class $[eqn]$ by the number of recruits produced per kilogram of spawning stock biomass:

$$[eqn]$$

The reproductive success of fish is expected to vary by age [37–39]. However, due to a lack of available evidence to quantify these effects, we made the simplifying assumption that each kilogram of spawning biomass contributed equally to the total reproductive output, regardless of age. This assumption would have a negligible impact on the overall $[eqn]$ estimates, but would have an effect on the relative contribution of the underlying demographic parameters (i.e., the fecundity and mortality contributions). We calculated the age-specific fecundity ( $[eqn]$ ) as the number of offspring produced by each spawner in the stock. The number of spawners in year $[eqn]$ and age class $[eqn]$ ( $[eqn]$ ) is the proportion of mature individuals in each age class ( $[eqn]$ ) multiplied by the number of individuals in each age class:

$$[eqn]$$

where annual fecundity becomes:

$$[eqn]$$

Fecundity is offset by one year, as the recruits observed in year $[eqn]$ enter into the Euler-Lotka formulation in year $[eqn]$ .

Survivorship-to-age.

The number of individuals surviving to each age (survivorship-to-age) depends on age-specific mortality rates (i.e., survival-at-age). Commercially exploited fishes experience fishing mortality ( $[eqn]$ ) once they grow large enough to become vulnerable to the fishing gear [40]. Fishing mortality ( $[eqn]$ ) and natural mortality ( $[eqn]$ ) in combination yield the total mortality ( $[eqn]$ ) affecting each age class in each year. Due to challenges with disentangling natural mortality and fishing mortality (known as separability) different assumptions are made about natural mortality during stock assessments. For example, natural mortality may be assigned a different value for each age (e.g., Bering Sea Gadus chalcogrammus), allowed to vary over time and by age (e.g., NAFO 5–6 GB Gadus morhua), input as a constant across ages (e.g., NAFO 5–6 Pollachius virens), or specified via another method. These natural mortality estimates are derived within the assessment model or estimated empirically using various methods [41]. We did not explore various ways of calculating M and used the estimates of natural mortality by age and year available in the data repository (representing outputs from each assessment), which we converted into annual rates at-age for input into the Euler-Lotka equation. Both natural mortality and fishing mortality were calculated as proportional values.

The data also contained the number of fish harvested by the fishery each year and for each age class. To estimate fishing mortality ( $[eqn]$ ) for each age class in each year we divided the catch-at-age data by the assessed numbers at-age.

$$[eqn]$$

Together, fishing mortality ( $[eqn]$ ) and natural mortality ( $[eqn]$ ) provide a total mortality estimate ( $[eqn]$ ) for each age in each year.

$$[eqn]$$

We defined survival-at-age as:

$$[eqn]$$

Next, we used survival-at-age to calculate cumulative survivorship-to-age used as an input for the Euler-Lotka equation:

$$[eqn]$$

The availability of catch-at-age time series in the data repository allowed for the separation of fishing and natural mortality rates, enabling estimates of $[eqn]$ and $[eqn]$ (as above). As a result, we were able to calculate two $[eqn]$ values, the realized $[eqn]$ from the reconstruction that accounted for $[eqn]$ , $[eqn]$ , and fecundity, and a theoretical $[eqn]$ value using the $[eqn]$ and fecundity estimates from the reconstruction. Given this second $[eqn]$ value excluded the impact of fishing, and we refer to this as the demographic lambda ( $[eqn]$ ). We recognize that the age structure in the underlying time series is still influenced by fishing pressure, and therefore $[eqn]$ does not represent the maximum intrinsic rate of growth in the absence of fishing [42]. The value is the potential growth of the population if the ‘surplus production’ was not harvested in a given year. Thus, $[eqn]$ represents the productivity that the stock could have achieved in each specific year, if there had not been fishing.

Estimation.

It is not straightforward to propagate observation error in each time series when calculating survivorship-to-age and fecundity-at-age. To ensure that variability did not unduly affect predictions for annual stock productivity, we fit the Euler-Lotka model using an iterative minimization process, where predicted abundance in the following year:

$$[eqn]$$

was compared to the abundance estimate from the stock assessment for the corresponding year ( $[eqn]$ ).

Where necessary, we adjusted the annual age-specific fecundity and natural mortality estimates until the $[eqn]$ estimate from the Euler-Lotka equation resulted in an abundance time series estimate ( $[eqn]$ ) that was within 5% of the stock assessment abundance time series (median difference was −0.1%). This adjustment was performed iteratively, with a step change (for most stocks each step was a change of 0.5%, but 0.05% was used for some stocks) first applied to the natural mortality estimates at each age. If the estimated abundance was not within 5% of the observed abundance estimate, then a step change in the fecundity was applied to each age. This procedure was repeated until the estimated abundance was within 5% of the observed abundance.

This approach ensured that our analyses reflected the same understanding of stock productivity, natural mortality and fecundity as the current assessments, by effectively treating the abundance estimates in the data repository as the observations. Although there are inherent limitations to treating stock assessment outputs as observed data [43], these assessments represent the best available knowledge on the realized population dynamics of each stock. Thus, abundance time series from the stock assessments are an appropriate benchmark against which to compare.

The minimization procedure resulted in the final estimates of $[eqn]$ being approximately 6.7% smaller than the initial estimates. Population growth rates from life-table analyses typically decline when variability in the input parameters is accounted for (e.g., [44]). However, it is also possible that the removals time series used in the stock assessments did not capture all sources of fishing mortality (e.g., unreported catch, bycatch and incidental mortality, [45], [46]). This would result in the estimates of $[eqn]$ being lower than reality, necessitating slight reductions in fecundity and/or increases in natural mortality during the iterative minimization process.

Density dependence

Life-table analyses make no assumptions about population regulation and thus do not incorporate a density-dependent relationship between spawners and recruits [47]. However, our estimated productivity parameters could be explored for evidence of density dependence by modelling each stock’s annual population growth rate against relative stock size (expressed as a proportion of maximum abundance) in the same year. Density independence occurs when $[eqn]$ shows no relationship with abundance. Negative density dependence is observed when $[eqn]$ declines as abundance increases; this is the type of compensatory response described by typical spawner-recruit relationships [12]. Positive density dependence occurs when $[eqn]$ declines with declining abundance, which is also called an Allee effect or depensation when it occurs at low abundance [14,48–50]. We evaluated both the realized productivity estimates ( $[eqn]$ ) and the theoretical demographic productivity values excluding the impact of fishing mortality ( $[eqn]$ ) for evidence of density dependence.

We used Generalized Additive Mixed Models [51,52] to describe potential relationships between $[eqn]$ and abundance:

$$[eqn]$$

Where $[eqn]$ is the proportion of the maximum observed abundance for the $[eqn]$ observation. The stock ( $[eqn]$ ) is assumed to be a random effect, which is normally distributed with a mean of 0 and a variance of $[eqn]$ . The response variable $[eqn]$ is assumed to follow a normal distribution on the log scale. The $[eqn]$ by $[eqn]$ interaction was modeled using the fs basis in the 'mgcv' package in R (version 4.5). This formulation treats the smooth curves as random effects; it is used when a factor has a large number of levels, and it enables each stock to have a unique smooth [52,53]. We explored various model formulations and alternative error distributions; however, all models had similar challenges overcoming either the underlying non-normality or the heteroskedasticity of the data. Given these challenges, we did not use the model to determine the statistical significance of the observed trends (e.g., we did not evaluate if the smooth terms are significant). This approach was adopted to capture the global relationship and provided a relative metric of variability among the stocks. We compared values of $[eqn]$ estimates for each stock at the minimum, maximum, and 40% of the maximum abundance; the latter value was chosen based on both theoretical expectations and observed population dynamics [e.g., 50,51,54]

For a life-table reconstruction such as this, the maximum value of $[eqn]$ is constrained at higher abundances because the maximum abundance is a known value (e.g., the maximum possible $[eqn]$ at 50% of the maximum abundance is 2, at 25% it is 4, and at 100%, it is 1). Conversely, the minimum value of $[eqn]$ is constrained at lower abundances (e.g., $[eqn]$ must be $[eqn]$ 1 when the stock is at its minimum observed abundance). As a result, at the maximum abundance $[eqn]$ will be $[eqn]$ 1, at the minimum $[eqn]$ will be $[eqn]$ 1, and higher variability in $[eqn]$ values at low abundances should be relatively common.

Recovery potential

Numerous demographic parameters can be derived from the outputs of the life-table reconstruction [29]. Here, we focus on three that are strongly linked to recovery potential: lifetime reproductive success, generation length, and doubling time [55].

Lifetime reproductive success.

Lifetime reproductive success, also called the net reproductive rate, describes the average number of female spawners produced by one female spawner throughout her lifetime. Values greater than 1 suggest a stock increased, values less than 1 suggest a stock declined, and a value of 1 indicates that the average female has produced sufficient offspring over her lifetime to replace herself [29]. Lifetime reproductive success ( $[eqn]$ ) was calculated for each cohort in which fecundity and natural mortality were available for complete cohorts (i.e., cohorts in which all available age classes are represented), assuming a 50:50 sex ratio.

$[eqn]$ was calculated both by using only the natural mortality component ( $[eqn]$ ) and with the combined effect of natural mortality and fishing mortality ( $[eqn]$ ), respectively. This provided a measure of the reproductive potential from observed stock characteristics alone, and how it changed with fishing pressure, recognizing that setting $[eqn]$ does not fully remove the effect of fishing pressure on the cohort. $[eqn]$ was calculated for each cohort as:

$$[eqn]$$

This formulation ignores any contribution from a plus group (i.e., the effect of individuals who are older than the oldest age in the life-table) by assuming that the last observed age in the data represents the longevity of the cohort.

Generation length.

Generation length is an alternate indicator for the relative reproductive rate of a stock, where shorter generation lengths are associated with faster population growth. Here, we calculated generation length for each cohort as the mean age of spawners that produced recruits [29]. Values of both the lifetime reproductive rate and generation length were derived for each complete cohort, which meant that we could also evaluate changes in reproductive capacity over time for each stock. Additionally, we compared theoretical reproductive potential and generation lengths from life history characteristics alone ( $[eqn]$ ) with those expected under current fishing mortality ( $[eqn]$ ).

We defined generation length as the mean age of the spawners that produced the observed recruits in a cohort ( $[eqn]$ , [29]).

$$[eqn]$$

Doubling time.

Comparing the median time frame over which abundance can double (i.e., doubling time) in the absence of fishing provides an indication of relative productivity among species and stocks. We applied a Monte Carlo approach that draws from the observed distribution of fecundity and natural mortality in a projection simulation to estimate the doubling time for each stock.

The annual fecundity and natural mortality values varied considerably across stocks in this study. Therefore, we used a simulated stock projection to account for this variability when evaluating recovery potential. Across 1000 simulations, we took a random sample of the annual values for fecundity and natural mortality, used the life-table reconstruction to estimate $[eqn]$ , and calculated abundance in the following year. Abundance of each stock was projected over 100 years, and the proportion of the simulations in which the stock size had doubled from the initial abundance was recorded each year. We considered stocks to have high recovery potential if 75% of simulated trajectories doubled after 10 years, medium potential if 50–74% of simulations doubled after 20 years, and low recovery potential if < 50% of simulated trajectories doubled after 20 years. For simplicity, we did not incorporate density-dependent effects or any correlation between fecundity and natural mortality. Once a population had doubled, it was considered to have remained at that abundance or greater for the remainder of the projection. This indicates that the metric reflects the earliest possible doubling time for each simulation in the absence of fishing mortality.

Results

Population growth rates (λ)

The median population growth rate when excluding the effects of fishing ( $[eqn]$ ) for the 77 stocks in this analysis was 1.091, representing an innate capacity for annual population growth of approximately 8.7%. Aggregating over years, the stock with the highest median theoretical population growth rate was Atlantic cod in the Irish Sea (median $[eqn]$  = 1.44) and the lowest was pollock in the Gulf of Alaska (median $[eqn]$  = 0.89; Table 2). When the effects of fishing were accounted for, the overall median population growth rate ( $[eqn]$ ) dropped to 0.971, corresponding to overfishing by 2.9% annually (Figs 2 and 3).

Table 2: A comparison of the λdem and λreal estimates for stock at the minimum and maximum abundance, along with the estimate at 40% of the maximum abundance. Values greater than 1 indicate that the population would tend to increase in size, while values less than 1 indicate the population would tend to decline in size.

Boxplot of theλ estimates for each stock of G. morhua, M. aeglefinus, P. virens, C. harengus, P. platessa, and S. solea. Blue are the λdem estimates and grey are the λreal estimates.The red stars indicate the median and the dashed vertical line represents a λ of one. The x axis is on the log scale.

Boxplot of the λ estimates for the remainder of the stocks in this analysis, plotted by the location of the stock.The abbreviated species name is provided on the y axis. Blue are the λdem estimates and grey are the λreal estimates. The red stars indicate the median and the dashed vertical line represents a λ of one. The x axis is on the log scale.

The realized range of population growth rates was much narrower, and the median $[eqn]$ was below 1 more than 50% of the time in 85.7% of the stocks (Figs 2 and 3). However, many of the stocks included in our analysis experienced infrequent large recruitment events that resulted in high $[eqn]$ estimates (Figs 2 and 3). For example, while the Gulf of Alaska pollock stock had the lowest median $[eqn]$ (0.89), it had a mean $[eqn]$ of 1.22, and when the effects of fishing were included, this stock had the $[eqn]$ highest mean $[eqn]$ (1.19, Figs 2 and 3).

Mean $[eqn]$ was lower than the median for just 8 of the 77 stocks, and the mean $[eqn]$ was less than one for 15 of the 77 stocks, indicating that years of high population growth resulted in a net positive productivity for most stocks. The majority of the stocks with a mean $[eqn]$ value less than one were in the central or NW Atlantic (N = 8), and included 5 of the 14 Atlantic cod (G. morhua) stocks.

Several species are harvested throughout the North Atlantic, including two relatively large and long-lived Gadidae species: Haddock (Melanogrammus aeglefinus) and Atlantic cod (Gadus morhua). These two species exhibited some of the highest variability in population growth rates both over time and among stocks (Fig 2). For both species, the variability in the population growth rate declined moving from the Eastern Atlantic to the Western Atlantic (Fig 4). The pattern was most pronounced for Atlantic cod, where stocks in the Eastern Atlantic rarely experienced high values of $[eqn]$ (Figs 2 and 4). Pollock (Pollachius virens) is another Gadidae species harvested across the North Atlantic. The variability in $[eqn]$ was lower for this species throughout its range (Figs 2 and 4). Atlantic herring (Clupea harengus) is a relatively small fast-growing species also found throughout the North Atlantic. However, there were no clear regional patterns in $[eqn]$ for this species and its variability in $[eqn]$ was lower than observed in several Atlantic cod and haddock stocks.

*The average standard deviation of λdem for Atlantic herring (C. harengus), Atlantic cod (G. morhua), haddock (M. aeglefinus), and pollock (P. virens) stocks in the western Atlantic (NAFO regions 2-5), Northern Atlantic (Iceland, Faroe’s, and the northeast Arctic), and Eastern Atlantic (Irish, Celtic, North, and Baltic Seas). The error bars represent one standard deviation.

Two relatively short-lived and fast-growing species (Sand eel [Ammodytes tobianus] and Norway pout [Trisopterus esmarkii]) in the North Sea exhibited high variability in $[eqn]$ , though this variability was smaller than that observed in the haddock and Atlantic cod stocks in the area (Fig 3). For the remainder of the other stocks in this analysis, years in which the value of $[eqn]$ were elevated (e.g., above $[eqn]$ 1.5) were rarer, but no clear taxonomic or regional patterns were evident, with whiting (Merlangus merlangus) in the Celtic Sea, summer flounder (Paralichthys dentatus) and scup (Stenotomus Chrysops) in the Western Atlantic, and Gulf of Alaska pollock (Gadus chalcogrammus) all having moderate variability in $[eqn]$ (Fig 3).

There was no clear temporal trend in the lambda estimates, however $[eqn]$ in the 2000s was low (1.08), and in the 2010s the standard deviation of $[eqn]$ was higher than in the other decades, resulting in a low median $[eqn]$ (1.07), but a mean $[eqn]$ estimate similar to what was observed before the 2000s (1.23; Table 1). The high standard deviation in the 2010s was somewhat impacted by a single stock (CS RA M. aeglefinus). When this stock was removed from the analysis the standard deviation remained the highest in the time series, but was not as extreme (sd of $[eqn]$ declined to 0.63), while the median and mean $[eqn]$ ’s remained similar to values observed before the 2000s (e.g., the median and mean $[eqn]$ were 1.07 and 1.2 respectively). This analysis also indicated that relative abundances (proportion of maximum observed) of the stocks have been in decline over the last 40 years, with the lowest values observed in the 2010s (Table 1).

Table 1: A comparison of the median, mean, and standard deviation of λdem and λreal by decade for all the stocks in this analysis. The table includes the number of stocks included in each decade (N) and the median proporiton of the maximum observed abundance (Prop of max).

Density dependence

Overall, both $[eqn]$ and $[eqn]$ increased as abundance declined. The changes in the rate of population growth were more pronounced at low and high abundance levels. The relationship between $[eqn]$ and the proportion of maximum abundance was relatively flat when abundance ranged between 40–70% of maximum (Fig 5). The $[eqn]$ remained above one when stocks approached their maximum observed abundance in more than 87% of the stocks (Tables 2 and 3 and S1 and S2 Figs). However, $[eqn]$ dropped below one when the stock abundance exceeded approximately 33% of the maximum abundance. At lower abundances, the model estimated $[eqn]$ increased as abundance declined more slowly than $[eqn]$ (Fig 5). This is likely due to reductions in fishing mortality at low abundance, but since the fitted values did not overlap, some fishing mortality remained even at very low abundance. For the vast majority of stocks, the lowest model-estimated values of $[eqn]$ (75 of 77) were found at the maximum abundance in the time series (Tables 2 and 3 and S1 and S2 Figs). However, two Atlantic cod stocks in the Northwest Atlantic had their lowest population growth rates observed at abundances markedly below the maximum. These Atlantic cod stocks seemed to have the least compensatory capacity at relative abundances near 40% of their maximum.

Table 3: A comparison of the difference in λdem and λreal for each stock at the minimum and maximum abundance, along with the change at 40% of the maximum abundance. Positive values indicate that the λ estimate has increased as abundance declines (consistent with negative density dependence), while negative values indicate the λ estimate has declined (consistent with positive density dependence).

Density dependence of λ.Panel (a) compares the GAM smooths for λdem and λreal and the GAM smooth and the GAM estimated 95% confidence interval. Panel (b) shows λdem against the proportion of maximum abundance. The points represent the estimate for each year and stock along with the GAM smooth and the GAM estimated 95% confidence interval. Panel (c) shows λreal against the proportion of maximum abundance. Panel (d) shows the data and GAM smooths for two Atlantic cod (G. morhua) stocks in the Northeast Atlantic and two stocks in the Northwest Atlantic.

There was generally an increase in the variability of the $[eqn]$ estimates as abundance declined, largely due to the increase in the magnitude of the maximum value of $[eqn]$ (Fig 5). While there was a consistent overall pattern, the estimated random effects demonstrated substantial variability among stocks via the strength of negative density dependence. For example, stocks in the Western Atlantic (particularly Atlantic cod) and Alaska had some of the lowest $[eqn]$ estimates at low abundance, while several stocks in the North Sea (including Atlantic cod) had among the highest $[eqn]$ estimates (Table 2). Of the 66 stocks that declined below 40% of the maximum abundance, 7 had $[eqn]$ estimates < 1 at 40% of their maximum abundance; when the effects of fishing were included, 42 of the stocks had $[eqn]$ below 1 (Table 2).

Recovery potential

When excluding fishing mortality, the overall median lifetime reproductive success (i.e., the number of female spawners per female spawner) was 2.8 (mean = 11.1; Figs 6 and S3 Fig). When fishing mortality was accounted for, the median reproductive success dropped to 1 (mean = 4.3). This difference between the mean and median for stocks experiencing exploitation again highlights the importance of occasional large reproductive events in sustaining these harvested stocks (Figs 6 and S3 Fig). For species found across the North Atlantic, the lifetime reproductive success when excluding fishing mortality tended to be higher in the east than the west, a pattern most evident in Atlantic cod (G. morhua) stocks. While the reproductive potential of cohorts within a stock could be highly variable, there were no clear temporal trends in lifetime reproductive success (S4 and S5 Figs).

𝐑0 Estimates for each stock of G. morhua, M. aeglefinus, P. virens, C. harengus, P. platessa, and S. solea.The blue values are calculated using natural mortality and grey includes both natural mortality and fishing mortality. The red stars indicate the median and the dashed vertical line represents a lifetime reproductive output of one. The x axis is on the log scale.

Generation length (calculated as the mean age of spawners contributing to a recruit cohort) was lower and tended to be less variable under fishing across all stocks (Figs 7 and S6 Fig). Similar to other metrics, there were regional and inter-specific patterns. For example, Atlantic herring (C. harengus) is typically considered to be a short-lived species, yet its generation length in the Central and Western Atlantic tended to be higher than the generation length of Atlantic cod (G. morhua) and haddock (M. aeglefinus) in the Eastern Atlantic (Fig 7). Similar to lifetime reproductive output, there were no obvious temporal trends in the generation length among cohorts for the stocks in this analysis (S7 and S8 Figs).

Generation time estimates for each stock of G. morhua, M. aeglefinus, P. virens, C. harengus, P. platessa, and S. solea.The blue values are calculated using natural mortality and grey includes both natural mortality and fishing mortality. The red stars indicate the median.

Based on simulated stock trajectories, the doubling time in the absence of fishing was relatively high for the majority of species (Fig 8). Of the 77 stocks in this analysis, 39 were classified with a high recovery potential, which increased to 68 if those having a moderate recovery potential after 20 years were included (S9 Fig). Unlike stocks with higher recovery potential, there were clear regional patterns among stocks with low recovery potential. Of the 9 stocks in this analysis classified as having low recovery potential, 6 were in Alaska and 3 were in NAFO Regions 2–4 of the Western Atlantic. This meant that 55% and 43% of assessed stocks in Alaska and NAFO Regions 2–4, respectively, were classified as having a low recovery potential. This contrasts to 0 of the 44 stocks in the Eastern Atlantic classified with low recovery potential in the absence of fishing (Fig 8).

Proportion of the simulations in which a given stock has doubled in size when there is no fishing mortality.Stocks are grouped by taxonomic order, with 39 Gadiform stocks, 2 Scombriform stocks, 1 Scorpaeniform stock, 20 Pleuronectiform stocks, 11 Clupeiform stocks, and 4 Perciform stocks. The x-axis represents the year of the simulation, with the first 20 years shown.

Doubling time was influenced by both the magnitude and variability of $[eqn]$ . Stocks with a median $[eqn]$ were classified as having moderate to high recovery potential whenever the mean $[eqn]$ was $[eqn]$ 1.06 or the standard deviation in $[eqn]$ was $[eqn]$ 0.29 (Fig 9). In contrast, stocks with a median $[eqn]$ near 1 and with low variability tended to be classified as having low recovery potential.

The doubling time recovery potential, the colors represent low (grey), moderate (white), or high(blue) for each stock after 20 years compared to the medianλdem and the standard deviation of λdem. The red dashed line on the left panel is a λdem value of one.

Discussion

By using life-tables to assess the recovery potential of 77 commercially harvested marine fish stocks, we found that most stocks have an innate capacity for population growth in the absence of fishing, with 88% displaying a moderate to high recovery potential. However, we found that these stocks declined in size in the majority of years (i.e., $[eqn]$ was < 1 in 56.8% of the years) which suggests that stocks were managed in such a way that large but infrequent recruitment events were required to sustain productivity over the long term. Quota-based fisheries management strategies often implicitly or explicitly aim for neutral risk (i.e., a 50% probability) of population decline when setting harvest advice (e.g., [56]; [57]), however, this analysis indicates that this was not achieved in practice. In addition, there was considerable variability in the observed demographics and recovery potential within and among species and stocks. For example, while some Atlantic cod stocks exhibited strong negative density dependence and potential for recovery, other stocks demonstrated a limited capacity for population growth. Collectively, these findings highlight the importance of developing context-specific recovery strategies that consider the demographics of individual stocks, as a reliance on only life-history characteristics will likely prove insufficient to predict recovery potential.

Moreover, when the effects of fishing were excluded, these results indicated these stocks retained an innate capacity for increase of approximately 8.7%. If fishing ceased, the stocks were likely to recover (i.e., double in size) within 20 years whenever the median $[eqn]$ was $[eqn]$ 1.06 or the standard deviation of lambda was $[eqn]$ 0.29. Widespread distributions, high egg production per female, and broadcast spawning strategies suggest that harvested marine fishes should possess very high potential for annual population growth compared to other vertebrates [21]. Therefore, when environmental variation results in favorable conditions, stocks are anticipated to exhibit high survival, resulting in large recruitment events or higher overall survival (e.g., [58], [59]). While we found evidence of this, stocks with high variability in $[eqn]$ also had high recovery potential, however, both intra- and inter-species variability did not conform to expectations from conventional life-history theory [25,60,61]. For example, Atlantic herring are a small, relatively short-lived, low-trophic-level pelagic species; their life-history characteristics are considered to result in highly variable population growth potential, which should allow them to exploit favorable environmental conditions and increase quickly in abundance [62,63]. However, the variability of $[eqn]$ for several Atlantic herring stocks was lower than for the higher trophic level and longer-lived Atlantic cod and haddock stocks in the same region. Furthermore, these three species co-occur throughout the North Atlantic and the timespan of data collection was similar, suggesting they likely experienced similar variability in oceanographic conditions. Other stocks of small, lower-trophic-level species, such as lesser sand eel (A. tobianus) and Norway pout (T. esmarkii), were characterized by higher $[eqn]$ values, which aligned with theoretical expectations for these species. However, these two stocks were also found in the North Sea, which appeared to be an environment in which high variability in $[eqn]$ is commonplace for harvested stocks.

These results suggest that equilibrium theory and general life-history characteristics (e.g., longevity, trophic level) are poor predictors of population growth rates and recovery potential. For stocks that have collapsed or are considered to be at very low abundance, it would be prudent to evaluate the observed demographic parameters for the stock prior to developing harvest advice. Our analysis demonstrates that several stocks had little capacity to increase in abundance even in the absence of fishing, which would guarantee further declines if any level of harvesting continued. The stocks with limited recovery potential belonged to relatively distinct taxonomic and regional categories. For example, 7 stocks had $[eqn]$ estimates below one at 40% of their maximum abundance, indicating these stocks had limited capacity for population growth at low abundance in the absence of fishing. Six of these stocks were in the Gadidae family, with four in the NW Atlantic (three of which were Atlantic cod stocks) and two Alaskan pollock stocks (G. chalcogrammus). These stocks exhibited characteristics similar to endangered populations, which can also have lifetime reproductive rates < 1 [64,65].

In our analysis, the most marked example of intra-species variability in productivity was Atlantic cod across the North Atlantic. We have little understanding of how large most of these Atlantic cod stocks would have been before being impacted by large-scale harvesting, but historical reconstructions suggest the stocks in Iceland and the Western Atlantic were substantially larger [66–68]. Similar to many wild salmon populations [69,70], the combination of potential extreme abundance declines and enviornental changes may have resulted in these cod stocks having little innate capacity for growth or recovery. Indeed, it has been suggested that the Southern Gulf of St. Lawrence cod stock will go extinct unless natural mortality declines significantly [71,72], a conclusion supported by our analysis. However, the low recovery potential of Atlantic cod was not observed in the Eastern Atlantic; the stocks around mainland Europe all had high recovery potential. These regional differences may be related to fisheries-induced evolution, which has been linked to substantial change in the life-history characteristics (e.g., age and size at maturity), the loss of large old spawners, and/or to environmental conditions within the Western Atlantic [73–78]. This indicates that regional variability must factor strongly into a species’ ability to withstand fishing pressure, irrespective of life-history. Given that variability in growth potential among populations could be substantially greater than among species (e.g., compare the variability of $[eqn]$ between Atlantic cod and Atlantic herring), life-history characteristics and/or trophic level appear to be poor predictors of how marine fishes will respond to harvesting. This problem would compound for data-limited stocks, wherein these assessments often rely on species life-history traits (e.g., [79]). Consequently, such assumptions are more likely to lead to inappropriate harvest advice. That life-history characteristics did not link strongly to expected productivity is extremely concerning, given the need for sustainable fisheries management of data-limited stocks to mitigate risks of overfishing across the globe [79,80].

While no clear temporal trend was observed in $[eqn]$ , both the 2000s and 2010s were unusual relative to the previous three decades. In the decades with sufficient information to compare, the $[eqn]$ estimates were lowest in the 2000s and only the 1970s had a lower $[eqn]$ when accounting for the effects of fishing. The median $[eqn]$ remained relatively low in the 2010s, however the variability of the $[eqn]$ estimates was higher than what had been observed previously, resulting in a mean $[eqn]$ similar to what was observed in earlier decades. Intriguingly, the abundance continued to decline relative to its maximum and reached the lowest level observed in the 2010s. It was approaching 40% of the maximum abundance, a level in which compensatory dynamics begin to be observed more clearly for many species [50]. Compensatory response should cause $[eqn]$ to increase. Thus, the increase in the variability in the 2010s may be related to stronger compensatory dynamics in some stocks (e.g., an increase in years with relatively large recruitment events). However, it is also possible that broad scale climate effects could increase variability, given that phenological changes (e.g., timing of reproduction, physiology, behaviour) brought on by environmental variability are expected to affect demographic rates and thus population dynamics [81]. Finally, estimates of cohort size from age structured models are the most uncertain at the end of the time series, because the most recent cohorts have not yet fully moved through the model. The apparent increase in variability of the $[eqn]$ in the 2010s could simply be related to additional variability in the estimated size of recent cohorts.

Using historical datasets intrinsically constrained the population growth rates estimated from the life-table analyses. For example, when the stock was at its highest abundance in the observed time series, $[eqn]$ , by definition, had to be $[eqn]$ 1. Conversely, at the lowest abundance observed $[eqn]$ had to be $[eqn]$ 1. There were also more subtle constraints on the $[eqn]$ estimates. For example, $[eqn]$ had to be $[eqn]$ 2 when a stock was at 50% of the maximum abundance, and the minimum value of $[eqn]$ was similarly constrained by the minimum observed abundance observed in the time series. These constraints all contributed to an increase in the variability of $[eqn]$ as abundance declined. However, they affected all stocks in the same manner and thus would not introduce systematic bias into our understanding of regional or species-specific patterns. In addition, we assumed the fishery removals in the stock assessments were correct and modified the productivity parameters to match the abundance time series from those stock assessment. This resulted in $[eqn]$ values that were 6.7% lower than those calculated using the initial productivity parameters found in the stock assessments. If this decline in $[eqn]$ was instead attributed to fishing mortality (i.e., unobserved removals), then recovery potential in the absence of fishing would have increased (e.g., the doubling time would have declined).

Our life-table reconstruction enabled the underlying demographic rates to be estimated for stocks with sufficient age-structured data and provides an opportunity to compare the dynamics of these stocks using concepts developed in conservation biology and ecology. In doing so, our optimization procedure assumed that both the abundance estimates and the fishery removals from the stock assessment were correct. Although these assumptions are unlikely to be fully met, harvest advice is being derived for these stocks under the same general assumptions. Thus, our analysis reveals that current assessments can incorporate vastly different understandings of productivity and the potential for population growth for the same species across different regions, even when incorporating similar elements to describe population dynamics (e.g., equilibrium theory). Stocks characterized by low variability in population growth rates do not appear to be overly influenced by variation in environmental conditions (or that there was little environmental variability in the ecosystem), while those with a high innate capacity for population increase could indicate the opposite. However, it is not immediately apparent which conditions (e.g., environmental, ecological and/or anthropogenic) might be the primary determinants of population growth rates for a given stock, species, or system. Our findings underscore the need for alternative frameworks for deriving advice, particularly those that incorporate consideration of inter-species relationships and/or environmental conditions under the umbrella of ecosystem-based fishery management [82,83]. They also underscore the need to collect stock-specific data on marine fishes, rather than relying on species-level information from global repositories (e.g., Fishbase, Fishlife; [84]) or assuming values from closely related species to derive harvest advice (e.g., in cases in which there is insufficient data to estimate stock productivity, 79,[85]). While our analysis suggests that management advice may be fairly accurate for stocks with moderate variability in population growth rates or those that retain high innate capacity for population increase, it is likely to be the least accurate when recovery potential is low.

Supporting information

S1 FigDensity dependence of $[eqn]$ for each stock in the analysis plotted regionally.(TIFF)

S2 FigDensity dependence of $[eqn]$ for each stock in the analysis plotted regionally.(TIFF)

S3 Fig $[eqn]$ Estimates for the remainder of the stocks in this analysis (x-axis), plotted by the location of the stock(panels).The abbreviated species name is provided on the y-axis.(TIFF)

S4 Fig $[eqn]$ Estimates for every cohort for each species excluding the effects of fishing.(TIFF)

S5 Fig $[eqn]$ Estimates for every cohort for each species including the effects of fishing.(TIFF)

S6 FigGeneration length estimates for the remainder of the stocks in this analysis (x-axis), plotted by the location of the stock (panels).The abbreviated species name is provided on the y-axis.(TIFF)

S7 FigGeneration length estimates for every cohort for each species excluding the effects of fishing.(TIFF)

S8 FigGeneration length estimates for every cohort for each species including the effects of fishing.(TIFF)

S9 FigThe number of stocks (-axis) which doubled in size after a given number of simulated years (x-axis) with no fishing mortality The dashed line represents the 77 stocks in the analysis.(TIFF)

S10 FigBoxplots of $[eqn]$ estimates from the analysis in which the abundances were back-calculated for stocks in which the youngest age-class data available were greater than 1.A subest of the results are shown here, with data for each stock of G. morhua, M. aeglefinus, P. virens, C. harengus, P. platessa, and S. solea. The blue stars indicate the mean and the dashed horizontal line represents a $[eqn]$ of one. The y axis is on the log scale.(TIFF)

S11 FigBoxplots of $[eqn]$ estimates from the analysis in which the abundances were back-calculated for stocks in which the youngest age-class data available were greater than 1.A subest of the results are shown here, with data for each stock of G. morhua, M. aeglefinus, P. virens, C. harengus, P. platessa, and S. solea. The blue stars indicate the mean and the dashed horizontal line represents a $[eqn]$ of one. The y axis is on the log scale.(TIFF)