Density-independent mortality at early life stages increases the probability of overlooking an underlying stock–recruitment relationship

,


Introduction
In the population dynamics of fish, recruitment is both a key process and a major source of uncertainty due to its large variability, often over several orders of magnitude (Shepherd et al., 1990;Hilborn and Walters, 1992).In their 1957 monograph, Ray Beverton and Sidney Holt treated recruitment as densitydependent process (Beverton and Holt, 1957), resulting in a seminal stock-recruitment relationship (SRR) that has since remained widely used in the modelling of fisheries.
SRRs such as the one derived by Beverton and Holt (BH-SRR) consistently fail to capture the recruitment variation observed in stock assessment data, which has given them a reputation of theoretical models that fail to fully reflect recruitment dynamics.We will show that the reason SRRs fail to capture recruitment dynamics is the complexity of recruitment processes and the high V C International Council for the Exploration of the Sea 2021.This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.mortality experienced by pre-recruit individuals over the course of their development from spawning to recruitment.These prerecruit stages can last up to several years during which early-life stages transit through several trophic levels and, typically, different ecological niches and habitats.Consequently, pre-recruits are subject to often very different sources of mortality across stage that may also exhibit differences in their variance.
Mortality in the earliest life stages, i.e. eggs and larvae, is particularly high because of their vulnerability to abiotic and biotic conditions, causing large recruitment variability (Shepherd and Cushing, 1980;Cushing, 1990;Myers and Cadigan, 1993).That is, individuals experience a wide range of environmental and trophic factors that include the physical environment (Stachura et al., 2014;Szuwalski and Hilborn, 2015), predation (Huse et al., 2008;Brosset et al., 2020), and competition (Fortier and Harris, 1989;Ricard et al., 2016;Somarakis et al., 2018).The drivers that best explain recruitment dynamics are often different between stocks of the same species (Trochta et al., 2020) and may vary over time.Consequently, there is often little to no visible relationship between the spawning stock and recruitment in empirical data (Szuwalski et al., 2015;Pierre et al., 2018;Zimmermann et al., 2019), leading some stock assessment practitioners to doubt the existence of an SRR.However, logic and biological reality demand that there is a link between spawning biomass and recruits.Under the common definition of a fish stock, where population dynamics are governed by intrinsic parameters and no relevant migration occurs, the number of eggs spawned depends at least to some extent on the mature component of the stock, and zero spawning biomass produces zero recruits (Haddon, 2010).
Despite the observation of (sometimes very large) variability in recruitment, the relationship between spawning biomass and recruits as major source of density-dependent population regulation has remained a paradigm in the study of population dynamics of marine fish (Zimmermann et al., 2018).A range of different definitions of the SRR exist, but the canonical ones are due to Beverton and Holt (1957) and Ricker (1954).All SRRs rely on the assumption that recruitment is a function of the productivity of the spawning stock and density-dependent mortality of prerecruits due to intra-or inter-cohort competition (Mangel, 2006;Haddon, 2010).In the commonly used forms of the BH-SRR, the parameters linking spawning stock biomass to recruitment are descriptive and may appear arbitrary, which may have misled many readers to overlook that Beverton and Holt derived their model from a mechanistic description of the mortality experienced during the early-life stages until recruitment.Specifically, the model included two distinct mortality parameters, one density-dependent and the other density-independent, which may represent different sources of mortality such as predation or feeding conditions (China and Holzman, 2014;Van Poorten et al., 2018).
Atmospheric and oceanic processes that govern the physical conditions of marine ecosystems [e.g. the North Atlantic Oscillation (Stenseth et al., 2003) and the Subpolar Gyre (Ha ´tu ´n et al., 2005) in the Atlantic Ocean], tend to show long-term oscillations and stochastic variability that subsequently affect recruitment of fish directly (Henderson and Seaby, 2005;Payne et al., 2012), and indirectly by controlling ecosystem productivity (Ha ´tu ´n et al., 2016) and spawning distributions (Miesner and Payne, 2018).Dynamics of predators are also highly variable because their abundance fluctuates in response to environment, food availability, and fishing, but also because prey preference may change or feeding occurs opportunistically.The latter implies that predation pressure is often the stochastic result of spatial overlap of predators and prey in time (Embling et al., 2012;Bartolino et al., 2017;Nikolioudakis et al., 2019;Brosset et al., 2020).
The BH-SRR in its original form with explicit mortality terms enables us to understand how the stage-specific mortality at early-life stages shapes recruitment dynamics and the SRR.It can, thus, provide a useful tool for demonstrating how variation in mortality over time manifests as variation in recruitment and potentially masks the SRR.In this study, we explored how variance in one or both mortality terms affects the recruitment dynamics emerging from the BH-SRR.This allowed us to determine at which point the variance in mortality of early-life stages is so high that it will blur a stock-recruitment relationship to the degree that it seemingly disappears.To compare these simulations with empirical data, we used a set of case studies of Northeast Atlantic fish stocks as baseline for our model parametrization and to explore which degree of variance creates the recruitment patterns that are most similar to a set of stock assessment time series.We complemented these results with empirical estimates of early-life mortality CVs in the literature to illustrate existing variance in various ecosystems and evaluate them against the threshold values emerging from the current study.

Methods
The origin of the Beverton-Holt stock-recruitment relationship Beverton and Holt (1957) derived their SRR from explicit cohort dynamics based on density-dependent survival for a period T from spawning until recruitment; m 2 represents mortality that depends on cohort abundance, e.g.intra-cohort competition for food or space, and m 1 cohort-independent mortality, e.g.predation or environment.Recruitment R is cohort size at time T from an initial cohort size N 0 integrated over the recruitment period T : We define a to be e Àm1T and b to be m2 m1 À Á 1Àe Àm1T ð Þ to obtain (Quinn and Deriso, 1999): Beverton and Holt (1957) defined the parameters differently, leading to a different form of the BH-SRR that can be obtained from Equation (3) by dividing the top and bottom of the righthand side by a; we will work with Equations ( 2) and (3) instead of their equation because of the easier interpretation of the parameters.Because initial cohort size corresponds to number of eggs spawned (or larvae for life-birthing species), recruitment can be directly linked to the spawning stock biomass, typically by assuming a linear relationship N 0 ¼ uB, with fecundity u and abundance or biomass of the spawning stock B. Subsequently, the BH-SRR is: where R(B) is the number of recruits produced by spawning biomass B and a and b are the stock-recruitment parameters.
In the present study, we (i) conducted numerical simulations to explore the role of temporal variability in the mortality parameters on detectability of the BH-SRR; (ii) compared simulated to existing stock-recruitment data from several Northeast Atlantic stocks assessed by the International Council for the Exploration of the Sea (ICES) to explore the type and degree of variability that explain the observed time series best; and (iii) compiled data on early-life stage rates of mortality (or proxies thereof) to provide an estimate of the variability in the rates of mortality.

Numerical simulations
We used numerical simulations to illustrate the dynamics that follow from Equation (2) and the consequences of variation in m 1 and m 2 .The simulation protocol is illustrated in Figure 1.
To simulate variation in the two mortality terms m 1 and m 2 , we assumed that they followed log-normal distributions e lþr with the parameters l ¼ logðx 1;2 Þ and standard deviation r ¼ logðdÞ, where x 1;2 are the mean values for m 1 and m 2 derived from empirical estimates, and the CV is determined by d ¼ 1.00, 1.05, . .., 9.00.The mortality parameters were modelled as a lognormal distribution to ensure strictly positive values, and to include a long tail on the right side of the distribution that reflects low probability events of very high mortality and, thus, the occasional occurrence of recruitment failure.
d ¼ 1 represents no variance, i.e. a CV of 0, whereas d ¼ 9 results in a high CV (with a mean $3 across all simulated time series).For each value of d and iteration of our simulations, we drew m 1 and/or m 2 randomly from these distributions in each year, generating a time series of time-varying mortality.We ran the simulations with four distinct scenarios: (i) variance in m 1 and fixed m 2 ; (ii) variance in m 2 and fixed m 1 ; (iii) variance in m 1 and m 2 (i.e.stochastic annual m 1 and m 2 are independent of each other); (iv) co-variance in m 1 and m 2 (i.e.stochastic annual of m 1 and m 2 are correlated with a Pearson's correlation coefficient of 0.9); and (v) exploration of the entire surface of CV variance in m 1 and m 2 .
Figure 1.Conceptual illustration of the simulation process: variance is introduced for the cohort-independent and cohort-dependent mortality terms according to different scenarios, determining the shape of the log-normal distribution from which annual parameter values are drawn.The emerging stochastic time series of m 1 and m 2 generate a corresponding stock-recruitment time series.This process is repeated 200 times for each scenario and level of variance, and summary statistics are calculated across the resulting set of replicates.
We used simplified population dynamics, where recruitment RðB t Þ represented the input into existing biomass B at time t þ 1, and biomass decreased due mortality, Z: (5) To emulate variation in total mortality due to changes in fishing pressure and natural mortality over time, Z was modelled as autocorrelated, log-normally distributed term with a mean of 0.3, from which stochastic annual mortalities were drawn: Z ¼ 0:3 Á e N 0;0:5 ð Þ ÁAR1 , with AR1 ¼ 0.5.For each set of parameters, we simulated population dynamics for a time period of 50 years, using a random draw from a log-normal distribution (l ¼ À1 and r ¼ 1) as starting value B 0 .Because random draws of annual m 1 , m 2 , and M, as well as the initial stock size B 0 cause stochastic outcomes for each run, we simulated 200 replicates for every value of d in each scenario.The cumulative means of results all started to level off around 50 replicates, thus 200 repetitions were considered sufficient to determine means and confidence intervals of deviation, and statistical significance.
To determine the probability of overlooking an underlying SRR, we used linear regression of log Rt Bt ¼ a þ bB t þ e as a measure of statistically detectable stock-recruitment relationship, assuming no observation error in B t or R t .This approach is a common (Hilborn, 1985;Walters, 1985;Ricard et al., 2016;Szuwalski et al., 2019) and straight-forward way to determine a potential SRR in data.Thus, we took the perspective of a fisheries scientist who does not know how the data were generated.The goal was not to optimize the fitting procedure but to apply a simple, easily comparable metric across a large number of simulated datasets.We used the statistical significance of slope parameter b to indicate the probability of finding a statistically significant SRR, conditioned on the parameters used to simulate the data.Additionally, we calculated the variation in the resulting recruitment with the CV of each recruitment time series, and the deviation of the simulated data from a standard Beverton-Holt model by fitting Equation (4) to the simulated data of each run as the sum of squared residuals (SSQ).
All simulations and analysis were conducted in R 4.0.2(R: Development Core Team, 2020) and the tidyverse packages (Wickham et al., 2019).

Comparison with stock-recruitment data in NE Atlantic stocks
To gauge suitable parameter values for m 1 and m 2 , we estimated the parameters for six commercially and ecologically important fish stocks with comparatively long assessment time series of estimated spawning stock biomass and recruitment for the case study.Parameter estimation was conducted by fitting Equation (2) to the stock-recruitment time series using maximum likelihood estimation, assuming lognormal error distribution.
In a second step, to compare simulated recruitment data with the assessment time series, we applied the same simulation protocol (Figure 1) as in the generic simulations, except the assessed SSB was used as input SSB in all simulations.Our approach was, thus, as follows: first, we fitted Equation (2) to the data to estimate m 1 and m 2 to determine the mechanistic Beverton-Holt SRR with no variation in rates of mortality.Second, we followed the same protocol as described above by using the estimates of m 1 and m 2 as mean value of log-normal distributions and combining them with r ¼ logðdÞ, where we sequentially increased d under the assumptions of scenarios 1 (variation in m 1 and fixed m 2 as estimated), 2 (variation in m 2 and fixed m 1 as estimated), 3 (independent variation in m 1 and m 2 ), and 4 (covariation in m 1 and m 2 ).Third, we selected for each scenario and d the run that deviated the least from the assessment time series based lowest sum of squared residuals.This allowed us to explore the entire range of CVs for scenarios and determine how variance in mortality affects the similarity between simulated and observed stock-recruitment data in these stocks.

Time-varying mortality in assessment time series
To compare the results from our generic simulations and case studies, we conducted a literature review to find existing independent time series of natural and predation mortality that are currently used in stock assessments, complemented with variance in space of predator abundance.

Results
Estimates of the cohort-independent mortality parameter m 1 from the stock assessment time series of six stocks were consistently higher than m 2 by, in average, a factor 10 (Table 1).The ratio between m 1 and m 2 ranged between 2.5 for Norwegian springspawning herring to 28 for North Sea cod.Based on the means of the estimates, the mean values for m 1 and m 2 in the subsequent generic simulations were set to m 1 ¼ 0:2 and m 2 ¼ 0:02.When assuming no variation in mortality, these parameter values result in ca. 4 percent of the initial cohort size remaining at recruitment.Trajectories for the selected parameter values are shown in Supplementary Figure S1 together with examples of higher/lower m 1 and m 2 .The experienced natural mortality rates M given the same set of parameter values are illustrated in Supplementary Figure S2.Increased variance in natural mortality blurs the stockrecruitment relationship Although higher CVs in mortality coefficients increased in general the probability of overlooking an underlying SRR, variation in the cohort-independent mortality m 1 had a much stronger effect on the statistical significance (linear regression slope b with mean p-value < 0.05) of the SRR than cohort-dependent mortality m 2 (Figures 2 and 3a).Statistical non-significance tended to be reached around a CV ¼ 0.75 to 1.00 for cohort-independent mortality, independent of whether m 1 alone was varied and m 2 remained fixed, or whether both m 1 and m 2 were varied.In the latter case, independence or covariance between the two mortality parameters did not change the results.Variance in cohortdependent mortality alone, the SRR remained detectable for most simulated time series for all levels of CV, resulting in now significant increase above the threshold of p ¼ 0.05 across all simulated time series (Figure 2).Conversely, the effects on total observed variation in recruitment and deviation from a deterministic SRR were more pronounced when variation in cohort-dependent mortality m 2 was introduced than with variation in m 1 (Figure 3b and c).Although a higher CV in m 1 or m 2 resulted in all scenarios in a gradual and clear increase in variation of the resulting recruitment (Supplementary Figure S3) and its deviation from the deterministic BH-SRR (Supplementary Figure S4), the impact was more strongly driven by cohort-dependent mortality m 2 .In contrast to the effects on probability of overlooking the underlying SRR, the results showed compounding effects when m 1 and m 2 are allowed to vary at the same time, and the deviation from the SRR was strongest when there was covariance between the two mortality parameters.
The same patterns are found for the entire surface of sum of squares deviation and statistical significance (Supplementary Figure S5): whereas the level of statistical significance was largely determined by the CV in density-independent mortality (m 1 ), quantitative changes in the deviation from the deterministic BH-SRR were dominated by the CV in density-dependent mortality (m 2 ).The latter also applied for the CV in the recruitment time series, underlining the stronger quantitative effect of m 2 on the recruitment variation, although the difference was less pronounced than for the residual sum of squares deviation.
The specific parameter values used influence the results quantitatively but not qualitatively.The mean values of m 1 and m 2 and the ratio between them are particularly important for the emerging SRR, since they scale total mortality between spawning and recruitment, and determine the distribution of the mortality parameters.Higher means of m 1 and m 2 decreased the survival of recruits and shifted the SRR towards lower stock and recruitment values, reducing detectability of the SRR (Supplementary Figure S6).However, the pattern of how the different scenarios affect the detectability remained independent of the specific mean mortalities.Differences in mean m 1 and m 2 had very little impact on the CV of the simulated recruitment time series (Supplementary Figure S7).Furthermore, the results remain also consistent when using equal values for m 1 and m 2 , i.e. a ratio of 1 (Supplementary Figures S8 and S9).This shows that the main results of our study are not sensitive to the specific parametrization of m 1 and m 2 .
In Supplementary Figure S10, we show examples of simulated stock-recruitment time series with no (CV ¼ 0), intermediate (CV ¼ 0.2-0.6), and high (CV ¼ 0.8-1.2) variance in m 1 and/or m 2 .Compared to the default case (CV ¼ 0), increasing CVs resulted in all cases in a broader distribution of recruitment values but also in a larger spread in the range of stock biomasses, especially regarding the upper biomass levels.The latter effect was, however, much more pronounced in scenarios 2 and 3 with variance in density-dependent mortality (m 2 ).On the other hand, variance in density-independent mortality (m 1 ) tended to generate a larger distribution of recruitment independent of biomass compared to the effects of m 2 on recruitment, which were more proportional to stock biomass where variance in recruitment increased with biomass.

Comparison with stock-recruitment data in NE Atlantic stocks
When comparing SRRs with increasing variance to six datasets of North Atlantic stocks, the results showed a clear tendency for simulated data to be most similar to observed data when CVs are at or above 0.5 (Supplementary Figure S11).In most cases, the range of CVs resulting in simulated recruitment most similar to the assessed time series aligned well among scenarios, and typically the quality of fit started to become worse than the zerovariance model when CV > 1. Notable exceptions are Norwegian spring-spawning herring and North Sea cod where higher CVs than 1 in m 1 and m 2 , respectively, did not result in a decrease of model performance but continued to generate output that was more similar to the data than most other models.North Sea cod, Barents Sea capelin, and sardine were the three stocks where one scenario (variation in m 2 ) outperformed the other scenarios.

Empirical variation in the rates of mortality
Although empirical time series on variation in mortality at earlylife stages are scarce, we were able to collate CV estimates for a range of different species and ecosystems (Table 2).The CVs of total natural mortality at age 0 used in stock assessments were overall relatively low, ranging from CV ¼ 0.04 (North Sea whiting) to CV ¼ 0.46 (Northeast Arctic cod).In comparison, the available CVs of specific cannibalism/predation mortalities used in stock assessments were substantially higher: 1.24 for cannibalism in Northeast Arctic cod and 1.42 for predation of Northeast Arctic cod on haddock.Existing time series of spatial variation in predator presence showed particularly high CVs, represented by various seabird species in the California current with CVs of 1.32-4.12.

Discussion
Although fundamental biological principles require that recruitment is linked to the spawning stock size (Myers and Barrowman, 1996;Vert-pre et al., 2013;Zimmermann et al., 2018), in most stocks the SRR explains only a minor proportion of recruitment variability (Cury et al., 2014;Munch et al., 2018).As a result, the empirical relevance of SRRs has been increasingly dismissed despite their continued conceptual relevance (Rose et al., 2001;Szuwalski et al., 2015;Szuwalski et al., 2019) [which has been reflected in stock assessments when, for instance, steepness is set at or close to 1 and thus recruitment is assumed to be entirely or largely independent of stock size, see e.g.Aires-da-Silva and Maunder (2007) and Maunder (2012)].Building on the SRR derivation in Beverton and Holt's (1957) monograph, our results provide a synthesis of these divergent perspectives.We demonstrated how the original derivation of the BH-SRR, with mortality separated into explicit cohort-dependent (driven by competition among the pre-recruits) and cohort-independent terms (mortality caused by environmental variation and predation) can help explain observed stock-recruitment patterns, particularly the frequent lack of a clear SRR.Using the BH-SRR, we showed that (i) an underlying stock-recruitment relationship is easily masked by variance in pre-recruit mortality and (ii) even a low to intermediate CV of 0.75 in cohort-independent mortality can conceal an existing SRR.Despite the potentially large recruitment variation introduced by variation in cohort-dependent mortality, we empirically estimated cohort-independent mortality to be higher by in average a factor 10 and found that variance in cohort-independent mortality had a stronger effect on the statistical significance of the stock-recruitment relationship, masking the relationship at lower levels of variation.This suggests that cohort-independent sources of mortality such as predation may be the main cause for the weak empirical evidence for SRRs and suggests a direction for future empirical work.
To provide a perspective on the levels of variance that might be expected from empirical data information, we compared simulations with added variation in density-dependent and densityindependent mortality in the BH-SRR coefficients with existing stock-recruitment data.In all six case studies, estimated cohortindependent mortality was higher than cohort-dependent mortality, however with a substantial range in the ratio between both mortalities.The difference between the two sources of mortalities was estimated to be lowest for the two herring stocks and NEA cod, suggesting that both cohort-dependent and -independent drivers play a relevant role for recruitment.In contrast, the ratio between the two mortality parameters was highest for Barents Sea capelin and North Sea cod, indicating that those two stocks might be strongly governed by external drivers of recruitment.This reflects largely existing knowledge about recruitment dynamics in these stocks that include some of the longest available time series of stock-recruitment data and a wealth of recruitment studies, although the precise causes for recruitment variation have remained surprisingly elusive.However, most studies link the observed recruitment variation foremost to a mixture of densitydependent and -independent variables, particularly the spatiotemporal dynamics of environmental conditions and predation experienced during larval drift, such as in Norwegian springspawning herring (Skagseth et al., 2015;Garcia et al., 2020;Tiedemann et al., 2020), North East Arctic cod (Yaragina et al., 2009;Endo et al., 2020), North Sea cod (Akimova et al., 2016), capelin (Gjøsaeter et al., 2016), and sardine (Garrido et al., 2017).The complexity of the underlying environmental processes, feedbacks between them and their possibly changing relevance over time poses a challenge for pinpointing clear causalities.This is reflected in the comparison of our BH-SRR generated stock-recruitment data with the observed stock-recruitment data, as our results indicate that CV levels close to the threshold that make an Natural mortality and predation CVs are based on modelled time series of annual natural mortality or predator consumption at age 0 or 1 using observations from stomach samples and a multispecies model for North Sea stocks (ICES, 2018e).Natural mortality is the sum of a fixed mortality component and annual estimates of variable predation mortality.We included only periods with time-variant natural mortality (thus, constant natural mortality in North Sea herring before 1975 was excluded).*Data by J. Santora (pers. comm.), see also Santora and Sydeman (2015) and (Sydeman et al., 2015).
SRR difficult to detect may result in simulated stock-recruitment pattern most similar to observed ones.Nevertheless, no clear patterns as to specific levels of CVs or the type of mortality emerged across stocks, preventing any clear conclusions.
In most of the 11 stocks where time-varying natural mortality was included in stock assessment (Table 1), the CVs of natural mortality were below the threshold for detection of an SRR based on our generic simulations.However, the time-varying natural mortalities used in these assessments are the sum of a constant base mortality and varying predation mortality, often just based on one specific predator or, in the case of the North Sea, estimated through a multispecies model.In addition, only time series for ages 0 and 1 are included in these assessments, omitting the earlier life stages that experience the highest mortality and, possibly, largest variation in mortality (Bogstad et al., 2016).Thus, the empirical CVs of natural mortality provide an incomplete picture of the total variance in mortality and may be underestimates.Empirical CVs of consumed pre-recruits by specific predators as well as spatial predator distribution were substantially higher, but it is not possible to determine how these indirect estimates of predation translate into actual variation in pre-recruit mortality.The specific impact of predation may also depend on the functional relationship between prey abundance and predator, which may range from density-independent to various forms of density dependence.Moreover, the impact on total mortality of densityindependent predation mortality such as from seabirds may be more significant when prey abundance is low (Saraux et al., 2020).The limited availability of mortality time series or estimates of the variance of mortality underlines that time-varying mortality is rarely estimated or used within assessment frameworks, and that there is substantial knowledge gap when it comes to empirical estimates of pre-recruit mortality.This applies especially to the earliest life stages where mortality is particularly high and variable (Bogstad et al., 2016) but empirical information difficult to obtain.However, in the spirit of Beverton and Holt, because something is difficult does not mean it should not be done.
Recruitment is a key determinant of the status of a stock and can directly affect changes in biomass and management success (Pierre et al., 2018;Zimmermann and Werner, 2019).Understanding recruitment and its relationship with the spawning stock is therefore a key component in stock assessment and management advice (Mangel et al., 2013;He and Field, 2019).SRRs play a key role for both tactical and strategic management (Needle, 2001) through short-term for forecasts for management advice and Management Strategy Evaluations to set reference points (Punt et al., 2016).Lacking better alternatives, most management approaches assume classic SRRs such as Beverton-Holt or long-term means of recruitment (or an ensemble of SRRs), sometimes complemented with some form of stochasticity.We demonstrated the large impact of variation in pre-recruit mortality, underlining that better knowledge of CVs of pre-recruit mortality could provide a more complete picture of the possible recruitment dynamics as well as the resulting uncertainty.Because the dynamics of many commercial stocks are determined by recruitment, a better incorporation of pre-recruit mortality into assessment frameworks could contribute to improved accuracy of forecasts and, thus, more sustainable management.
The interplay of cohort-independent and cohort-dependent mortality on the detectability of an SRR aligns with other concepts that attempt to explain recruitment variability through spatial and temporal overlap of early-life stages with environmental conditions and food availability, such as the match-mismatch (Cushing, 1990) or optimal window (Cury and Roy, 1989) hypotheses.Density dependence and independence are often intertwined, as density-dependent processes such as food competition depend on density-independent processes such as ecosystem productivity.Notably, predation can be both density-independent, e.g. when predation mortality is caused by spatio-temporal overlap of an opportunistic predator, or partially density-dependent when a predator responds to the density of its prey.Moreover, predation pressure depends on the density and distribution of the predator, which vary over time and space and contribute to variation in mortality (Akimova et al., 2019).The relevance of densitydependent and -independent processes may also vary over time in response to stock size, and thus the state of a stock may partially determine the degree to which the different mortality sources affect overall recruitment variation.This implies that fishing pressure can also impact recruitment dynamics by reducing the spawning stock to low levels where density-dependent processes become less relevant.Nevertheless, our results have consistently shown that although cohort-dependent mortality may shape the overall recruitment variability and act as major determinant of cohort strength, the resulting recruitment variability remains largely proportional to spawning stock size.Density-independent sources of mortality, on the other hand, masked the link between spawning stock and recruitment more strongly, even if their total quantitative impact may be smaller.The key reason is that density-independent mortality acts in a constant per capita manner while density-dependent mortality declines as cohort size declines, explaining the relative importance of the two kinds of mortality.
Typically, both cohort-dependent and -independent sources of mortality are important for recruitment dynamics, so that their effects may be compounded and confounded.The relative importance of cohort-dependent and -independent mortality may vary over time, i.e. inter-annually as well as from stage to stage between spawning and recruitment.Detailed case studies of NEA cod and haddock, for instance, show that the variation in mortality rates and whether mortality is mainly density-dependent or density-independent differs among stages of pre-recruits (Bogstad et al., 2013;Langangen et al., 2014).These complex layers of time-varying variation may also explain why significant relationships between recruitment and environmental drivers are found in some periods and not in others (Bogstad et al., 2013).Robust statistical evidence for clear links between environment and recruitment has therefore remained elusive (Myers, 1998).
In our study, we did not account for observation and process errors in stock-recruitment data and their effects on the detectability of an SRR.Typically, stock-recruitment data, such as we used for the empirical analysis, do not contain direct observations but are model outputs constrained by their own assumptions and uncertainties (Brooks and Deroba, 2015).This becomes particularly relevant when comparing datasets generated with the different assessment models (Dickey-Collas et al., 2015).Here, we assumed no observation errors when determining the detectability of SRRs.However, observation errors may add another layer of uncertainty, as stock and recruitment data may be unreliable and contain errors that distort the SRR (Walters and Ludwig, 1981).Recruitment data are difficult to estimate and largely backcalculated from cohort sizes at later ages (sometimes complemented with indices from egg, larvae or 0-group surveys), resulting in large uncertainty.In addition, spawning stock biomass may be an inaccurate proxy for spawning output (Marshall et al., 1998;Kell et al., 2016), because quantity and quality of eggs spawned are affected by other factors than just biomass of mature individuals (Subbey et al., 2014), such as demographic structure and body size (Hixon et al., 2014;Stige et al., 2017;Barneche et al., 2018), energy allocation and maternal effects (Berkeley et al., 2004;McBride et al., 2015), and skipped spawning (Skjaeraasen et al., 2012(Skjaeraasen et al., , 2020)).This will further distort of the stock-recruitment relationship.Another often overlooked driver of recruitment variation may be the degree to which a stock and pre-recruits tend to concentrate spatially (Iles and Beverton, 2000), determining the vulnerability to both external drivers as well as density-dependent processes.
Beverton and Holt's SRR used in our study contributes a more mechanistic understanding of recruitment processes, but omits important life-history dynamics such as the importance of growth for survival and the potentially large differences in causes and extent of mortality experienced across different pre-recruit stages.
Here we focused solely on mortality of pre-recruits as direct driver of variability, although survival in these life stages depends strongly on larval size that defines the possible prey spectrum as well as the predation risk.This means that mortality from low food availability and predation can both be linked to larval growth.An extension of our analysis could therefore account explicitly for growth as a function of its environment (temperature, food availability), which would allow for more mechanistic perspective on the causes of mortality (e.g.Mangel et al., 2010).In this paper, we demonstrated how the mechanistic version of the BH-SRR can be useful to explore the effects of different sources of pre-recruit mortality, illustrating how they can distort and mask the underlying stock-recruitment relationship.Our results highlight the relevance of a multi-stage perspective on recruitment that accounts explicitly for the time between spawning and recruitment and the different (variable) mortalities the sequential life stages encounter.Multi-stage recruitment models are thus an important tool to improve our knowledge on recruitment dynamics and SRR estimates (Brooks et al., 2019).Doing so may ultimately improve forecasts of stock dynamics and therefore management advice, since predicting recruitment has remained a key cause of error and uncertainty in fisheries management.Incorporating variation in natural mortality and predicting recruitment with models that capture the multi-stage nature (Brooks et al., 2019) and nonlinearity (Munch et al., 2018;Sguotti et al., 2020) of pre-recruit dynamics will contribute to better assessment estimates and management advice.
In conclusion, Beverton and Holt's seminal work in their 1957 monograph continues to provide important, novel insights and useful tools for these future advancements in fisheries research.Our study highlights that the SRR developed by Beverton and Holt and the underlying principles governing population dynamics remain valid and important, even more so when an ecosystem-based approach to management accounts more explicitly for the multiple drivers of recruitment dynamics.

Figure 2 .
Figure 2. Statistical significance (p-value) of linear regression slope for different levels of variance in cohort-independent mortality m 1 (a), cohort-dependent mortality m 2 (b), in m 1 and m 2 independently (c), and in m 1 and m 2 with positive covariation (d).Each CV represents the variation of one simulated time series of m 1 and/or m 2 .Dots represent each a simulation run with a specific CV value, solid lines are fitted GAM smoothers with 95% confidence intervals (shaded areas, too small to be visible for a wide range of values), dashed lines indicate the significance threshold at p ¼ 0.05.

Figure 3 .
Figure 3. Boxplots summarizing the probability of overlooking the underlying SRR (a), the CV of recruitment time series (b) and the residual SSQ deviation from the deterministic BH-SRR (c) for three values of CVs of the time-varying mortality coefficients (0.5, 1, and 2) and four scenarios [variance in cohort-independent mortality (m 1 Þ, in cohort-dependent mortality (m 2 ), both mortality coefficients (m 1 and m 2 ) independently (without covariance), and with covariance between the coefficients].Each data point represents one simulated time series with a CV of m 1 and/or m 2 and the corresponding summary statistic.Boxplots show median (solid black line), first and third quartiles (lower and upper bounds of boxes), 1.5 times the interquartile range (whiskers), and the outliers outside of this range (dots), over all bootstrapped runs for each CV and scenario.The dashed line indicates the significance threshold at p ¼ 0.05.

Table 1 .
Parameter values for m 1 and m 2 estimated from stockrecruitment time series of six stocks and corresponding ratio between m 1 and m 2 .

Table 2 .
CVs of mortality or predation pressure. *