# Changing recruitment capacity in global fish stocks

See allHide authors and affiliations

Edited by James A. Estes, University of California, Santa Cruz, CA, and approved November 3, 2015 (received for review March 7, 2015)

### This article has a Letter. Please see:

- Relationship between Research Article and Letter - March 08, 2016

## Significance

Marine fish stocks play an important role in marine ecosystems and provide a source of protein for billions of people worldwide. Recent environmental changes have affected the distribution of many stocks, but it is yet unclear whether their productivity is affected as well. We show that recruitment capacity (the ability of stocks to produce surviving offspring) has been significantly altered by both environmental changes and biological changes brought about by overfishing. In total, these effects have reduced recruitment capacity by 3% of the historical maximum per decade, on average. This paper helps us to understand and track previously unrecognized changes in fish stock productivity during the early stages of their life cycle.

## Abstract

Marine fish and invertebrates are shifting their regional and global distributions in response to climate change, but it is unclear whether their productivity is being affected as well. Here we tested for time-varying trends in biological productivity parameters across 262 fish stocks of 127 species in 39 large marine ecosystems and high-seas areas (hereafter LMEs). This global meta-analysis revealed widespread changes in the relationship between spawning stock size and the production of juvenile offspring (recruitment), suggesting fundamental biological change in fish stock productivity at early life stages. Across regions, we estimate that average recruitment capacity has declined at a rate approximately equal to 3% of the historical maximum per decade. However, we observed large variability among stocks and regions; for example, highly negative trends in the North Atlantic contrast with more neutral patterns in the North Pacific. The extent of biological change in each LME was significantly related to observed changes in phytoplankton chlorophyll concentration and the intensity of historical overfishing in that ecosystem. We conclude that both environmental changes and chronic overfishing have already affected the productive capacity of many stocks at the recruitment stage of the life cycle. These results provide a baseline for ecosystem-based fisheries management and may help adjust expectations for future food production from the oceans.

Human well-being is closely linked with the productivity of marine fisheries, which provide a significant source of protein for more than half of the world’s population (1). However, ongoing environmental and biological changes may impact productivity through a variety of mechanisms, including larger habitat areas for temperate species (2), altered body sizes (3), food availability (4), and increased exposure to oxygen-depleted and acidic waters (5). Recent research has documented marked changes in the distributional patterns of marine species that are consistent with climate forcing (6, 7). However, the net effect of these changes on global fish stock productivity is not clearly understood. In particular, documented environmental changes (4, 8, 9) and the long-term consequences of overfishing (10, 11) all impose relevant but poorly constrained effects. Here we help address this issue by evaluating the evidence for empirical trends in the relation between the size of the reproductively mature population (or “spawning stock”) and the annual production of juvenile offspring (“recruits”) using a recently synthesized global database of stock-recruit time series (12). We then test the relation between empirical recruitment trends and regional environmental variables associated with temperature, phytoplankton abundance, and historical overfishing.

Recruitment is modeled by relating the size of the spawning stock biomass to the annual production of recruits. The magnitude of annual recruitment is highly variable (13), yet it provides the basis for population growth and stock productivity by determining the initial number of fish that may grow, die, or be harvested by the fishery (14) (i.e., total productivity is the combination of recruitment, individual growth, and mortality). As such, the stock-recruit relationship has been identified as “the most important and generally most difficult problem in the biological assessment of fisheries” (14). The simplest commonly used recruitment function is the Ricker model_{MAX} is a biomass-independent measure of maximum recruitment and does not depend on current stock size. This property of the measure is attractive as it allows comparison of both abundant and heavily depleted stocks, but it also means that *SI Appendix*), we adopt it as a simple and comparable metric of fish stock productivity at the recruitment stage of their life cycle.

## Results

When recruitment models are fitted to data (Fig. 1 *A–F*), there is often considerable structure in the residual variation (Fig. 1 *G–I*) that suggests that biological productivity may have changed significantly over time. Trends can be observed as directed declines (Fig. 1*G*), threshold-like dynamics (Fig. 1*H*), or regime shifts (Fig. 1*I*; note that the observed shift coincided with the 1977 reversal of the Pacific Decadal Oscillation) (15). We evaluated evidence for changes in recruitment by performing model selection with respect to static or time-varying biological parameters within the Ricker model (i.e., *Methods*) and estimated changes through time where parameters are indeed found to vary. We then summarized trends in recruitment as follows. For an individual stock, we computed the linear slope of *k*. Groupings are made on basis of individual large marine ecosystems (LMEs) and major taxonomic groups. We adopt the LME definition as a simple, ecologically meaningful (16) and management-relevant (17) way to spatially categorize individual stocks. To relate recruitment trends to the environment, we use multiple regression to model _{MSY}). Environmental variables ΔSST and ΔCHL were computed from quality-controlled, publically available databases consistent with the time window covered by stock assessments within individual LMEs, and B:B_{MSY} was calculated as the mean values across all stocks within each LME. See *Methods* and *SI Methods* for full details.

We found that stock-recruitment data supported time-varying recruitment capacity (R_{MAX}) for 71% (*n* = 186) of stocks according to model selection (Fig. 2). Of these, 69% (*n* = 128) showed negative trends (Fig. 2). For all stocks combined, *P* < 0.001; Fig. 2*D*). However, there was a broad-scale divergence in values between the North Pacific and North Atlantic oceans, with the North Atlantic showing steeper declines. In contrast, the North Pacific showed approximately neutral trends across four LMEs, each with a relatively large number of stocks. Across all LMEs, we estimated that 31 of all 39 LMEs (79%), and 20 of 27 LMEs with more than three assessed stocks (74%), showed negative *B*).

There was significant variation associated with different taxa. Groundfish (bottom-associated species such as flatfishes, Pleuronectiformes, and cod-like Gadiformes) showed the most negative *C*). At the species level, the most negative values were observed for several North Atlantic species such as American plaice (*Hippoglossoides platessoides*), European plaice (*Pleuronectes platessa*), common European sole (*Solea vulgaris*), and Atlantic cod (*Gadus morhua*). In the North Pacific, however, many groundfish species showed opposite patterns, with stocks of rex sole (*Glyptocephalus zachirus*), flathead sole (*Hippoglossoides elassodon*), and arrowtooth flounder (*Atheresthes stomias*) trending positively. Pelagic (open-water) species such as herring (*Clupea harengus, C. pallasii*) and swordfish (*Xiphias gladius*) often showed

In general, we found individual stock-recruit parameters changed in a way that resulted in stronger density-dependent processes and reduced maximum reproductive rates. Of individual stocks with negative ΔR_{MAX}, 71% displayed more negative _{MAX} values correlated strongly (*r* = 0.82; *P* < 0.001), suggesting that the observed trends are robust to stocks having variable time series length. We also found that ΔR_{MAX} was generally independent of the assumed form of density dependence in the stock-recruit model or to whether the model let _{MAX}. Likewise, using an alternative metric of recruitment success (expected recruitment at the median historically observed biomass) we found no major change in the resulting trends (see *SI Methods* for details on these sensitivity analyses).

Importantly, average trends in recruitment capacity in each ecosystem were found to be significantly related to environmental and fishing-related variables (ΔCHL and B:B_{MSY,}) across all LMEs (Fig. 3). Considering all species together (Fig. 3*A*), *A*), which accounted for 38% of the total variance. Again, an interesting contrast emerged when isolating the heavily exploited groundfish (combining orders Pleuronectiformes and Gadiformes; Fig. 3*B*). Here, the history of overfishing emerged as the most important predictor and explained 58% of the total variance. Analysis of the pelagic Perciformes and Clupeiformes revealed a positive effect of ΔCHL and negative effect of ΔSST, but these were marginally insignificant. We also investigated patterns of recruitment variation with respect to maximum body size but found no significant relationships.

## SI Methods

Here we provide additional detailed regarding the analysis. For compact description of the estimation, we define the recruitment parameter vector *t* with variance

### Maximum Likelihood Estimation.

Maximum likelihood estimation (MLE) is based on the normally distributed one-step ahead prediction errors of the filtering algorithm (34), termed the innovations. The innovations for the dynamic regression are given by*c* is a constant. Log L is maximized using standard nonlinear optimization, yielding MLE estimates of

### Model Selection.

To determine whether individual recruitment time series have stationary or nonstationary parameters, we applied model selection using various parameterizations of the matrix *i*) static stock-recruit relationship (all elements of *k* = 1); (*ii*) time-varying maximum reproductive rate, static density dependence [element (1,1) nonzero, all others zero; *k* = 2]; (*iii*) static maximum reproductive rate, time-varying density dependence [element (2,2) nonzero, all other zero; *k* = 2]; and (*iv*) time-varying maximum reproductive rate and density dependence [elements (1,1) and (2,2) nonzero, all others zero; *k* = 3]. Note we did not test the full covariance matrix with nonzero off-diagonal elements due to the inherent statistical dependency between the two stock-recruitment parameters (14). Under each parameterization, the Kalman filter/smoother algorithm yields a likelihood that can be used to compute the BIC. {Note that we adopt the BIC over the more common AIC to be more conservative in model selection in the sense that BIC favors fewer parameters due to its stricter penalty term [i.e., for

### ΔR_{MAX} and ΔR M A X k .

The dynamic linear regression analysis yields time series of recruitment parameters *i* to denote an individual stock). We then summarize the trend using a linear slope based on least squares regression*i,* *t,* and *i*) Sample from the uncertainty distribution of the inferred time-varying recruitment parameters to compute bootstrapped time series of *ii*) Compute the standardized linear slope *iii*) Repeat 500 times. This procedure yields bootstrapped distributions of

To combine the results of the individual trend analyses, we perform a random effects meta-analysis at the taxonomic, regional, and global scales. For any particular grouping, the random-effects model is written*rmeta* (32). All meta-analytic results for LMEs and taxa are given in *SI Appendix*.

### Sensitivity Analyses.

#### Sensitivity of R_{MAX} trends to alternative forms of density dependence.

The Ricker function is only one of several possible parameterizations for the recruitment function. Other common recruitment functions include the Beverton–Holt model, *y* axis, altering the *y* intercept, but leaving the slope unchanged. This transformation resulted in constant

#### Sensitivity of R_{MAX} trends to parameterization of Q.

The BIC model selection algorithm was used to choose the appropriate structure of the covariance matrix **Q** by repeating the model selection two times under a restricted set of possible **Q** parameterizations: (*i*) static stock-recuit vs. time-varying *ii*) static stock-recruit vs. time-varying **Q** was parameterized to represent time-varying

#### Sensitivity to an alternative definition of recruitment capacity.

For some heavily depleted stocks, or those exhibiting extremely weak density dependence, the majority of the historically observed biomass may be below the biomass that produces maximum recruitment; therefore, R_{MAX} may not be a practical metric for management. To test the robustness of trends to the specific definition of recruitment capacity, we repeated the trend analysis with an alternative definition; specifically, recruitment evaluated at the median observed biomass, _{MAX} at the regional and global scales, with correlations of 0.85 and 0.92, respectively. Magnitudes of _{MAX} but were generally estimated with smaller bootstrap resampled variances so resulting in very similar *P* values at the global and regional scales.

#### Analysis over the years 1980:2000.

To test the sensitivity of our results to variable time series length, we repeated the time-varying recruitment analysis by subsetting time series for the period 1980:2000 and excluded stocks that do not cover this period. This subset excluded six stocks from the analysis. The scatterplot between

## Discussion

Taken together, these results provide empirical context for understanding contemporary changes in the productivity of exploited marine fish stocks. To date, future forecasts of fisheries productivity have varied in their predictions; for example, the productivity of temperate species has been projected to increase 30–40% based on expansion of fish habitat and increased primary productivity (2), whereas models of individual fish metabolism predict shrinking body sizes with warming oceans (3) that could affect fecundity and productivity. A recent global study of fisheries time series demonstrated that the relationship between adult biomass and total yield can be highly nonstationary (18), but the forcing of such changes has remained unclear. Here we focused our empirical analysis on stock recruitment dynamics and related observed nonstationary patterns to changes in plankton abundance and the history of overfishing. The observed changes in productivity at the recruitment stage of the life cycle may provide a partial explanation for nonstationary patterns observed in fisheries yield for the affected stocks.

We caution that these trends in recruitment biology represent broad-scale spatial and temporal patterns when averaging over many stocks and regions. These patterns should be combined with other model-based forecasts that weigh factors related to habitat quantity and quality to more fully determine expected change in biomass distribution and the productivity of individual stocks. We further note that the drivers of recruitment capacity identified here likely vary in importance among stocks and regions. Bottom-up changes in plankton concentration and top-down effects of overfishing are all known to affect recruitment in complex ways, including effects at both the adult (e.g., maternal effects on recruitment) (19) and larval stages (e.g., food availability) (20). Our results, however, make neither assumptions nor inferences regarding stock-specific mechanisms. Finally, correlations in recruitment may also be important for inferring long-term trends and patterns of shared responses to environmental changes and fishing at the regional scale. We emphasize that a more detailed hierarchical approach that accounts for recruitment correlations (21, 22) and species interactions is needed to fully resolve regional patterns and drivers and thus provide direct management guidance for individual stocks within individual LMEs.

At larger scales, the apparent divergence in productivity among the North Pacific and North Atlantic provides an interesting contrast, possibly linked to divergent ecological histories. The North Pacific experienced a large oceanographic regime shift in the 1970s (15), which resulted in relatively flat long-term environmental trends (Fig. 4). Observed patterns suggest that recruitment capacity may have tracked this variability (Fig. 1*I*), resulting in small *Deepwater Horizon* (24) spill in 2010]. It is also important to note that the database is most representative of North American and European stocks due to the relative scarcity of stock assessments in tropical oceanic regions of the world (Fig. 4*C*) (12, 25) where earth system models predict that plankton productivity will decline more strongly than in the coastal and temperate regions that dominate the stock recruitment database (26). This historical bias in spatial coverage limits our understanding of global fish populations as a whole (25).

In addition to impacting the productivity of marine fish stocks, observed changes in recruitment parameters may also have consequences for the stability of populations. Recent theoretical work has linked observed patterns of population stability (27) to changes in stock recruitment parameters (28) due to age-selective fishing. It was hypothesized that population stability has decreased in stocks due to increases in the mean and variance of the maximum reproductive rate

In summary, empirically estimated trends in recruitment capacity (Figs. 1 and 2) provide evidence for environmental- and fishing-related changes in the productivity of marine fish stocks (Fig. 3). Although far from uniform at the stock level, observed trends were significantly related to ongoing environmental and biological change at the ecosystem scale, specifically changes in phytoplankton biomass and the history of stock biomass depletion (Fig. 4 *B* and *C*). The reality of time-varying biological parameters requires managers to revisit the common assumptions of fixed maximum sustainable yields (30) and emphasizes the need for ecosystem-based management strategies that investigate and account for observed environmental and fishing-related impacts on the long-term productive capacity of fish stocks. Such strategies are enabled by the methods presented here, in that the complex effects of environmental changes can be tracked within a reasonably simple assessment framework. Accounting for such changes is a prerequisite for the successful rebuilding and sustainable harvesting of fisheries resources in a rapidly changing environment.

## Methods

### The RAM Legacy Stock Assessment Database.

All stock recruitment data were extracted from the RAM Legacy Stock Assessment Database (12), which is a global, quality-controlled database, available publicly at ramlegacy.org/. Stock assessments provide estimates of both spawning stock biomass (kilograms) and recruitment (no. individuals). We analyzed 262 of the ∼420 time series available in the database based on (*i*) whether a recruitment relationship was already assumed in generating the stock assessment estimates (12) and (*ii*) whether the spawning stock biomass and recruitment time series were estimated directly, as opposed to indirect proxies such as spawner egg abundance. All series were then normalized to unit variance for easy comparison across stocks and regions. A list of species used in the analysis, along with their designated LME, can be found in Table S1. Frequency histograms of the start and end dates of the stock recruitment time series are shown in Fig. S1 and tabulated in Table S2.

### Nonstationary Recruitment Model.

The Ricker model can be linearized by reexpressing recruitment as log survival

This model can be fitted to data as a linear regression. To model nonstationary recruitment relationships, we let the recruitment parameters vary in time (21, 31) by specifying the following linear Gaussian state space model

where the recruitment parameters are treated as dynamic latent states. Note that **Q**, which are estimated by the method of maximum likelihood. Details of the estimation can be found in *SI Methods*, including the model selection algorithm based on the Bayesian information criterion (BIC). Model selection was used to determine whether variance parameters should be zero or nonzero, thus determining whether the data support static or time-varying recruitment parameters (*SI Methods*). For stocks with at least one time-varying parameter, the trend in *SI Appendix*. Model selection results are also given in Table S3.

### Meta-Analysis.

The nonstationary recruitment analysis and subsequent trend analysis were applied to each stock individually, and then regional and taxonomic patterns were summarized using a random effects meta-analysis model. The random effects model is written

where *k*, *SI Methods*). The meta-analysis model was implemented in the R package *rmeta* (32). All meta-analytic results for LMEs and taxa are given in the *SI Appendix*, which gives group-specific slopes and contributions from individual stocks. Three sensitivity analyses were also performed and are documented in *SI Methods*. These analyses included robustness tests against (*i*) alternative forms of density dependence (Figs. S2–S4); (*ii*) BIC model selection algorithm (Fig. S5); (*iii*) the choice of alternative metrics of recruitment success (Fig. S6); and (*iv*) the impact of variable time series length (Fig. S7).

### Global Scale Correlates of Recruitment.

To correlate recruitment trends to environmental change and overfishing intensity, we fit multiple regression models of the form

where _{MSY} is an index of historical overfishing, representing the mean historical ratio of annual biomass to target biomass levels as extracted from the stock assessments (12),

The multiple regression analysis was fit three times on three sets of species. The first included all species in each LME, and two more subsets (within LMEs) were made on the basis of taxonomic order. One taxonomic grouping included Gadiformes and Pleuronectiformes (generally bottom-associated species) and other included Clupeiformes and Perciformes (pelagic, open-water species). These orders do not occur in all LMEs; therefore, the regression analysis on the subsets included fewer data points. The Gadiformes and Pleuronectiformes occurred in 21 LMEs and the Perciformes and Clupeiformes occurred in 23 LMEs.

## Acknowledgments

We thank all the stock assessment scientists whose dedicated work made this research possible. We also thank C. Minto, R. Mohn, D. Ricard, J. Hutchings, J. Mills-Flemming, and S. Ambrose for discussions; D. Hively, R. Hilborn, and D. Boyce for data; S. Ambrose for data processing; and four anonymous reviewers whose suggestions strengthened the paper. We gratefully acknowledge the Natural Sciences and Engineering Research Council of Canada and the Sobey Fund for Oceans for funding.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: gbritten{at}uci.edu.

Author contributions: G.L.B., M.D., and B.W. designed research; G.L.B. performed research; G.L.B. analyzed data; and G.L.B., M.D., and B.W. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1504709112/-/DCSupplemental.

## References

- ↵.
- FAO

- ↵
- ↵
- ↵
- ↵
- ↵.
- Pinsky ML,
- Worm B,
- Fogarty MJ,
- Sarmiento JL,
- Levin SA

- ↵
- ↵
- ↵
- ↵.
- Conover DO,
- Munch SB

- ↵
- ↵
- ↵
- ↵.
- Hilborn R,
- Walters CJ

- ↵
- ↵.
- Sherman K

- ↵
- ↵.
- Vert-pre KA,
- Amoroso RO,
- Jensen OP,
- Hilborn R

- ↵
- ↵
- ↵
- ↵.
- Thorson JT,
- Stewart IJ,
- Taylor IG,
- Punt AE

- ↵.
- Worm B, et al.

- ↵
- ↵.
- Costello C, et al.

- ↵
- ↵
- ↵.
- Shelton AO,
- Mangel M

- ↵
- ↵
- ↵.
- Peterman RM,
- Pyper BJ,
- Macgregor BW

- ↵.
- Lumley T

- ↵
- ↵.
- Durbin J,
- Koopman SJ

*Time Series Analysis by State Space Methods*(Oxford Univ Press, Oxford) - ↵
- ↵.
- Myers RA

- ↵.
- Myers RA,
- Barrowman NJ

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Ecology

- Social Sciences
- Sustainability Science

### See related content:

- Ecological income inequality in the ocean- Mar 08, 2016