Summing the Parts: Improving Population Estimates Using a State-space Multispecies Production Model

Conference article

J. Northw. Atl. Fish. Sci., Vol. 57, Special Issue: 1–15

Paul M. Regular^*, Mariano Koen-Alonso, M. Joanne Morgan, Pierre Pepin^†, Rick M. Rideout

Fisheries and Oceans Canada, Northwest Atlantic Fisheries Centre, 80 East White Hills, St. John’s, Newfoundland and Labrador, A1C 5X1, Canada
*Corresponding author; E-mail: Paul.Regular@dfo-mpo.gc.ca
†Current address; Three Dog House, 1023 Indian Meal Line, Portugal Cove - St. Philip’s,
Newfoundland and Labrador, A1M 3C4, Canada

 

Regular, P.M., Koen-Alonso, M., Morgan, M.J., Pepin, P., and Rideout, R.M. 2026. Summing the parts: Improving population estimates using a state-space multispecies production model. J. Northw. Atl. Fish. Sci., 57(2): 1–15. https://doi.org/10.2960/J.v57.m750

Abstract

Carrying capacity is a fundamental concept in ecology that has inspired the development and application of a broad range of population models. In the context of fisheries science, production models have been employed globally to calculate carrying capacity and guide the sustainable use of fish populations. Production models have, however, been criticized for failing to account for species interactions and environmental effects. We aim to fill some of these gaps by introducing a novel state-space multispecies production model. We apply our extended model to commercially important demersal fish species off the east coast of Canada to assess its ability to reveal species interactions and the relative impacts of fishing and environmental effects. Our results indicate that accounting for species interactions increases the accuracy of biomass estimates for species within a community. The model also revealed strongly correlated process deviations, unrelated to fishing or density-dependent effects, which unexpectedly indicates that widespread collapses were primarily driven by a common environmental driver rather than fishing. Such inferences indicate that this may be a promising avenue for producing more holistic and accurate assessments for multiple species with relatively minimal data requirements (time-series of landings and fisheries-independent indices). Finally, this approach may serve as a stepping stone towards an ecosystem-based approach to fisheries management.

 

Keywords:Multispecies production model, State-space framework, Ecosystem-based fisheries management, Demersal fish populations, Environmental drivers

PDF, Supplemental Materials

Introduction

The concept of carrying capacity has long been foundational in applied population ecology, being widely used in the management of renewable resources (Chapman and Byron, 2018; Hilborn et al., 1995). The understanding that populations produce more offspring than an environment can sustain led to the notion that the ‘surplus’ can be harvested sustainably (Pauly and Froese, 2021). These ideas are exemplified by the concept of maximum sustainable yield (MSY) in the classic single-species surplus production modeling framework (Schaefer, 1954). In one equation, this framework attempts to explain interannual changes in biomass using fisheries landings and estimates of intrinsic growth rate and carrying capacity. Given the theoretical elegance of the approach, it has been both widely adopted and scrutinized. Estimates of carrying capacity, and resultant derivations of MSY, are frequently criticized for being time-invariant, which ignores the ubiquity of environmental variation (Del Monte-Luna et al., 2004). Moreover, traditional surplus production models tend to focus on single-species dynamics and often disregard the complexities arising from species interactions within ecosystems (Gamble and Link, 2009).

Recognizing the limitations of single-species approaches, there have been calls to move towards an ecosystem-based approach to fisheries management (EBFM; Latour et al., 2003). EBFM acknowledges the intricate web of ecological interactions and aims to ensure the sustainability and integrity of marine ecosystems while supporting viable fisheries (Pikitch et al., 2004). To successfully implement EBFM, it is crucial to develop models that account for species interactions and the dynamics of multiple species within the ecosystem. Substantial progress has been made in the development of multispecies models and a spectrum of approaches have been developed, ranging from complex models that attempt to account for all parts of marine ecosystems (e.g., Fulton et al., 2011) to multispecies age-structured assessment models (e.g., Albertsen et al., 2018) to multispecies surplus production models (e.g., Bundy et al., 2012; Gamble and Link, 2009; Mueter and Megrey, 2006). However, application of these approaches to fisheries management have often been hindered by data limitations (e.g., age data are frequently not available) and knowledge gaps (e.g., incomplete understanding of food-web interactions). There is therefore a need for methods to help bridge the gap between single-species and multispecies assessment in data or information poor circumstances.

In this paper, we borrow concepts from single-species surplus production modeling (Millar and Meyer, 2000) and multispecies modeling (Albertsen et al., 2018) to construct a model that incorporates the impacts of fishing on single-species populations and accounts for species interactions within the ecosystem. The data requirements of this model are relatively minimal, requiring species-specific landings and survey indices of biomass. As a case study, we apply this model to commercially important demersal fish species off the east coast of Canada, specifically the Grand Banks of Newfoundland. This case study is particularly germane due to the widespread collapse of most stocks in the area during the early 1990s (Lear, 1998), leaving the relative contributions of fishing and environmental impacts uncertain (Pedersen et al., 2017). Using this case study, we aim to reveal species interactions, distinguish the impacts of fishing from environmental effects, and demonstrate the utility of this modeling approach. The subsequent sections of this paper will present the conceptual framework of our model, describe its application, and discuss the implications of our findings for ecosystem-based fisheries management.

Methods

Model formulation

Trends in fish populations have frequently been described using state-space production models of the form

where By is biomass at the start of year y, Cy is the catch through year y, gB is production as a function of biomass, Iy,s is an index of relative abundance in year y from survey i, qi is the time-invariant catchability coefficient for survey index i, δ is process error, and ε is observation error. Statistical challenges aside, the most difficult aspect of this model to parameterize is the production function as it needs to capture changes caused by growth, recruitment, and natural mortality. Schaefer (1954) proposed a solution by applying the logistic equation to describe self-limiting growth,

where r is the maximum per-capita rate of change and K is the carrying capacity. That is, a populations’ intrinsic ability to grow (rB) is limited by the size of the current population relative to the maximum biomass the system can support (1-B/K). While this formulation offers an elegant description of single-species population dynamics, it assumes that density-dependent effects are solely caused by intraspecific competition and ignores the potential effects of other species inhabiting the same ecological area, competing for the same resources. We present an extension of equation (3) that models the growth of a focal species s while accounting for intra- and interspecific competition. Here, density-dependent effects are incurred when the total biomass of all species, indexed by s', exceeds the capacity of the system, K_Σ,

While intrinsic rates of growth may vary across species, this formulation implies that the growth of all species is ultimately limited by the finite amount of energy in a region (i.e., as the total population of all species in the system increases towards KΣ, year-over-year growth of all species slows). Combining equations (1), (2), and (4), our model becomes

One simplifying assumption of this formulation is that each unit of biomass, regardless of species, contributes equally to limiting the realized growth of any given species. This assumption is partially relaxed through the covariance structure, which allows positive or negative associations among species beyond the shared carrying capacity. Though such associations may arise from the observation process, we assume that most covariance stems from ecosystem processes. We therefore treat observation errors as independent and normally distributed deviations such that εy,i,s∼N0,τi,s, where the standard deviation parameter τi,s represents species- and survey-specific levels of observation error.

To account for temporal dependencies driven by ecological processes, a more flexible error structure is used to describe the process errors. Species interactions may drive positive or negative population responses resulting from direct or indirect associations. Deviations from expected production may also display temporal dependence if the factors contributing to the process errors change gradually over time. Such inertia may cause positive or negative process errors to persist over several years. A first-order autoregressive (AR1) process was therefore applied to account for temporal dependence. Both sources of dependence are modeled using a multivariate AR1 process where

 

 

Here δy=δy,1,…,δy,S⊤ is the vector of process deviations across species. The parameter ϕ controls temporal autocorrelation, where values near 0 indicate weak persistence and values near 1 indicate strong persistence. Although different species may exhibit different levels of temporal correlation, a single ϕ is assumed to avoid over-parameterization. Nevertheless, effective differences in temporal correlation among species can still emerge through the covariance matrix Σ, which captures cross-species dependence:

 

 

Species-to-species correlations are described by ρs,s' (s≠s'), with values between -1 and 1 representing negative to positive association. This flexible structure allows the testing of alternative hypotheses that process errors are independent through time or across species (ϕ=0 or ρs,s'=0). The possibility that process errors are similarly correlated across all species may also be tested by estimating only one ρ parameter. Finally, the magnitude of the process error deviations is controlled by the species-specific standard deviation parameters, σs.

In implementation, Σ is parameterized as D R D, where D=diagσ1,…,σS and R=LL⊤ is a correlation matrix obtained from a Cholesky factor L constructed from unconstrained parameters. This formulation ensures that Σ is positive definite.

Minor extensions of the formulation also permit the fitting of covariates which may describe an underlying linear effect. Two options were implemented, one that affects the process errors by substituting eδy,s in equation (5) with eβδXy,sδ+δy,s, and another that affects the carrying capacity by substituting KΣ in equation (5) with KΣeβKXyK, where β parameters capture linear effects of covariates included in design matrices, X. The idea is that some factors may affect positive or negative changes in the populations while others may affect change in the total carrying capacity of the system. A covariate option for intrinsic rates of increase was not implemented since one goal of this model is to obtain estimates of this vital rate, which is not expected to change rapidly as it is shaped by natural selection (Hutchings, 2021). Equation (5) was also modified to enable the fitting of a single-species Schaefer production function by dropping the summation of biomass and estimating species-specific carrying capacities, Ks, rather than a system-level carrying capacity, KΣ (i.e., applying equation (3) separately to each species).

Statistical framework

This model was implemented using template model builder (TMB; Kristensen et al., 2016), which is a R (R Core Team, 2021) package that enables the fitting of complex nonlinear random effects such as the latent B variable in state-space production models (equation (1)). Such variables are not directly measured but are inferred indirectly via observed values. Data fitting is accomplished using a combination of Laplace approximation and automatic differentiation to evaluate the joint likelihood (Kristensen et al., 2016). Like the production model described by Pedersen et al. (2017), both frequentist and Bayesian inference of model parameters are possible. In development, we found that estimation of parameters was generally more successful when vaguely informative priors were specified, and in some cases, not estimable when unconstrained.

Case study

Data

The multispecies production model described above requires two basic inputs for each species: a time-series of catch (Cy,s in equation (5)), and an index of population size (Iy,s in equation (6)). The Northwest Atlantic Fisheries Organization (NAFO) and Fisheries and Oceans Canada (DFO) have been collecting and curating such information for multiple fish populations along the shelves of Newfoundland and Labrador (NL) since the 1970s. The communities inhabiting these shelves can be divided into several regions with distinct productivity i.e. ecosystem production units; Pepin et al. (2014). For our case study, we tallied catch data and calculated survey indices of multiple demersal fish populations from three regions (Fig. 1): 1) the Northeast NL Shelf (NAFO divisions 2J3K), 2) the Grand Bank (NAFO divisions 3LNO), and 3) Southern NL (NAFO sub-division 3Ps).

Catch data were extracted from NAFO’s STATLANT 21A database (https://www.nafo.int/Data/STATLANT-21A, accessed 2022–01–21) and aggregated by region, species, and year. Survey indices were derived from the standardized, stratified random bottom-trawl surveys conducted each spring and fall by DFO; one of the largest fisheries-independent surveys in the world, covering more than 500 000 km² annually (roughly the size of Sweden or the Yukon, Canada) to depths up to 1 500 m. Since the inception of this program in 1971, survey platforms and protocols have undergone a series of changes that affect the continuity of the data collected in each region and season. A Yankee then Engel otter trawl, with nets designed to catch large demersal fish, were used from 1971 to 1994. Starting in the fall of 1995 survey gear was changed to a Campelen shrimp trawl with a small mesh codend, which allowed a broader range of species and size groups to be captured (Chadwick et al., 2007). Within each era of the survey (Yankee, Engel, or Campelen) and for each season and region, samples used in this study were limited to strata that were sampled most years (> 80%) and to species found across more than 10% of these core strata. This often resulted in the exclusion of strata >750 m because these areas have been inconsistently sampled by the survey. Stratified analyses (Smith and Somerton, 1981) were then conducted on the remaining species to obtain indices of total biomass. To minimize bias introduced by inconsistent survey coverage, indices from years where more than 20% of the biomass was likely missed, inferred from time averaged percent occupancy within strata, were excluded from our analysis [sensu NAFO guidelines; page 10, NAFO (2019)]. Finally, species were ranked by cumulative commercial catch and limited to the seven most commonly caught species, or species group when catch was not consistently distinguished by species, within each region. On the Northeast NL Shelf, the included species were redfish spp. (Sebastes fasciatus and Sebastes mentella combined), wolffish spp. (Anarhichas lupus and Anarhichas minor combined), witch flounder (Glyptocephalus cynoglossus), American plaice (Hippoglossoides platessoides), Greenland halibut (Reinhardtius hippoglossoides), Atlantic cod (Gadus morhua), and skate spp. (Amblyraja radiata and Malacoraja senta combined). On the Grand Bank, the included species were redfish spp., yellowtail flounder (Myzopsetta ferruginea), American plaice, Greenland halibut, haddock (Melanogrammus aeglefinus), and Atlantic cod. Finally, along Southern NL, the species included were redfish spp., witch flounder, American plaice, white hake (Urophycis tenuis), haddock, Atlantic cod, and skate spp.

Priors

For simplicity, all priors were normally distributed and, in most cases, upper and lower inflection points (μ±σ; ~68% of the total area under the curve) of each normal prior were defined using values in log or logit space. Upper and lower values were based on previous research or knowledge to impose fairly generic and vaguely informative priors. Interactive dashboards accompanying this paper provide visualizations of the priors and posteriors for core parameters (Supplement 1); here we summarize their parameterization and report the ranges used to define each prior.

Intrinsic growth rate, rs, and carrying capacity, K

To capture a broad range of intrinsic growth rates (Thorson, 2020), 0.1 and 1 were chosen as the lower, rl, and upper, ru, values; this translates to a normal prior with a mean of -1.15 and a standard deviation of 1.15 on the log-scale. The log-scale KΣ prior was informed by total levels of catch and rl and ru,

where KΣ,l and KΣ,u are the lower and upper values. On the lower end, this prior imposes the assumption that the fishery is unlikely to have caught the equivalent of KΣ and, on the upper end, it assumes that the maximum observed catch represents a portion of KΣ (sensu Froese et al., 2017). The division by the lower and upper values for r also accounts for the potential range of productivity. Note that the maximum time-series catches are species specific when estimating species-specific carrying capacities.

Starting biomass, B1,s

Like KΣ, catch was used to constrain the plausible range of the biomass at the beginning of the time series, B1,s, which we will denote as B0s to simplify notation. Specifically,

where B0l,s and B0u,s are the lower and upper values that are log transformed to define the normal prior for logB0s. Again, adjusting for the potential range of productivity from rl to ru, this assumes that the fishery did not catch all of the biomass in the first year and it assumes that the catch represents a portion of the biomass in the first year. A potential flaw with the upper value for this prior is that a lack of market demand may contradict the assumption that landings are coarsely proportional to stock size. However, the upper range chosen was considered reasonable for this case study as there was an active fishery for each species examined here at the start of the time series.

Process error variance, σs2, temporal correlation, ϕ, and species-to-species correlation, ρs,s

Considerable process variability may arise from variable recruitment, natural mortality, and/or growth. For instance, species in the Scorpaenidae family, which includes the Sebastes genus, have notoriously variable recruitment which frequently results in “spasmodic” stock dynamics (Cadigan et al., 2022; Licandeo et al., 2020). Also, there is evidence that heightened and variable natural mortality contributed to the collapse and slow recovery of Atlantic cod along the Northeast NL Shelf (Cadigan, 2015; Regular et al., 2022). To account for a wide range of possible standard deviations, σs, a vague prior with 0.01 and 1 as the lower and upper values were chosen. In log-space, this translates to a normal distribution with a mean of -2.3 and sd of 2.3.

There is little knowledge to determine the potential level of temporal or species-to-species dependence present in the focal systems; vague priors were therefore defined for ϕ and ρs,s. For ϕ, 0.1 and 0.9 were logit transformed [logx/1-x] to define the upper and lower inflection points for a normal prior for the logit of ϕ. For ρ, the logit transformation was shifted [log1+x/1-x] to capture negative and positive lower (-0.9) and upper (0.9) values. When these parameters are estimated, these priors give most credence to moderate temporal correlation (0.5) and no species-to-species correlation (0) but still allow the possibility of high levels of correlation (0.9).

Catchability, qi,s, and observation error variance, τi,s2

The catchability of the spring and fall surveys conducted by DFO likely changed over time given the shift in gear from a Yankee to Engel to Campelen trawl as well as spatial shifts in survey coverage of the strata in each region. Survey catchability qi,s is therefore indexed by season and gear, i, and species, s. Moreover, the lower inflection points for the qi,s priors were informed by the average survey coverage by gear and season. Survey coverage within each region was computed by dividing the average spatial coverage of strata across years by the total area of all strata (i.e., average area covered / area of survey domain). To account for the fact that the proportion of a stock accessible to the survey may be lower than the nominal coverage, this ratio was further adjusted: multiplied by 0.2 for deep-water species (Greenland halibut, Atlantic halibut, witch flounder, redfish spp., white hake, silver hake, and monkfish) or by 0.5 for shallower shelf-associated species. These multipliers reflect prior expectations that gear selectivity (e.g., escapement under the footgear; Walsh, 1992) reduces availability for all species, while deep-water species require an additional adjustment because a portion of the stock may occupy depths beyond the survey domain. The upper inflection point was set to 1 as it is possible that the survey indices represent overestimates of the true population size in some instances.

A prior for observation error variance was informed by unbiased design-based estimates of survey variance associated with annual biomass estimates (Smith and Somerton, 1981). These estimates were used to calculate the coefficient of variation (CV) for each survey, year, and species; specifically,

where I represents the design-based indices of biomass used in equation (2). These CVs were log transformed and survey, i, and species, s, specific means and standard deviations were used to inform the prior for τi,s2 in log-space. Rather than treat the estimates of survey CV as a perfect indicator of total observation error, the prior was widened by multiplying the standard deviation of log CV by two to account for observation variance introduced by potential distributional shifts outside of the survey domain. Another two-time multiplier was applied for deep-water species because they have more scope for shifting outside the survey domain.

Model selection

The degree to which density-dependence of the demersal fish community inhabiting the NL shelves is driven by intra versus interspecific competition is unknown. Nor is the degree to which species interactions affect population dynamics. It is, however, well known that the community has been exploited for more than 500 years and this history was punctuated by a collapse of several stocks, most notably cod, in the early 1990s (Lear, 1998). This collapse was not isolated to commercial stocks, as the biomass of several non-commercial species collapsed at the same time. This was followed by a reorganization of the community, which implies that the system experienced a regime shift (Pedersen et al., 2017). Given this context, and the structure of the multispecies surplus production model described above, a series of hypotheses were tested:

Full:For this hypothesis, it is posited that population dynamics within regions are governed by an overarching system-level carrying capacity that affects all focal species as the aggregate biomass approaches system limits (i.e., inter-specific density-dependent effects are assumed). For this hypothesis and all others, annual reported landings of each species are accounted for in the production equation (5). Residual variations not explained by intrinsic growth, density-dependent effects or landings are described by process errors that are 1) correlated across species and said correlations are assumed to be unstructured, meaning relationships can be of differing strengths — positive, neutral, or negative — for each species-to-species pair; and 2) assumed to be temporally correlated, following an AR1 structure. Finally, a shift covariate was applied to the carrying capacity parameter to enable the estimation of different system limits before and after the community-wide collapse to assess support for a regime shift in each systems’ capacity for the focal demersal species.

Just shift:Model structure is the same as the full model, except the process errors are assumed to be independent and identically distributed random variables (iid) to assess the hypothesis that temporal and species correlations are explained by the shift.

Just correlation:Model structure is the same as the full model, except population dynamics are assumed to be affected by a common and time-invariant carrying capacity. The shift covariate was not applied.

Shared correlation:Model structure is the same as the just correlation model, except species-by-species correlations in the process errors are assumed to be the same across all pairs. This structure implies that there is a common but unknown environmental variable affecting the population dynamics of all species. The shift covariate was not applied.

Just species correlation:Model structure is the same as the shared correlation model, except the process errors at each time step are assumed to be iid. This structure implies that there is a common but unknown environmental process affecting all species, but the process is noisy with no temporal dependence. The shift covariate was not applied.

Just temporal correlation:Model structure is the same as the shared correlation model, except species-by-species correlations are assumed to be iid. This structure implies that environmental processes affect each species differently, however, there may be carry-over effects from one year to the next. The shift covariate was not applied.

No correlation:Model structure is the same as the shared correlation model, except both temporal and species correlations are set to zero (i.e., process innovations are iid across time and species). That is, population dynamics are thought to be affected by a common time-invariant carrying capacity and residual variations not explained by intrinsic growth, inter-specific density-dependent effects or landings are noisy and independent across time and species. The shift covariate was not applied.

Single-species:Model structure is the same as the no correlation model, except population dynamics are assumed to be governed by species-specific carrying capacities (i.e., intra-specific density-dependent effects are assumed). This formulation is analogous to standard state-space Schaefer production models. The shift covariate was not applied.

The predictive ability of each of these models was tested using two cross-validation approaches: 1) leave-one-out cross-validation (LOO-CV), and 2) hindcast cross-validation (Hindcast-CV). LOO-CV is a form of exhaustive cross-validation where the model is repeatedly conditioned on a training set missing one observation (i.e., one-fold) until the number of model folds equal the number of observations in the data. The missing observations are predicted at each fold, permitting assessments of the models’ ability to predict the actual value that was left-out at each fold. The hindcast-CV approach is similar, but it focuses on the models’ ability to predict the future. Under this approach, the model is repeatedly conditioned on a training set missing observations from the terminal year such that each fold excludes an increasing number of years worth of data from the tail of the time series (Kell et al., 2016). We folded back 20 years and, for each fold, predicted survey indices were compared to the observed survey indices (e.g., observed indices from 2020 were compared with predicted survey indices for 2020 from the model conditioned on data from 1978–2019). For both approaches, we denote predicted survey indices as Ij^ and the left-out observations Ij, where j represents all unique combinations of years, species, and survey indices present in the left-out data. LOO-CV and Hindcast-CV prediction error scores for each model for each region were calculated by taking the mean squared error,

These scores were also averaged across methods (LOO-CV, Hindcast-CV) and region (Northeast NL Shelf, Grand Bank, Southern NL) to obtain an overall score of predictive ability of each model. A well-fitting model will result in predicted values that are close to the excluded values and, therefore, result in lower scores. Ultimately, these scores enable a model-comparison approach to hypothesis testing. These scores are comparable even for the single-species model as it compares species-specific data points to predictions.

Code and Data Availability

All code used for the analyses in this study is publicly available in a research compendium on GitHub: https://github.com/PaulRegular/multispic. The compendium is organized as an R package, enabling straightforward replication of the results presented here and allowing the same methods to be applied to datasets from other regions.

Results

Both cross-validation metrics (LOO-CV and Hindcast-CV scores) indicate that most multispecies production model formulations outperform a single-species production model when applied to seven species within three ecosystem production units (Northeast NL Shelf, Grand Bank, and Southern NL) off the east coast of Canada (Fig. 2). Focusing on overall scores, the performance of the single-species production model was similar to a multispecies formulation that assumes there is no correlation in the process errors across species or time. There tends to be more notable decreases in the scores as temporal and species correlations are introduced, indicating an improvement in the predictive ability of these models. The “shared correlation” and “just correlation” formulations, in particular, tended to receive the lowest scores, and dropping the species and temporal correlations in lieu of a shift covariate (“just shift” formulation) resulted in a deterioration of predictive ability. Scores were improved when temporal and species-to-species correlations were introduced along with the shift covariate (“full” formulation); however, the fit of the “full” model tended to be poorer than the “just correlation” formulation, which further indicates that the “shift” covariate degraded the predictive ability of the model. Subsequent plots focus on the best fitting formulation, “just correlation”, to demonstrate model predictions.

The multispecies production model with unstructured species-to-species correlation and AR1 temporal correlation (“just correlation”) offered an explanation of the trends in survey indices of focal species across three ecosystem production units with little signs of systematic bias (see residual plots included in model dashboards; Supplement 1). Predicted indices track observed values and, by estimating catchability parameters by species and survey, indices from temporally fragmented surveys are stitched together and their trends are used to inform a continuous underlying trend in biomass (Fig. 3). The earlier Yankee and Engel eras of the Canadian surveys tended to receive lower catchability estimates than the Campelen era survey; indices since 1996 therefore tend to be closer, in relative terms, to the underlying estimates of biomass from the model.

Isolating residual changes in biomass not explained by reported fisheries landings or the production function (equation (4)) highlights deviations that can be either additions (positive residuals) or subtractions (negative residuals) to total biomass (Fig. 4). Substantial additions occurred throughout the 1980s, followed by sharp subtractions in the early 1990s across all three ecosystem production units. On the Northeast NL Shelf, biomass exceeded the estimated carrying capacity throughout much of the 1980s, and residual subtractions exceeding 500 kt per year precede the abrupt declines of the early 1990s. Most species have since recovered on the Northeast NL Shelf, with total biomass approaching carrying capacity while community composition has remained largely unchanged. On the Grand Bank, total biomass only briefly exceeded carrying capacity in the 1980s and, following declines in the early 1990s, subsequent increases were dominated by redfish spp. Off Southern NL, biomass remained below carrying capacity until the late 2010s, with smaller declines in the early 1990s and later increases largely driven by redfish spp.

Further inspection of the process errors reveals common patterns across all focal species across three ecosystem production units (Fig. 5). Like the exponentiated and unstandardized process errors (“residual change”) shown in Fig. 4, the standardized process errors highlight substantive subtractions in the early 1990s, representing time-series lows for 19 out of 21 populations. The standardized values also reveal parallel and periodic increases and decreases within each region. The species-to-species correlations in the process errors estimated by the “just correlation” model indicate that 84% (53 of 63) of pairs were positively correlated. Also note that the “shared correlation” model estimates of a common species-to-species correlation parameter were 0.79 (95% CI: 0.70, 0.86), 0.68 (95% CI: 0.58, 0.76), and 0.81 (95% CI: 0.66, 0.90) for the Northeast NL Shelf, Grand Bank, and Southern NL ecosystem production units, respectively. This model received similar cross-validation scores as the “just correlation” model which, taken together, further supports the inference that the process errors are primarily positively correlated across species within each production unit.

Discussion

Rather than assuming that populations are primarily regulated by species-specific carrying capacities, our multispecies production model assumes that both intra and interspecific competition stunts growth as total biomass approaches the environment’s maximal load. This assumption is conceptually similar to aggregate production models (Bundy et al., 2012; Fogarty et al., 2012; Mueter and Megrey, 2006) as it is rooted in the idea that total production, and consequently system-level MSY, is limited by the amount of resources available in a given ecosystem. In contrast to aggregate production models, we also attempt to capture the dynamics of species within a community. Fisheries landings, competitive interactions, predation, and prey availability all affect species-level production (Lotka, 1925; Schaefer, 1954; Volterra, 1926). While our model explicitly accounts for landings, species interactions are implicitly accounted for by estimating species-to-species correlations (sensu Albertsen et al., 2018; see Gamble and Link, 2009 for a more explicit approach). Finally, by utilizing a state-space framework akin to single-species state-space production models (e.g., Millar and Meyer, 2000; Winker et al., 2018), we attempt to differentiate population processes from noise and bias from surveys of fish populations. The overall structure of the model allows species-specific dynamics to be captured while avoiding the assumption that the dynamics of each species is isolated and independent from other species sharing the same ecosystem production unit and potentially competing for the same resources.

Our case study focuses on the population dynamics of commercially important demersal fish stocks off the east coast of Canada. Stocks in the area collapsed in the early 1990s (Lear, 1998), and the relative contribution of fishing and environmental impacts has been highly debated (Pedersen et al., 2017). We attempt to disentangle the impacts of fishing from environmental effects using our multispecies production model and, in doing so, we provide empirical evidence that environmental factors played a non-negligible role in the changes observed in the region. First, we found general support for models with a system-level carrying capacity, consistent with the expectation that species within the same ecosystem production unit are constrained by a finite amount of available energy (Pepin et al., 2014). Second, we found evidence for synchronous changes in the demersal fish community, which implies that a common bottom-up driver is impacting the dynamics of these species (see also Bundy et al., 2012). This is supported by species-specific studies on capelin (Mallotus villosus) and Atlantic cod in the region which highlight the influence of bottom-up drivers (e.g., Buren et al., 2014b, 2014a; Koen-Alonso et al., 2021; Regular et al., 2022). Taken together, evidence is mounting that fishing was not the sole cause of the collapses observed in the early 1990s.

Our inference that environmental factors were a key driver of stock collapses was unexpected given the compelling narrative that fishing activity was the primary driver (e.g., Gomes et al., 1995; Hutchings, 1996). Since our model utilizes reported fisheries landings, a portion of these losses may be attributed to illegal fishing activity. However, it seems unlikely that the industry had the capacity to extract the amount needed to match the estimated losses. For instance, annual catches in the late 1980s across the Northeast NL Shelf and the Grand Bank totaled ~450 kt while residual losses estimated by the model in the early 1990s was ~1 000 kt. The fishing industry would have had to covertly double its efforts to explain the declines. It follows that the decline must at least in part, if not primarily, be due to an unknown environmental driver. This contention is not new (Atkinson, 1994; see, for example, Morgan et al., 2002; Pedersen et al., 2017), however, it remains contentious and perplexing as we lack specific causal explanations. While increasingly cold conditions through the 1980s and early 1990s undoubtedly affected the distribution of multiple species (Montevecchi and Myers, 1997; Robertson et al., 2021; Rose et al., 2000), it is not yet clear whether shifting temperatures was the primary driver of the collapse and, if it was, the causal pathway has yet to be determined.

Regardless of the environmental driver behind the 1990s collapse, it is possible that increasingly industrialized and intense fishing activity through the 1960s and 1970s reduced population diversity and, consequently, hampered the ability of the species within the community to buffer subsequent environmental changes (sensu the portfolio effect, Schindler et al., 2010). Yet, a recent study found no evidence of genetic diversity loss in heavily exploited species like Atlantic cod (Pinsky et al., 2021). Another hypothesis is that fishing activity bounded the safe operating space of the system, triggering an alternate stable state (Scheffer et al., 2015). While not a perfect test of chaotic dynamics, we did assess the possibility of a systematic shift in system-level carrying capacities and found little support for this hypothesis.

Although we found little support for a structural shift, clear changes in community structure were evident. These shifts may have emerged from the combined effects of ocean climate fluctuations (Cyr et al., 2025), interspecific competition, and shifting energy pathways. It is well known that the dominant forage species in the area shifted from capelin to shrimp (Dawe et al., 2012) and this change was detrimental for cod (Link and Sherwood, 2019; Mullowney and Rose, 2014; Regular et al., 2022) and perhaps other piscivorous species that rely on capelin. Shrimp are an important prey item for redfish species (Brown-Vuillemin et al., 2022), so it is possible that the increasing shrimp population helped support concurrent recruitment pulses of redfish. We admit that this conjecture is speculative and acknowledge that the scale and dominance of redfish appears unrealistic; however, we add it as a simple example of how bottom-up forces may be driving the observed changes in the community. The reality is obviously more complex and the observed restructuring of the communities may be akin to the “paradox of plankton” where the continuous interaction of ecological and environmental factors give rise to “oscillations and chaos, with a continuous wax and wane of species within the community” (Scheffer et al., 2003). In any case, our results suggest that the observed productivity changes in these ecosystems may be better explained by underlying environmental drivers than by a sudden structural shift to an alternate stable state.

Like all models, our multispecies production model is an imperfect abstraction of nature and while it may be useful in some contexts, it is important to consider its limitations when interpreting results. First, it is important to remember that there may be a spatial mismatch in the structure and function of the populations included in this study as some stock boundaries differ from the regions used in this study. For instance, Atlantic cod in NAFO divisions 2J3KL are considered a separate stock from cod in divisions 3NO (Templeman, 1962) and here we split 2J3K (Northeast NL Shelf) and 3LNO (Grand Bank) into distinct regions. Assuming that our results are comparable to previous results, it is peculiar that they indicate that total biomass in the Northeast NL Shelf production unit was above the carrying capacity of that region through the 1980s. This finding contradicts historic records that suggest populations such as Atlantic cod in 2J3KL were at substantially higher levels in the 1970s and earlier (Rose, 2004; Schijns et al., 2021), which implies that the carrying capacity should be higher than estimated by the model presented here. Still, it is possible that the 1970s represents a period of unusually high productivity, where the system may have exceeded the carrying capacity (see also Cyr et al., 2025). In our model, carrying capacity should be viewed as a long-term average, with environmental variability creating transient conditions that can raise or lower realized productivity.

Results from the Grand Bank and Southern NL also indicate that the demersal fish community is currently dominated by redfish and, consequently, the system’s biomass appears to be approaching its carrying capacity. Though redfish are currently rebounding in parts of eastern Canada (Cadigan et al., 2022), the implication that they dominate the demersal community seems unrealistic. This may be an artifact of low estimates of survey catchability or the model’s inability to properly account for year effects. Observation errors are assumed to be lognormally distributed; however, extreme catch events or black swan events in space can introduce ‘year effects’ that may be better captured by a distribution with heavier tails, such as the t-distribution (Anderson and Ward, 2019). It is also possible that treating all species’ biomass as equivalent in limiting growth has inflated the apparent role of redfish. This highlights the need to better evaluate how species differ in their contributions to density-dependence.

Despite these limitations, the approach holds practical value because it relies on inputs that are routinely collected in many regions (fishery‐independent surveys and landings), making it broadly applicable in data-limited contexts. In turn, the multispecies production framework provides management-relevant outputs for ecosystem-based fisheries management (EBFM): it estimates community-level constraints (shared carrying capacity) and cross-species coherence, which can be used to design reference points and harvest control rules that recognize system-level limits and trade-offs among species (sensu Moffitt et al., 2016; Mueter and Megrey, 2006). This type of application could complement the ecosystem reference point-based framework currently implemented in NAFO to inform on the risk of ecosystem overfishing (Koen-Alonso et al., 2022; NAFO, 2022a, 2022b, 2022c). When residual dynamics become more synchronous, managers can apply stronger precautionary buffers to account for weakened portfolio effects (Schindler et al., 2010). The model can also support projections conditioned on community dynamics, useful for testing rebuilding trade-offs. Because it is lightweight relative to end-to-end ecosystem models (e.g., Fulton et al., 2011) yet consistent with key principles of EBFM (Pikitch et al., 2004), this approach offers a tractable option to embed multispecies dynamics and environmental considerations into existing assessment and strategy-evaluation workflows.

Conclusion

Practitioners are becoming increasingly aware of the need to apply EBFM (Pikitch et al., 2004); however, progress has been impeded by a series of myths (Patrick and Link, 2015). One myth is that “EBFM can only be implemented in regions with copious data, and the corollary, doing EBFM requires models that are too complicated” (Patrick and Link, 2015). Our results support the arguments of Patrick and Link (2015) that data requirements need not be prohibitive, models need not be overly complex, and perfect knowledge of every process is not required to make meaningful progress. Using standard survey and catch data, the multispecies production model presented here offers a tractable means of estimating community-level constraints and multispecies trends, as well as producing species-specific projections conditioned on recent community dynamics. As such, it provides a framework for developing reference points and testing harvest control rules that recognize trade-offs among species. Echoing earlier arguments (Pikitch et al., 2004), EBFM can and should be advanced even where data are limited. With minimal inputs, this method may help translate ecosystem considerations into operational advice, offering pragmatic support for the implementation of EBFM.

Acknowledgements

We thank the many ships’ crew and research staff involved in leading the surveys and collecting the data used in this study. These surveys were supported by Fisheries and Oceans Canada (DFO). An earlier draft version of the abstract, introduction, and discussion sections were written with the assistance of ChatGPT (March 14, 2023 version).

References

Albertsen, C. M., Nielsen, A., and Thygesen, U. H. 2018. Connecting single-stock assessment models through correlated survival. ICES Journal of Marine Science, 75: 235–244. https://doi.org/10.1093/icesjms/fsx114

Anderson, S. C., and Ward, E. J. 2019. Black swans in space: modeling spatiotemporal processes with extremes. Ecology, 100: e02403. https://doi.org/10.1002/ecy.2403

Atkinson, D. B. 1994. Some observations on the biomass and abundance of fish captured during stratified-random bottom trawl surveys in NAFO Divisions 2J and 3KL, autumn 1981–1991. NAFO Sci. Coun. Stud., 21: 1–12.

Brown-Vuillemin, S., Chabot, D., Nozères, C., Tremblay, R., Sirois, P., and Robert, D. 2022. Diet composition of redfish (Sebastes sp.) during periods of population collapse and massive resurgence in the Gulf of St. Lawrence. Frontiers in Marine Science, 9. https://doi.org/10.3389/fmars.2022.963039

Bundy, A., Bohaboy, E. C., Hjermann, D. O., Mueter, F. J., Fu, C., and Link, J. S. 2012. Common patterns, common drivers: comparative analysis of aggregate surplus production across ecosystems. Marine Ecology Progress Series, 459: 203–218. https://doi.org/10.3354/meps09787

Buren, A. D., Koen-Alonso, M., Pepin, P., Mowbray, F., Nakashima, B., Stenson, G., Ollerhead, N., and Montevecchi, W. A. 2014a. Bottom-up Regulation of Capelin, a Keystone Forage Species. PLoS One, 9: e87589. https://doi.org/10.1371/journal.pone.0087589

Buren, A. D., Koen-Alonso, M., and Stenson, G. B. 2014b. The role of harp seals, fisheries and food availability in driving the dynamics of northern cod. Marine Ecology Progress Series, 511: 265–284. https://doi.org/10.3354/meps10897

Cadigan, N. G. 2015. A state-space stock assessment model for northern cod, including under-reported catches and variable natural mortality rates. Canadian Journal of Fisheries and Aquatic Sciences, 73: 296–308. https://doi.org/10.1139/cjfas-2015-0047

Cadigan, N. G., Duplisea, D. E., Senay, C., Parent, G. J., Winger, P. D., Linton, B., and Kristinsson, K. 2022. Northwest Atlantic redfish science priorities for managing an enigmatic species complex. Canadian Journal of Fisheries and Aquatic Sciences, 79: 1572–1589. https://doi.org/10.1139/cjfas-2021-0266

Chadwick, E., Brodie, W., Colbourne, E., Clark, D., Gascon, D., and Hurlbut, T. 2007. History of annual multi-species trawl surveys on the Atlantic coast of Canada. Atlantic Zonal Monitoring Program Bulletin, 6: 25–42.

Chapman, E. J., and Byron, C. J. 2018. The flexible application of carrying capacity in ecology. Global Ecology and Conservation, 13: e00365. https://doi.org/10.1016/j.gecco.2017.e00365

Cyr, F., Adamack, A. T., Bélanger, D., Koen-Alonso, M., Mullowney, D., Murphy, H., Regular, P., and Pepin, P. 2025. Environmental control on the productivity of a heavily fished ecosystem. Nature Communications, 16: 5277. https://doi.org/10.1038/s41467-025-60453-6

Dawe, E., Koen-Alonso, M., Chabot, D., Stansbury, D., and Mullowney, D. 2012. Trophic interactions between key predatory fishes and crustaceans: comparison of two Northwest Atlantic systems during a period of ecosystem change. Marine Ecology Progress Series, 469: 233–248. https://doi.org/10.3354/meps10136

Del Monte-Luna, P., Brook, B. W., Zetina-Rejón, M. J., and Cruz-Escalona, V. H. 2004. The carrying capacity of ecosystems. Global Ecology and Biogeography, 13: 485–495. https://doi.org/10.1111/j.1466-822X.2004.00131.x

Fogarty, M., Overholtz, W., and Link, J. 2012. Aggregate surplus production models for demersal fishery resources of the Gulf of Maine. Marine Ecology Progress Series, 459: 247–258. https://doi.org/10.3354/meps09789

Froese, R., Demirel, N., Coro, G., Kleisner, K. M., and Winker, H., 2017. Estimating fisheries reference points from catch and resilience. Fish and Fisheries, 18: 506–526. https://doi.org/10.1111/faf.12190

Fulton, E. A., Link, J. S., Kaplan, I. C., Savina-Rolland, M., Johnson, P., Ainsworth, C., Horne, P., Gorton, R., Gamble, R. J., Smith, A. D., and Smith, D. C. 2011. Lessons in modelling and management of marine ecosystems: the Atlantis experience. Fish and Fisheries, 12: 171–188. https://doi.org/10.1111/j.1467-2979.2011.00412.x

Gamble, R. J., and Link, J. S. 2009. Analyzing the tradeoffs among ecological and fishing effects on an example fish community: A multispecies (fisheries) production model. Ecological Modelling, 220: 2570–2582. https://doi.org/10.1016/j.ecolmodel.2009.06.022

Gomes, M. C., Haedrich, R. L., and Villagarcia, M. G. 1995. Spatial and temporal changes in the groundfish assemblages on the north-east Newfoundland/Labrador shelf, north-west Atlantic, 1978–1991. Fisheries Oceanography, 4: 85–101. https://doi.org/10.1111/j.1365-2419.1995.tb00065.x

Hilborn, R., Walters, C. J., and Ludwig, D. 1995. Sustainable exploitation of renewable resources. Annual Review of Ecology and Systematics, 26: 45–67. https://doi.org/10.1146/annurev.es.26.110195.000401

Hutchings, J. A. 2021. A Primer of Life Histories: Ecology, Evolution, and Application. Oxford University Press. https://doi.org/10.1093/oso/9780198839873.001.0001

Hutchings, J. A. 1996. Spatial and temporal variation in the density of northern cod and a review of hypotheses for the stock’s collapse. Canadian Journal of Fisheries and Aquatic Sciences, 53: 943–962. https://doi.org/10.1139/f96-097

Kell, L. T., Kimoto, A., and Kitakado, T. 2016. Evaluation of the prediction skill of stock assessment using hindcasting. Fisheries Research, 183: 119–127. https://doi.org/10.1016/j.fishres.2016.05.017

Koen-Alonso, M., Lindstrøm, U., and Cuff, A. 2021. Comparative Modeling of Cod-Capelin Dynamics in the Newfoundland-Labrador Shelves and Barents Sea Ecosystems. Frontiers in Marine Science, 8: 579946. https://doi.org/10.3389/fmars.2021.579946

Koen-Alonso, M., Pepin, P., Fogarty, M., and Gamble, R. 2022. Review and Assessment of the Ecosystem Production Potential (EPP) model structure, sensitivity, and its use for fisheries advice in NAFO. NAFO SCR Document, 22/002: 1–52.

Kristensen, K., Nielsen, A., Berg, C., Skaug, H., and Bell, B. 2016. TMB: Automatic Differentiation and Laplace Approximation. Journal of Statistical Software, 70: 1–21. https://doi.org/10.18637/jss.v070.i05

Latour, R. J., Brush, M. J., and Bonzek, C. F. 2003. Toward Ecosystem-Based Fisheries Management: Strategies for Multispecies Modeling and Associated Data Requirements. Fisheries, 28: 10–22. https://doi.org/10.1577/1548-8446(2003)28[10:TEFM]2.0.CO;2

Lear, W. H., 1998. History of Fisheries in the Northwest Atlantic: The 500-year Perspective. Journal of Northwest Atlantic Fishery Science, 23. https://doi.org/10.2960/J.v23.a4

Licandeo, R., Duplisea, D. E., Senay, C., Marentette, J. R., and McAllister, M. K. 2020. Management strategies for spasmodic stocks: A Canadian Atlantic redfish fishery case study. Canadian Journal of Fisheries and Aquatic Sciences, 77: 684–702. https://doi.org/10.1139/cjfas-2019-0210

Link, J., and Sherwood, G. 2019. Feeding, Growth, and Trophic Ecology. In: Rose, G. A. (Ed.), Atlantic Cod: A Bio-Ecology. John Wiley & Sons, pp. 219–286. https://doi.org/10.1002/9781119460701.ch6

Lotka, A. J. 1925. Elements of physical biology. Williams and Wilkins.

Millar, R. B., and Meyer, R. 2000. Non-Linear State Space Modelling of Fisheries Biomass Dynamics by Using Metropolis-Hastings within-Gibbs Sampling. Journal of the Royal Statistical Society: Series C (Applied Statistics), 49: 327–342. https://doi.org/10.1111/1467-9876.00195

Moffitt, E. A., Punt, A. E., Holsman, K., Aydin, K. Y., Ianelli, J. N., and Ortiz, I. 2016. Moving towards ecosystem-based fisheries management: Options for parameterizing multi-species biological reference points. Deep Sea Research Part II: Topical Studies in Oceanography, 134: 350–359. https://doi.org/10.1016/j.dsr2.2015.08.002

Montevecchi, W., and Myers, R., 1997. Centurial and decadal oceanographic influences on changes in northern gannet populations and diets in the north-west Atlantic: implications for climate change. ICES Journal of Marine Science, 54: 608–614. https://doi.org/10.1006/jmsc.1997.0265

Morgan, M. J., Brodie, W., and Kulka, D. 2002. Was over-exploitation the cause of the decline of the American plaice stock off Labrador and northeast Newfoundland? Fisheries Research, 57: 39–49. https://doi.org/10.1016/S0165-7836(01)00331-9

Mueter, F. J., and Megrey, B. A. 2006. Using multi-species surplus production models to estimate ecosystem-level maximum sustainable yields. Fisheries Research, 81: 189–201. https://doi.org/10.1016/j.fishres.2006.07.010

Mullowney, D. R., and Rose, G. A. 2014. Is recovery of northern cod limited by poor feeding? The capelin hypothesis revisited. ICES Journal of Marine Science, 71: 784–793. https://doi.org/10.1093/icesjms/fst188

NAFO, 2022a. Report of the NAFO Commission and its subsidiary bodies (STACTIC and STACFAD), 44th annual meeting of NAFO, 19–23 September 2022, Porto, Portugal. NAFO COM Document, 22/27: 168 p.

NAFO, 2022b. Report of the NAFO Joint Commission–Scientific Council working group on ecosystem approach framework to fisheries management (WG-EAFFM) meeting, 11–12 August 2022. NAFO COM-SC Document, 22–02: 1–22.

NAFO, 2022c. Report of the Scientific Council meeting, 03–16 June 2022. NAFO SCS Document, 22/18: 241 p.

NAFO, 2019. Report of the Scientific Council, 31 May – 13 June 2019, Halifax, Canada. NAFO SCS Doc, 19/20.

Patrick, W. S., and Link, J. S. 2015. Myths that Continue to Impede Progress in Ecosystem-Based Fisheries Management. Fisheries, 40: 155–160. https://doi.org/10.1080/03632415.2015.1024308

Pauly, D., and Froese, R. 2021. MSY needs no epitaph—but it was abused. ICES Journal of Marine Science, 78: 2204–2210. https://doi.org/10.1093/icesjms/fsaa224

Pedersen, E. J., Thompson, P. L., Ball, R. A., Fortin, M.-J., Gouhier, T. C., Link, H., Moritz, C., Nenzen, H., Stanley, R. R., Taranu, Z.E., Gonzalez, A. Guichard, F., and Pepin, P. 2017. Signatures of the collapse and incipient recovery of an overexploited marine ecosystem. Royal Society Open Science, 4: 170215. https://doi.org/10.1098/rsos.170215

Pepin, P., Higdon, J., Koen-Alonso, M., Fogarty, M., and Ollerhead, N. 2014. Application of ecoregion analysis to the identification of Ecosystem Production Units (EPUs) in the NAFO Convention Area. NAFO SCR Document, 14/069, 1–13.

Pikitch, E. K., Santora, C., Babcock, E. A., Bakun, A., Bonfil, R., Conover, D. O., Dayton, P., Doukakis, P., Fluharty, D., Heneman, B., Houde, E. D., Link, J., Livingston, P. A., Mangel, M., McAllister, m. k. Pope, J., and Sainsbury, K. J. 2004. Ecosystem-based Fishery Management. Science, 305(5682). https://doi.org/10.1126/science.1098222

Pinsky, M. L., Eikeset, A. M., Helmerson, C., Bradbury, I. R., Bentzen, P., Morris, C., Gondek-Wyrozemska, A. T., Baalsrud, H. T., Brieuc, M. S. O., Kjesbu, O. S., Godiksen, J. A., Barth, J. M. I., Matschiner, M. Stenseth, N. C., Jakobsen, K. S., Jentoft, S., and Star, B. 2021. Genomic stability through time despite decades of exploitation in cod on both sides of the Atlantic. Proceedings of the National Academy of Sciences, 118: e2025453118. https://doi.org/10.1073/pnas.2025453118

R Core Team. 2021. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Regular, P. M., Buren, A. D., Dwyer, K. S., Cadigan, N. G., Gregory, R. S., Koen-Alonso, M., Rideout, R. M., Robertson, G. J., Robertson, M. D., Stenson, G. B., Wheeland, L. J., and Zhang, F. 2022. Indexing starvation mortality to assess its role in the population regulation of Northern cod. Fisheries Research, 247: 106180. https://doi.org/https://doi.org/10.1016/j.fishres.2021.106180

Robertson, M., Gao, J., Regular, P., Morgan, M., and Zhang, F. 2021. Lagged recovery of fish spatial distributions following a cold-water perturbation. Scientific Reports, 11: 9513. https://doi.org/10.1038/s41598-021-89066-x

Rose, G. 2004. Reconciling overfishing and climate change with stock dynamics of Atlantic cod (Gadus morhua) over 500 years. Canadian Journal of Fisheries and Aquatic Sciences, 61: 1553–1557. https://doi.org/10.1139/f04-173

Rose, G., DeYoung, B., Kulka, D., Goddard, S., and Fletcher, G. 2000. Distribution shifts and overfishing the northern cod (Gadus morhua): A view from the ocean. Canadian Journal of Fisheries and Aquatic Sciences, 57: 644–663. https://doi.org/10.1139/f00-004

Schaefer, M. B. 1954. Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. Inter-American Tropical Tuna Commission Bulletin, 1: 23–56.

Scheffer, M., Barrett, S., Carpenter, S., Folke, C., Green, A. J., Holmgren, M., Hughes, T., Kosten, S., Van de Leemput, I., Nepstad, D., van Nes, E. H., Peeters, E. T. H. M., and Walker, B. 2015. Creating a safe operating space for iconic ecosystems: Manage local stressors to promote resilience to global change. Science, 347: 1317–1319. https://doi.org/10.1126/science.aaa3769

Scheffer, M., Rinaldi, S., Huisman, J., and Weissing, F. J. 2003. Why plankton communities have no equilibrium: solutions to the paradox. Hydrobiologia, 491: 9–18. https://doi.org/10.1023/A:1024404804748

Schijns, R., Froese, R., Hutchings, J. A., and Pauly, D. 2021. Five centuries of cod catches in Eastern Canada. ICES Journal of Marine Science, 78: 2675–2683. https://doi.org/10.1093/icesjms/fsab153

Schindler, D. E., Hilborn, R., Chasco, B., Boatright, C. P., Quinn, T. P., Rogers, L. A., and Webster, M. S. 2010. Population diversity and the portfolio effect in an exploited species. Nature, 465: 609–612. https://doi.org/10.1038/nature09060

Smith, S., and Somerton, G. 1981. STRAP: A User-Oriented Computer Analysis System for Groundfish Research Trawl Survey Data. Canadian Technical Report of Fisheries and Aquatic Sciences, No. 1030.

Templeman, W. 1962. Divisions of cod stocks in the northwest Atlantic. ICNAF Redbook, 3: 79–123.

Thorson, J. T. 2020. Predicting recruitment density dependence and intrinsic growth rate for all fishes worldwide using a data-integrated life-history model. Fish and Fisheries, 21: 237–251. https://doi.org/10.1111/faf.12427

Volterra, V. 1926. Fluctuations in the Abundance of a Species considered Mathematically. Nature, 118: 558–560. https://doi.org/10.1038/118558a0

Walsh, S. J. 1992. Size-Dependent Selection at the Footgear of a Groundfish Survey Trawl. North American Journal of Fisheries Management, 12: 625–633. https://doi.org/10.1577/1548-8675(1992)012%3C0625:SDSATF%3E2.3.CO;2

Winker, H., Carvalho, F., and Kapur, M. 2018. JABBA: Just Another Bayesian Biomass Assessment. Fisheries Research, 204: 275–288. https://doi.org/10.1016/j.fishres.2018.03.010