The dynamics of fish populations at low abundance and prospects for rebuilding and recovery

Abstract

Previous meta-analysis of spawner–recruit relationships suggested that depensatory behaviour is uncommon, and stocks pushed to low abundance are unlikely to suffer decreases in recruitment more severe than would be expected based on the decline in spawning stock. Using an updated database that has over 100 stocks that were depleted to less than 20% of their maximum observed stock size, we tested for depensatory behaviour in both total surplus production and recruitment and we also examined the probability of stock increase as a function of stock size and fishing pressure. The number of stocks that showed a significant improvement with depensatory models was less than that expected by chance. Hierarchical meta-analysis showed that the majority of the evidence was for no depensatory behaviour but could not rule out depensation at very low stock sizes. Stocks that are depleted to low abundance are expected to rebuild when fishing pressure is reduced if the environment has not changed but there is considerable evidence that the majority of fish stocks are impacted by changes in productivity regimes. Nevertheless, if stocks are very heavily depleted and fishing pressure is not reduced to quite low levels, the expected recovery time is both uncertain and long. Very low abundance should clearly be avoided for many reasons and the range of abundance where depensation cannot be ruled out is well below commonly adopted limit reference points.

Introduction

It is 100 years since the publication of Johan Hjort's classic work “Fluctuations in the great fisheries of northern Europe viewed in the light of biological research” (Hjort, 1914). In this monograph, Hjort considered what was known at the time about the causes of fluctuations in fish stocks and specifically addressed what was known about the formation of year-class strengths in relation to spawning stock abundance. Hjort found little relation between the two and stated “A rich spawning may produce a year class poor in numbers, while a large year class may have its origin in a year when the spawning was at its lowest.” In the last 100 years, we have acquired considerably more information about the behaviour of populations at low abundance and are able to expand considerably upon Hjort's understanding.

Perhaps, the most influential fisheries event in the last 50 years was the collapse and closure of the northern cod fishery in Eastern Canada in the early 1990s (Hutchings and Myers, 1994; Rice, 2006). Tens of thousands of people were put out of work, and the developing system of fisheries assessment and science was brought into question. A key issue raised by the cod collapse is the dynamics of stocks at low densities, and the potential for multiple stable states, where once pushed to a low abundance, populations might be stuck in that condition despite reduced fishing pressure. Walters and Kitchell (2001) used trophic models to show how such multiple stable states could arise, and Worm and Myers (2003) documented many cases of where the depletion of demersal fish species led to ecosystem changes often characterized by increases in invertebrate abundance. Others have shown that ecosystem state changes associated with fishing on other parts of the ecosystem (Anderson and Piatt, 1999; Utne-Palm et al., 2010).

The 1990s saw a renewed interest in the relationship between spawning stock and abundance, and Myers et al. (1994) used a compilation of over 100 spawner­–recruit relationships to argue that there was good evidence that lower spawning stocks produced lower average recruitments. At about the same time, there was work on whether rates of increase declined at low densities (known as depensation), a necessary condition for multiple stable states, and Myers et al. (1995) and Liermann and Hilborn (1997) found little evidence of depensation. Liermann and Hilborn (1997) state “we … found that, for all of the taxa, the most likely values fell close to or within the range of no depensation”. Using the same database as Myers et al. (1994), these two papers showed that per capita recruitment rates continually increased as the density of spawning stock was reduced.

The dynamics of populations at low densities play a crucial role in understanding extinction processes, as well as in the long-term dynamics of populations and ecosystems (Dulvy et al., 2003; Reynolds et al., 2005). Density-dependent rates of increase can be classified into two types, compensatory rates in which the per capita growth rate of a population increases as density decreases and depensatory rates where the per capita growth rate decreases as density decreases. Common mechanisms that lead to compensation include more food per individual at low densities, each individual more likely to find refuge habitat from predators at low densities and lower disease transmission rates at low densities. Compensatory rates lead to population stability, although very strong compensation, when combined with discrete time dynamics, can produce chaotic behaviour (May, 1976). Compensatory rates of increase have been identified in a wide range of species and indeed almost all discussion of density-dependence assumes compensatory processes.

Depensatory processes can include predation if the number killed by predators is largely independent of the abundance of the population (Roemer et al., 2002), group facilitation where at low densities such mechanisms as predator detection are less effective, habitat conditioning where individuals improve habitat quality for others, or the reduced probability of finding mates at low densities (often called the Allee effect; Liermann and Hilborn, 2001). Depensatory processes can lead to much more complex dynamics and have been the object of considerable speculation and exploration within the theoretical literature (Courchamp et al., 1999; Stephens and Sutherland, 1999). Multiple stable states generally require the population rate of increase to go negative at low densities, a condition known as “critical” depensation.

Both the studies of depensation (Myers et al., 1995; Liermann and Hilborn, 1997) used a dataset that contained estimates of the population size and the subsequent recruitment for hundreds of fish stocks. Myers et al. tested for the presence of depensation using a tradition p-value hypothesis test, comparing a model with and without depensation. They found little evidence for depensation, with no more stocks showing significance than would be expected by chance and concluded that “depensatory dynamics are not apparent for fish populations at the levels studied”. Liermann and Hilborn used a Bayesian approach to calculate the intensity of depensatory processes and for the dataset as a whole found little support for depensatory processes.

A major limitation of the data used in both studies was the few number of populations that were observed at low abundance. Of the 128 datasets examined, only 26 had high statistical power. The dataset used in these two studies has recently been updated to include 15 more years of data for fish stocks that were often driven to low densities (Ricard et al., 2012), thus providing a more powerful dataset to determining the frequency and intensity of depensation. In the intervening 17 years since the publication of the first paper, much more attention has been directed towards concerns about fish stock abundance and the impact of fishing on individual stocks and marine ecosystems.

The collapse of eastern Canadian cod (and other groundfish) also led to a focus on the frequency and intensity of fish stock collapses. Meta-analysis of changes in stock abundance showed that large declines in fish stock abundance were quite common and that some taxa (particularly clupeids) recovered about as frequently as they declined (Hilborn, 1997; Hutchings, 2000). Other taxa (especially gadids) showed much less frequent recovery. Hutchings and Reynolds (2004) argued that “reductions in fishing pressure, although clearly necessary for population recovery, are often insufficient” and heightened concern that populations pushed to low densities may not recover even when fishing pressure is reduced. More recent papers (Hutchings et al., 2012a, b) have explored life history correlates with extinction risk and recovery failure.

Mace (2004) criticized Hutchings' (2000) conclusions, pointing out that he had not taken into account whether the fishing pressure had been reduced enough to expect recovery and that the dataseries terminated in the early 1990s which did not allow enough time to see if recovery would actually take place if fishing pressure were sufficiently reduced.

Using an updated version of Myers' 1990s database (Ricard et al., 2012) that generally extends time-series at least 10 years beyond those in Myers' database, a number of results have emerged that provide a much better understanding of the dynamics of fish stocks. Between the early 1990s and the mid to late 2000s, we saw many stocks pushed to low abundance and many stocks recover from overfishing.

Keith and Hutchings (2012) used a non-parametric model to explore evidence for depensatory and compensatory dynamics in spawner–recruit data using the updated dataset and largely confirmed the earlier analysis. “…there is strong evidence of an Allee effect for only 1 of 104 species at their lowest recorded SSB … there is weak evidence of an Allee effect in another three species.”

Neubauer et al. (2013) used the updated database and found that stocks did indeed recover if fishing pressure was reduced and that once the exploitation rate was reduced to the level that would produce long-term maximum sustainable yield (often called FMSY), recovery was generally expected within 20 years, although populations with unusually low intrinsic rates of increase or those depleted to very low levels are expected to take longer (Figure 3B of their paper). This is a considerably different conclusion than the earlier work by Hutchings and Reynolds (2004). Stochastic changes in productivity also appear to play a large role in the recovery process. The chance of recovery to B_MSY within 10 years for a mildly depleted stock (40% of B_MSY) is nearly 50% even at fishing mortality rates up to 1.5 times F_MSY (Figure 2D of their paper), though of course a stock fished at this rate would not be predicted to remain above B_MSY for long.

The conclusion of Myers et al. (1994) that low spawning stock led to low recruitment was challenged by Gilbert (1997), who argued that recruitment in marine fish stocks is largely driven by environmental conditions that come in regimes of good and bad. Gilbert suggested that when environmental conditions move from good to bad, recruitment declines, and as a result spawning stock declines, thus giving the impression that lower spawning stocks lead to lower recruitment, when the causal mechanism is instead that lower recruitment leads to lower spawning stocks. When environment changes from bad to good, recruitment increases, spawning stock increases, and again it looks like larger spawning stocks lead to more recruitment.

Vert-pre et al. (2013) and Vert-pre (2013) explored the Gilbert hypothesis by examining whether changes in recruitment and surplus production were more tightly linked to biomass changes or shifting environmental regimes. They defined regimes as discrete changes in average surplus production or recruitment. For both measures of stock productivity, they found far more evidence that productivity changes are driven by periodic environmental shifts than by changes in adult biomass. For only 15% of stocks were changes in recruitment best explained by changes in biomass and for 24% of stocks changes in surplus production were best explained by changes in biomass. Table 1 summarizes these results with respect to four hypotheses—that recruitment or surplus production is driven by environmental regimes, that it is driven by stock biomass, that it is a result of a mix of environmental regimes and biomass, or that it is totally random.

Table 1.

The proportion of fish stocks for which different mechanisms best explain recruitment and surplus production from Vert-pre et al. (2013) and Vert-pre (2013).

Causal mechanism	Recruitment (%)	Surplus production (%)
Environmental regimes	50	27
Stock biomass	15	24
A mixed hypothesis	17	45
Random	18	4

Although the mixed model does include an effect of stock biomass, they found that 80% of the variability in surplus production in the mixed-model cases was due to changes in regime, and only 20% from changes in biomass. Their results provide strong evidence that biomass only weakly affects the recruitment or production and that the assumption that stock productivity or recruitment will rebuild by increasing stock size is poorly supported.

Furthermore, the changes in productivity are quite large. Figure 1 shows the frequency of regime changes, scaled so that +1 represents cases where the productivity increased by the mean productivity value over the time-series and −1 represent cases where productivity decreased by the mean productivity.

Figure 1.

Frequency of productivity changes from either regime or mixed models. The gap around zero on the x-axis occurs because the algorithm to detect regime changes does not identify small changes. From Vert-pre et al. (2013).

In this paper, we recast the question of low density dynamics in several ways. First, we formulate a depensation model that can be used for hierarchic meta-analysis and is derived from first principles. Second, we extend this approach from stock–recruitment (the focus of the work by Myers et al., 1995 and Liermann and Hilborn, 1997) to surplus production. This extension recognizes the possibility that depensatory processes may operate not only through changes in juvenile survival (i.e. through the stock–recruitment relationship), but also through other vital rates which contribute to population growth, including adult survival and somatic growth. Both models are applied to data from the updated RAM Legacy Stock Assessment Database (Ricard et al., 2012). Finally, we explore the empirical relationship between stock size, fishing pressure, and probability of stock size increasing.

Material and methods

Models of depensatory recruitment

Previous models of depensation used by Myers et al. (1995) and Liermann and Hilborn (1997) used ad hoc models that modified compensatory models to reduce the rate of increase at lower densities. We derived a model based on first principles, using the Allee effect as an example of a depensatory process, but we believe that this formulation can capture a range of possible processes. In the Allee effect, the probability of being mated or eggs being fertilized declines at low abundance because at low abundance females may not encounter a male. We assume that contact between individuals of the two sexes is a random encounter process, with the probability of any individual female encountering any individual male being p. If N is the number of males in the population then the probability that an individual female will encounter no males and thus not be mated is the product of a Bernoulli trial with each male and this will in turn be approximately Poisson-distributed. Thus, the probably of not being mated can be assumed to be the zero class of this Poisson distribution [Equation (1)].

$$\begin{array}{cc}Pr(\text{no}\text{mating})& =\frac{\text{exp(}-pN\text{)(}pN{\text{)}}^0}{0!}\\ & =\frac{\text{exp(}-pN\text{)}1}{1}=\text{exp(}-pN\text{).}\end{array}$$

(1)

The probability of being mated is simply 1 minus the probability of not being mated [Equation (2)].

$$Pr(\text{being mated})=1-\text{exp(}-pN\text{)}$$

(2)

This function is parameterized to define the depensation parameter as the population size at which 50% of the population is mated, N₅₀, which is used to produce the depensation parameter q that measures N₅₀ relative to the maximum population size [Equation (3)].

$$0.5=1-\text{exp(}-p{N}_{50}\text{)}\to p=-\frac{\text{log}(0.5)}{N_{50}}$$

$$\begin{array}{cc}Pr(\text{being mated})& =1-\text{exp}\left(\text{log(}0.5\text{)}\frac{N}{N_{50}}\right)\\ & =1-\text{exp}\left(\frac{\text{log(}0.5\text{)}N}{q\text{max(}N\text{)}}\right).\end{array}$$

(3)

The choice of 50% is totally arbitrary and has no impact on the interpretation of the analysis. If we had chosen 20% then the posterior distributions would have shifted farther to the right, but when actually applied in a model the choice of N₅₀ could have been replaced with any other level with no difference in conclusions about the intensity of depensation.

In most marine fish stocks, the parental stock in spawner–recruit analysis is normally referred to as the spawning-stock biomass (B_t), calculated as the number of individuals at age times the weight of females at age, and is often adjusted by the maturity-at-age. We thus define a depensatory impact that reflects the fraction of females who are mated [Equation (4)]. Although it would be better to use a biomass reference point such as the unfished equilibrium spawning population instead of the largest observed stock size, this reference point is often not available or very poorly estimated.

$$D_t=1-\text{exp}\left(\frac{\text{log(}0.5\text{)}{B}_t}{q\text{max(}B\text{)}}\right)$$

(4)

Although this formulation of depensation is derived from the Allee effect and finding mates, it captures a range of depensatory processes with a single parameter q that has a straightforward biological interpretation and is amenable to hierarchical meta-analysis. Scaling depensation to the maximum observed biomass is done to scale N₅₀ relative to the scale of the population (which may run from tens of tonnes to hundreds of thousands).

The relationship between the spawning biomass of the stock and the subsequent production of offspring, known as the stock–recruitment relationship in fisheries, has commonly used one of the two models. The Beverton–Holt model is derived from the assumption that compensatory processes from birth until measured as recruitment (often 1 year) is continuous. That is, the instantaneous rate of survival of a cohort depends linearly on its densities. The Ricker model assumes, instead, that the rate of survival depends on the initial cohort size. The derivation of both models is shown in Hilborn and Walters (1992), and Walters and Korman (1999) demonstrate that the Beverton–Holt model follows directly from arena foraging theory that evaluates foraging strategies in light of predation risk. Although the best fit for these two models for many datasets is often quite similar, there are consequential differences for population dynamics and management reference points such as maximum sustainable yield.

A more general model is the Deriso stock–recruitment function [Equation (5)], introduced in Deriso (1980), where B_t is the biomass at time t and R_t is the number of recruits generated from B_t. The Deriso model has a flexible range of shapes and is able to recreate the Beverton-Holt curve with g = 1 and the Ricker curve with g → ∞.

$$R_t=\frac{a{B}_t}{{(1+b{B}_t)}^g}\text{exp(}{ϵ}_t\text{)}$$

(5)

To include a depensatory model with the Deriso model, we multiply the biomass terms in the Deriso curve by the depensation process [Equation (4)] to obtain the depensatory stock recruitment curve shown in Equation (6). The depensation parameter q is bounded between (0,1), and q ≈ 0 indicates no depensatory dynamics since recruitment simplifies to the Deriso model.

$$\begin{array}{c}{R}_t=\frac{a{B}_t{D}_t}{{(1+b{B}_t{D}_t)}^g}\text{exp(}{ϵ}_t\text{)}\\ \\ D_t=1-\text{exp}\left(\frac{\text{log(}0.5\text{)}{B}_t}{q\text{ max(}B\text{)}}\right).\end{array}$$

(6)

We estimated the parameters of the Deriso model with and without the depensation function and compared model estimates of recruitment with the observed recruitment using AICc. We used a lognormal likelihood function

$$\begin{array}{c}L(a,b,g;R,B)=\prod _t\frac{1}{R_t\sqrt{2\pi {\sigma }^2}}\\ \exp \left(-\frac{1}{2{\sigma }^2}\left(\log (R_t)-\log \left(\frac{a{B}_t{D}_t}{{(1+b{B}_t{D}_t)}^g}\right)\right)\right)\end{array}.$$

(7)

Because the depensation function [Equation (4)] is new, we also did the same analysis using the depensatory function used by Myers et al. (1995) shown in Equation (8).

$$R_t=\frac{a{B_t}^{\delta }}{1+b{B_t}^{\delta }}.$$

(8)

In this equation, δ is the depensation parameter, and when δ > 1 the population exhibits depensation. The parameters a and b are the same as those in the Deriso model when g = 1.

Models of depensatory surplus production

Surplus production is the change in biomass that would occur in the absence of harvesting and results from the combined effects of recruitment, somatic growth, and natural mortality. A number of papers have explored surplus production patterns in fish stocks (Hilborn, 2001; Jacobson et al., 2001; Walters et al., 2008).

Surplus production is defined as the net change in biomass, plus harvest.

$$S_t=B_{t+1}-B_t+C_t.$$

(9)

Where S_tis the surplus production over year t, B_t the stock total biomass at time t, and C_t the catch removed between times t and t + 1.

A Pella–Tomlinson surplus production model (Pella and Tomlinson, 1969) was fitted to the data with and without depensation. The Pella–Tomlinson model was chosen rather than the more well-known Schaefer (logistic) model to allow more flexibility in the shape of the production function and because a recent meta-analysis has determined that the shape of the productivity vs. biomass relationship is variable (Thorson et al., 2012). The Pella–Tomlinson model can be written as:

$$\begin{array}{c}{S}_t=\gamma m\left(\frac{B_t}{B_{\mathrm{\infty }}}\right)-\gamma m{\left(\frac{B_t}{B_{\mathrm{\infty }}}\right)}^n\\ \gamma =\frac{n^{n/n-1}}{n-1}\end{array}.$$

(10)

Where S_t is the predicted surplus production over year t, B_∞ the carrying capacity, m the maximum sustainable yield, and n a parameter that determines the shape of the production function. To add depensation, we multiplied the predicted surplus production times the depensatory Equation (4). Because surplus production can take on both positive and negative values, we used a normal likelihood.

Hierarchical Bayesian analysis of depensation

We treated the depensation parameter for each stock s, q_s, as random variables distributed about a population mean and conducted a Bayesian analysis which provides the posterior predictive distribution of q for any new stock s that can be used as a prior distribution.

Using the Deriso model with depensation from Equation (6), we have parameters a, b, g, q, and σ for each individual stock. The key to the hierarchical model formulation for the depensation parameter q_i for any stock is that we assume that it is drawn from a beta distribution that has a two parameters α and β.

The recruits in year t for stock s are given by

$$\begin{array}{c}{R}_{ts}=\frac{a_s{B}_{ts}{D}_{ts}}{{(1+b_s{B}_{ts}{D}_{ts})}^g}\text{exp(}{ϵ}_{ts}\text{)}\\ \\ D_{ts}=1-\text{exp}\left(\frac{\text{log(}0.5\text{)}{B}_{ts}}{q_s\text{ max(}{B}_s\text{)}}\right),\end{array}$$

(11)

where for each stock s, the ε_ts are independently and identically normally distributed with a zero mean and a variance ${\sigma }_s^2$⁠.

$$q_i\sim Beta(\alpha ,\beta ).$$

(12)

The parameters a_s, b_s, g_s σ_s, α, and β were assumed to be uniformly distributed with a broad uniform prior (Table 2). The posterior distributions on all the parameters were found using JAGS software, with 100 000 MCMC runs. We used the following diagnostics: the Geweke diagnostic, which compares the first 10% and last 50% of the Markov chain to see if they come from the same distribution; the Raftery and Lewis diagnostic, which calculates the number of iterations and burn-in necessary to generate accurate results with 95% probability; the Heidelberg and Welch diagnostic, which tests the null hypothesis of the Markov chain being from a stationary distribution. We also examined the Gelman­–Rubin diagnostic and it indicated no problems with convergence. Autocorrelation and trace plots (for mixing) were also evaluated.

Table 2.

Priors used for Bayesian analysis of spawner–recruit data.

Parameter	Prior
as	Uniform 0–10 times the maxR/maxB
bs	Uniform 0–10 times max (R)
Gs	Uniform 0–100
σs	Uniform 0–5
α	Uniform 0–10
β	Uniform 0–100

We conducted the Bayesian analysis separately for four taxonomic groups, invertebrates, pelagics, sharks, rays, and skates, and demersal species. For the recruitment analysis, the number of stocks in each group was 5, 33, 2, and 73, respectively. For the surplus production analysis, the number of stocks was 16, 28, 1, and 66.

Data

The data are gathered from the RAM Legacy Stock Assessment Database version 1.0 (Ricard et al., 2012). Only stocks with time-series for biomass (either total biomass above a specific age as presented in the assessment or spawning-stock biomass) and recruits are included in the analysis. For the spawner–recruit analysis, we are using total biomass as a surrogate for spawning stock when spawning stock size was not available. For the surplus production analysis, one could argue that total biomass would be more appropriate, since somatic growth and natural mortality affect total biomass. Stocks are also excluded which include data that clearly exhibit that the recruitment or surplus production were deterministic model output. Deterministic model output was indicated if either the surplus production or recruitment followed a smooth curve when plotted against biomass. For each stock, years are only included if present in both the biomass and recruits time-series. This results in data compiled from 242 stocks. Of these 242 stocks, we eliminated any stocks that were not observed to have fallen below 20% of maximum observed biomass, leaving 111 stocks for surplus production analysis and 113 stocks for spawner recruit analysis.

Empirical rates of change

The key question is how likely are stocks to rebuild as a function of the stock abundance and the exploitation rate. The logistic growth model is the simplest formulation of population dynamics that addresses this, and the well-known result (Hilborn and Walters, 1992) is that if we plot the rate of increase with biomass on the x-axis and the exploitation rate on the y-axis, there is an isocline going from (0,r/2) to (k/2,0) that represents the division between the parameter space where the stock increases (below and to the left of this line) and the region where stock size decreases (above and to the right). The axis can be transformed into dimensions B/BMSY and U/UMSY in what is now commonly referred to as a Kobe plot (Worm et al., 2009, Figure 3b).

In the RAM Legacy database, we currently have 358 stocks with estimates of B/BMSY and U/UMSY. For each year for each stock, we can locate it in the space B/BMSY and U/UMSY and determine the direction change. We calculate the proportion of stocks that increased when in that region of B/BMSY and U/UMSY space and plot these probabilities. If stocks show strictly compensatory rates of change then we would expect that stocks would be as likely to increase as decrease around the isoclines connecting (0,UMSY) and (BMSY,0) and be more likely to increase the more the stock was in the lower left hand corner, and more likely to decrease towards the upper right hand corner. If depensation was present then the probability of increase would decline as biomass decreased.

Results

Depensation in stock and recruitment

Maximum likelihood

Of the 113 stocks we analysed, only 4 showed an improved AICc when depensation was added to the Deriso model. The Myers model also showed only four stocks had AICc improvements and an estimated value of δ that is greater than 1. The Deriso model provides quite a bit more flexibility in fitting the data, especially allowing for declining right hand limbs, and this flexibility appears to facilitate improvements in fit when adding the depensation term. Figure 2 shows the fitted spawner recruit data for the Deriso model for the four stocks that showed AICc improvement. Values of the estimated depensation parameter q for these stocks were 0.04, 0.06, 0.30, and 0.97 with AICc improvements of 5.4, 1.7, 0.68, and 0.08, respectively. The stock with an estimated of q = 0.97 had very broad confidence intervals.

Figure 2.

Stock–recruitment data and fits with and without depensation for the four stocks with improved AICc scores for the depensatory model. x-axis is spawning stock and y-axis recruitment. The dashed line is the model fit without depensation and the solid line with depensation.

The three stocks the Myers model in our analysis identified as having significantly improved AICc by adding depensation were North Sea Herring, Greenland Halibut 23KLMNO and NAFO Cod 3NO. The North Sea Herring and NAFO Cod 3NO were identified as the first and third biggest improvement in AICc using the Deriso model, but the Greenland Halibut 23KLMNO, did not show a significant depensation effect with the Deriso model. The Deriso model without depensation fits the data for this stock much better because it has a declining right hand limb that the Myers model cannot capture.

Bayesian analysis

The Bayesian analysis also showed little evidence for depensation. The posterior predictive distributions of the depensation parameter (in essence the prior for q for a stock with no data) are shown in Figure 3 for the Deriso model. For all four taxa, the mode is near zero, and for pelagic and demersal with good sample sizes, there is almost no weight for q > 0.04. For invertebrates with a much smaller sample size, there is some density out to 0.10, and for sharks, rays, and skates with only two stocks, really no inference about depensation can be made. For invertebrates, higher qs are less likely but the data are not particularly informative. These results overall reflects the fact that there is little evidence for depensatory dynamics, but because so few stocks have been to very low stock size, there is not enough data to rule out depensation at very low stock sizes. The high values of q estimated for some of the stocks in the maximum likelihood analysis have almost no weight in the Bayesian analysis, because the data for these stocks can be fit almost equally well with very low values of q (the best AICc improvements was 0.08) and many stocks with low q make these low q values much more probable.

Figure 3.

Posterior predictive distribution of the depensation parameter for the Deriso spawner–recruit model for four taxa.

Depensation in surplus production

Maximum likelihood

In the analysis of surplus production using the Pella–Tomlinson model, 8 of 109 stocks showed improvement in the AICc by allowing for depensation. The data and model fits for these eight are shown in Figure 4.

Figure 4.

The surplus production data and fits for the eight stocks with improved AICc scores for the depensatory model. The dashed line is the model without depensation and the solid line is the model with depensation.

Bayesian analysis

The posterior distribution of the depensation parameter (Figure 5) shows a similar pattern to that found for recruitment. The pelagic and demersal stocks show no support for values greater than 0.04 and also have a modal value near zero. The invertebrates with a smaller sample size do not rule out values below 0.1. There is only a single population in the sharks, rays, and skates, and whereas the mode is near q = 0, the data are not particularly informative.

Figure 5.

Posterior predictive distribution of the depensation parameter from the surplus production model for four taxa.

All the Markov chains had matching distributions on the ends in the Geweke diagnostic. The samples were again larger than that suggested in the Raftery and Lewis diagnostic calculation. The Heidelberg and Welch diagnostic indicated that each Markov chain came from a stationary distribution. There was virtually no autocorrelation and the distribution was well mixed.

Sensitivity to recruitment assumptions

One concern about using recruitment and spawning-stock biomass estimated from models as inputs to a statistical analysis is they are both model outputs, and in some cases, there was a spawner–recruit relationship used in estimating the recruitment. To test this, we repeated our Bayesian spawner–recruit analysis using only the datasets from Europe (24 stocks) where the stock assessments do not make any assumption about the spawner–recruit relationship. We found no significant difference in the results for these stocks compared with the entire sample.

Empirical rates of change

Figure 6 shows the probability of stock increase (or decrease) as a function of relative stock size and fishing pressure. The results look remarkably like that expected from the logistic growth model. When stocks are down and to the left of the equilibrium isocline, populations will increase more often than decrease and as the abundance goes lower, and the fishing pressure is decreased, the probability of increase becomes near 1. For example, the lower left hand corner represents U/UMSY < 0.2 and B/BMSY < 0.2, and there were 84 occurrences of stocks in that region and 82% of the time the stock increased.

Figure 6.

Empirical probability of increase as a function of fishing pressure (y-axis) and relative stock size (x-axis).

Discussion

Depensation

These results largely confirm the earlier work of Myers et al. (1995), Liermann and Hilborn (1997), and Keith and Hutchings (2012) that there is little evidence for depensation in stock and recruitment, and this is supported by the surplus production analysis where again little support for depensation was found. Although some stocks do have AICc improvement with the depensatory models, some misclassification due to chance alone would be expected. If we were to cast the analysis in terms of hypothesis testing using the likelihood ratio test, then at the p > 0.05 level of significance, two stock–recruit datasets and four surplus production datasets show a significant depensation parameter. Thus, we observed fewer stocks with significant depensation than we would expect by chance alone.

The frequency with which stocks are at low abundance is important for both the ability to detect depensation and the concern about depensation. Neubauer et al. (2013) suggest that the amount of time stocks are at low stock size is also important. If stocks rarely go to low stock size, we would be unlikely to detect depensation, but we would also be less concerned about depensatory dynamics. Within the range of spawning stock or total biomass below 20% of the maximum, there was a reasonably uniform distribution of the lowest level observed, and for both recruitment and surplus production, we found some occurrences of stocks as low as 1% of the highest abundance.

Empirical rate of increase analysis

The analysis of rates of increase as a function of fishing pressure and stock size provides the same overall picture. When stocks are at low abundance and when they are exploited at a lower rate, they are more likely to increase. At low fishing pressure and low abundance, stocks almost always increased in abundance just as predicted by non-depensatory models.

Low-density dynamics

Concern about low-density dynamics is imbedded in the concept of stock collapse. Francis and Shotton (1997) defined collapse as reduction to low abundance and failure to recover despite reductions in fishing pressure. Worm et al. (2009) and Neubauer et al. (2013) defined collapse as <0.2 BMSY. Many fisheries agencies have management objectives are to maintain stocks at or above the level that produces maximum sustainable yield, which is generally in the realm of 26–46% of the unfished stock size (Thorson et al., 2012). Thus, management policies that are successful in maintaining this level of abundance are unlikely to push stocks into a realm where depensation might occur. Our work suggests that falling below a depensatory population threshold should not be a major concern for management of stocks, but there are many reasons low population sizes should be avoided including lost potential yield, economic profitability, and employment.

Limitations

The ability to detect depensation depends to a great extent on the amount of data taken from low abundance populations. There certainly are many more stocks that have been observed at low abundance available to us now than there were in the analyses done in the 1990s, but still our sample is small. There is also the concern that populations that might show strong depensation were fished to commercial extinction and do not appear in our dataset. Stocks need to have been commercially important in the last 30 years to have assessments available to us.

We chose to use only stocks that had been depleted at some point to below 20% of their maximum observed stock size. In an AICc analysis, there is a trade-off between the number of stocks and the power of the test, choosing a lower threshold such as 10% would have made the probability of detecting depensation if it is present higher but with a smaller sample size. In the Bayesian analysis, the cut-off level is unimportant, and stocks that were never depleted have no influence on the final result.

The model we used was derived by modelling the dynamics of fertilization assuming random encounters between males and females, and it provides one functional form of a model that shows depensation but certainly other forms are possible as evidenced by the forms used in Myers et al. (1995) and Liermann and Hilborn (1997). Given that all three depensatory models lead to similar conclusions, we believe that our conclusions are robust to the specific model form assumed, but it would be useful to model a broad range of depensatory mechanisms from first principles and see if the inference made differs.

The structure of some of the stock assessments used in this analysis may influence recruitment estimates and thus impact meta-analyses such as this one. The original stock–recruitment database developed by Myers et al. (1994) was composed largely of outputs from backward projection models (virtual population analysis) in which no stock–recruitment relationship is specified. In contrast, more than half of the assessments in the RAM Legacy Stock Assessment Database that include recruitment estimates come from forward-projection statistical catch-at-age models (Ricard et al., 2012). Such models typically specify the form of the stock–recruitment relationship, but provide substantial flexibility for annual recruitment to differ from the prediction. Nevertheless, recruitments from such models are not free of influence from the assumed model and this influence must be considered when interpreting the results of meta-analysis. The fact that recruitments from the RAM Legacy Database are generally better explained by a regime shift model than a stock–recruitment relationship (Vert-pre, 2013) suggests that the influence of the specified stock–recruitment relationship is minimal, i.e. the estimated recruitments are largely independent of the model specified in the stock assessment.

Conclusions

Within the range of stock sizes for which we have data, there is little evidence that depensatory processes operate. Few populations for which we have data have been pushed to very low levels (1% of unfished biomass) and we cannot rule out depensation at such stock sizes. We have a very small sample of sharks, rays, and skates and using those data alone we really can say nothing about depensation in those taxa. If we believe that their dynamics are closely related to the taxa with large sample sizes then depensation would likely be of little concern, but their life histories are quite different and they may not have similar dynamics to pelagic and demersal fish.

All evidence points to fishing mortality being the key to stock recovery, when fishing pressure is reduced stocks almost always increase in abundance. However, all stakeholders must recognize that environmental changes can have major impacts on stock recovery and expected time and amount of rebuilding.

Acknowledgements

This research was supported by the National Science Foundation and NOAA through the CAMEO (grant numbers 1041570 and 1041678), the Walton Family foundation, and the University of Washington. We owe a large debt to the many people who have contributed to the RAM Legacy Database.

Author notes

Handling editor: Howard Browman