## Abstract

The abundance of a species' population in an ecosystem is rarely stationary, often exhibiting large fluctuations over time. Using historical data on marine species, we show that the year-to-year fluctuations of population growth rate obey a well-defined double-exponential (Laplace) distribution. This striking regularity allows us to devise a stochastic model despite seemingly irregular variations in population abundances. The model identifies the effect of reduced growth at low population density as a key factor missed in current approaches of population variability analysis and without which extinction risks are severely underestimated. The model also allows us to separate the effect of demographic stochasticity and show that single-species growth rates are dominantly determined by stochasticity common to all species. This dominance—and the implications it has for interspecies correlations, including co-extinctions—emphasizes the need for ecosystem-level management approaches to reduce the extinction risk of the individual species themselves.

## 1. Introduction

Assessment of extinction risk and biodiversity loss is a central problem in ecology, which has direct implications for ecosystem management practices [1] and policies for the exploitation of natural resources [2]. It is estimated that currently up to 0.1% of all known species go extinct every year, which is over one thousand times above the background extinction rate observed in fossil records [3]. Whether caused by habitat degradation, interspecies competition, climate change, overexploitation or the introduction of exotic species, the majority of such extinction events have not been anticipated. Existing inventories of the global conservation status, such as the IUCN Red List [4], are believed to include only a fraction of all endangered species and, even for those, significant uncertainty remains on their actual extinction risk. The central difficulty is that wild populations generally do not exist in a steady state from which minute deviations or trends could be detected. Instead, they tend to exhibit fluctuations over time [5–9].

Temporal variations in population abundance may be caused by species interactions, environmental changes, migration patterns, intrinsic nonlinearities and/or human exploitation [10–12], and are sometimes sufficiently irregular to be regarded as stochastic. Such irregular time dependence, when combined with unavoidably imperfect sampling of the population, poses considerable challenges for population forecast. Significant previous research has focused on determining the frequency composition (so-called ‘noise colour’) of such fluctuations and their correlations with the environment [13,14]. Despite their immediate implications for sustainable ecosystem exploitation and management, much less understanding has been generated about the factors that influence the growth rate of individual species.

Perhaps not surprisingly, historically there has been a large number of both false positives and false negatives in the assignment of extinction status and extinction risk. For instance, the New Zealand's bird *takahē*, which was considered extinct by the end of the nineteenth century, was rediscovered in the wild 50 years later [15]. This is one of now many known examples of a Lazarus taxon [16], in which a sparse population passed undetected for an extended period of time. The *passenger pigeon* in North America, on the other hand, went from an abundant population to functional extinction in less than 20 years—which then led to actual extinction in the early twentieth century [17]. Taken together, the picture that emerges is one in which the survival of a species appears to depend very subtly on the species' own abundance [18]. This picture is further complicated by the possibility of co-extinctions [19,20], whose actual role beyond directly dependent species (such as predator–prey and parasite–host) remains elusive [21–23]. Compared to the terrestrial case, extinctions of marine species caused by habitat loss or human exploitation have been rare, in part due to stabilizing effects such as changes in fish catchability at low population densities. Yet, the pace of marine defaunation is likely to accelerate dramatically as the strength and scope of human impact on marine ecosystems grow [24]. This underscores the need for a quantitative understanding of the growth dynamics of marine populations, including extinction risks.

## 2. Results

Here, we examine the dependence of growth rate on population abundance and stochastic factors, both when only the temporal variability of individual species is observed and when the dynamics of the entire ecosystem is taken into consideration. We base our analysis on data of the Large Marine Ecosystems (LME) portion of the Sea Around Us Project (http://www.seaaroundus.org) [25], which is a global-scale database on species abundance in 66 marine ecosystems. For each ecosystem, the data consist of the annual quantities (in tonnes) of its 12 most abundant species caught by fisheries over the period 1950–2006. Such landing data are a widely used tool for making inferences about marine populations [26–29], even if such use is occasionally controversial [30], as factors other than population abundance can influence commercial catches. We will show, however, that our results hold true for marine stock assessments, which integrate research surveys and catch data with independent information (such as mortality rates and population size/age structure) to obtain a more accurate estimate of population abundance (electronic supplementary material, §2, figure S6). Therefore, we present our analysis based on the more widely available landing data [25], followed by supplementary validation using stock assessment data.

For each species *i* in a given ecosystem, we assume that the annual abundance in year *t* can be approximated—up to a scaling factor—by the reported landing, denoted by although we will show that our results do not depend critically on this assumption (electronic supplementary material, §1, figures S1–S3). Our object of study is the year-to-year growth rate, To avoid ill-defined log functions associated with zero landings, we add 1 to each data point, so that the minimum value of is 1. Without any further assumptions, the year-to-year change in the species' population abundance can always be written as
2.1where represents the growth rate in year *t* decomposed into the average and standard deviation *σ*^{(i)} of the growth rate (calculated over the entire time series) and a time-dependent factor at time *t.* The term is thus the normalized growth rate fluctuation. For example, equation (2.1) reduces to the classical Ricker model [31] if is a deterministic decreasing linear function of the abundance which, as shown below, does not hold true for the marine ecosystems we consider. In contrast with previous studies, here we will make no *a priori* assumptions on the growth rate, instead deriving its properties directly from the data.

Figure 1 presents empirical properties of for all 66 ecosystems we consider. We first note that while the autocorrelation of the (log) population abundances is significant and close to one for successive years, as expected (figure 1*a*), the autocorrelation of the corresponding growth rates is comparatively low, with a majority of species having autocorrelation less than 0.2 in magnitude. This surprising observation indicates that the time-dependent component of the growth rate, can be regarded as a random variable drawn from an appropriate distribution that is nearly stationary. The data show that, to an excellent approximation, this distribution is given by a double-exponential function
2.2also known as the Laplace distribution (figure 1*b*–*c*). This is itself an important, novel finding, which further allows us to devise a statistical model, as discussed below. We have verified that our addition of 1's to the zero-landing years does not affect this distribution, which incidentally also governs the fluctuations of growth rate for the population abundance of each ecosystem as a whole (figure 1*d*). In addition, we have rigorously confirmed the fit of the individual species' growth fluctuations to the Laplace distribution versus the null hypothesis of a normal distribution, according to standard goodness-of-fit tests such as the Kolmogorov–Smirnov test and Akaike information criterion. These tests show that the former distribution is a significantly more plausible explanation than the latter distribution for the growth rates of a large majority of species (electronic supplementary material, figure S4).

Having established that the normalized growth rate fluctuations obey a stationary stochastic process described by a parameter-free Laplace distribution, we can propose equations (2.1) and (2.2) themselves as a stochastic model of population abundance. This model is expected to be appropriate away from the *floor* abundance However, the growth rate at the floor is much smaller than at any other abundance, with probability 89% of being zero, as shown in the inset of figure 1*c*. Naturally, because abundance is measured in integer units of landing tonnes, only indicates that the population is very low but not necessarily that it is zero. For this reason (and possibly due to migration and sample biases), the apparent ‘extinctions' considered here are local in space and in time. Indeed, the populations generally recover to detectable abundances in the systems under consideration. However, they do so more slowly than predicted by the overall growth rates. This reduced growth rate at low population density is consistent with the Allee effect (or depensatory population dynamics), which is a scenario previously observed for a number of species in diverse ecosystems [32]. Further analysis would be needed to conclusively address that particular issue in this case, which falls outside the scope of this paper. It should be noted, nevertheless, that the apparent reduced growth rate at low population abundances is counter to most existing models, including the classical Ricker model [33]. A recent study of exploited marine species stocks included in the RAM Legacy Stock Assessment Database [34] indicates that this missing element is in fact the probable explanation for the slow recovery of depleted stocks when compared with predictions from models commonly used in fisheries management [35].

We incorporate the floor effect into our model by taking the probability distribution of to be
2.3and *F*_{2}(*ξ*) = 0 otherwise. The parameter *p*^{(i)} is a measure of the recovery probability once the species reaches the floor abundance, *δ* is the Dirac delta function, and is defined so that when (for simplicity in equation (2.3), we used the approximation *ξ*^{(i)} ≥ 0 instead of the exact condition , since —and hence —is typically close to zero). Here, the average growth rate and standard deviation *σ*^{(i)} are estimated for and the recovery probability *p*^{(i)} is estimated for (see section ‘Parameter estimation’). The model defined by equations (2.1)–(2.3) offers projections on future population abundance based on past abundance, which in turn can be used to assess risk of *pseudo-extinctions*, defined as crossings below a given fraction *θ* of the historical maximum population.

Figure 2 validates our model against empirical data (see section ‘Model predictions and validation’). It shows, in particular, that the floor effect is crucial for the excellent agreement found between the pseudo-extinction predicted and the ones actually observed over the same period (figure 2*a*). We have tested that the heightened agreement with the empirical data holds regardless of the exact value used for the floor abundance (which also represents the limit of resolution in the data), insofar as it is not too large (figure 3). This is significant because current approaches for population variability analysis, including specialized commercial software, generally do not account for the exceptional case of growth at very low population densities (even though minimal viable population sizes are often assumed [36]). Neglecting the floor effect not only mis-predicts pseudo-extinctions observed in the past, but also tends to severely underestimate the risks of future pseudo-extinctions (figure 2*b*). For example, the popular deepwater redfish (*Sebastes mentella*) in Baffin Bay appears to have recovered from very low population abundances, but historical data indicate that when the abundance of this species approaches zero it has a high probability of remaining extremely low for extended periods (figure 2*c*). Neglecting this reduced growth leads to a significant underestimation of the pseudo-extinction risk. This underestimation is even more pronounced for a number of other species, such as the flatfishes (*Pleuronectiformes* spp.) in the Antarctic (figure 2*d*). Interestingly, for some species neglecting the floor effect may actually lead to an overestimation of the pseudo-extinction risk. For example, the blue mussel (*Mytilus galloprovincialis*) in the Iberian Coast reached the floor abundance in the year 2004, but started recovering immediately (figure 2*e*). Because the species exhibited positive growth the only time it reached the floor, predictions that take this into account naturally lead to an estimate of pseudo-extinction risk that is smaller and more reliable than predictions that do not.

Central to our analysis is the fluctuation of the growth rate, modelled as a stochastic term But what is the relation between the stochasticity of different species *i* in the same ecosystem? This question can be addressed by making the ansatz that are not independent but are instead drawn from an appropriate *joint* distribution in which the (marginal) distributions for the individual species follow equation (2.3) while having correlation *α*^{2} with those of the other species (see section ‘Estimation of common stochasticity’). We then consider a range of values of *α*, which represents the portion of stochasticity common to all species, and analyse the integrated impact on the population abundance of each ecosystem as a whole. To give all the species comparable weight, we focus on the average over the individual species' population weighted by the inverse of their average abundance (as in figure 1*d*). The occurrence of pseudo-extinctions is systematically underestimated when stochasticity is dominantly species-specific, as illustrated in figure 4*a* for *α* = 0. In fact, the agreement of the model with empirical data is significantly better when stochasticity is taken to be dominantly common to all species (figure 4*b*), and the agreement becomes excellent for *α* = 0.8 (figure 4*a*).

This corroborates the conclusion that variations in growth rate are largely synchronized within each ecosystem. We speculate that this synchronization is partially rooted in external, environmental fluctuations, which are known to play a crucial role in the dynamics of terrestrial populations [11,37]. In marine systems, likely environmental fluctuations include the El Niño Southern Oscillation, the North Atlantic Oscillation and riverine flood pulses. For fishes, the impacts of these fluctuations can be both direct and indirect—such as via externally driven changes in the lower food web or via ecosystem regime shifts [38]. Other potential sources of synchronized growth fluctuations include interspecies interactions as well as correlations with human activity—for example, shifts in overall fishing effort due to weather or changes in regulations. But regardless of the source of the common stochasticity observed here, the direct implication is an increased risk of the otherwise unlikely concurrent collapse of multiple species.

## 3. Discussion

Our findings should be compared with the case studies of the American breeding bird populations [39] and Hinkley Point's fish community [40], for which growth rates have been analysed. For both systems, the aggregated non-normalized growth rates were found to follow power-law distributions, but independent analysis of some of these data has shown that after rescaling (by the species standard deviation) to a variable equivalent to the normalized growth rate fluctuation the distribution becomes normal [41], which corroborates the conclusion that the growth rate distributions of individual species are short-tailed. This is consistent with the previous analysis of 544 long-term time series from the global population dynamics database [42], including both aquatic and terrestrial populations, which demonstrated that the abundances of most species are either lognormally distributed or shorter tailed than lognormally. If one neglects the (significant) year-to-year correlations in abundance, lognormal distributions for individual population abundances imply normal distributions for individual species' growth rates.

The results presented here, on the other hand, show that normalized growth rate fluctuations follow Laplace distributions and this is confirmed to remain true for individual species in the ecosystems we consider (electronic supplementary material, figure S4). Had we not normalized the growth rates to eliminate heterogeneity across species, the resulting aggregated non-normalized growth rate fluctuations would be more fat-tailed than what we observe, without necessarily revealing a simple scaling behaviour (electronic supplementary material, figure S5). Importantly, we have verified that the observed Laplace distributions for the normalized growth rates are not artefacts of the landing data used here as a proxy for population abundance, remaining valid for stock assessment data. This is demonstrated in electronic supplementary material, figure S6 for the RAM Legacy Stock Assessment Database [34], which is the most complete marine stock assessment catalogue available. The estimates of population abundance therein were acquired under controlled settings, using a variety of methodologies designed to avoid systematic biases, and in a variety of ecosystems, which substantiates the conclusion that Laplace statistics underlie the growth of marine populations, in general. Note that different types of growth rate distributions are consistent with the lognormally distributed abundances observed in previous studies; by equation (2.1), the (log) population abundance is (up to a constant) the sum of the growth rates in all previous years. Thus, by the Central Limit Theorem, *any* growth rate distribution (including Laplace) will eventually lead to normal distributions for the log abundances, provided the growth rates are independent and identically distributed with finite variance. However, direct analysis of the marine population abundances in this study reveals that their distributions are no closer to lognormal than to power laws (electronic supplementary material, figure S7). Altogether, our results are significantly different from those suggested by previous studies, and they do not follow from existing ecological models.

Our demonstration that the growth rates of marine species are governed by the Laplace distribution, which has a heavier tail than a normal distribution with the same standard deviation, has important implications for the analysis of empirical data. Because the likelihood of pronounced fluctuations is larger than expected from normal distributions, large short-term fluctuations (over the period of few years) are not necessarily a sign of abnormality as they may well be a natural property of the system. This, combined with the observed stickiness to the floor abundance, poses additional challenges to the identification of abnormal population dynamics. In particular, we have shown that neglecting the floor effect alone already leads to substantial underestimation of pseudo-extinction risks. Finally, since Laplace distributions have been previously identified in the growth of companies [43], our study establishes a new parallel between ecological and socio-economical networks, both of which are characterized by growth and competition in the presence of limited resources. As such, our results may also provide new insights into common stability mechanisms that govern otherwise disparate systems—as recently proposed in ref. [44].

## 4. Methods

### 4.1. Parameter estimation

To estimate the model parameters from a given time series of population abundance, we disregard the initial consecutive 1's, if any, as they often correspond to years for which no data are available. Using the resulting time series {*x _{t}*}, we compute the associated growth rate time series {

*r*} from the definition

_{t}*r*

_{t}= ln

*x*

_{t+1}− ln

*x*

_{t}, where we now omit the superscript (

*i*) for simplicity (a convention also adopted in the figures). The parameter estimation depends on whether or not we consider the floor effect in the model. For the model without floor effect (equations (2.1) and (2.2)), and

*σ*are simply the sample mean and standard deviation of {

*r*}. For the model with floor effect (equations (2.1)–(2.3)), we have and

_{t}*p*= 1 − (|

*I*

_{00}|/|

*I*

_{0}|), where

*I*

_{1}= {

*t*|

*x*

_{t}> 1},

*I*

_{0}= {

*t*|

*x*

_{t}= 1},

*I*

_{00}= {

*t*|

*x*

_{t}= 1 and

*x*

_{t+1}= 1}, and |·| denotes the number of elements in the set.

### 4.2. Model predictions and validation

For each time series {*x _{t}*} of length

*T*, we use its first

*T*

_{0}data points as the training set for parameter estimation and the remaining

*T*−

*T*

_{0}data points for validation. We calibrate our model both without and with floor effect, and simulate both variants for 1000 independent runs. Each run starts with the same initial condition and is computed for a total of

*T*−

*T*

_{0}time steps, with

*x*reset to 1 whenever its estimated value is below 1. For both the emprical and simulated data, we calculate the pseudo-extinction risk as the probability that , where is the pseudo-extinction threshold and The special case of pseudo-extinction risk at the floor level is defined as the fraction of years for which the population abundance is at the floor (i.e.

_{t}*x*= 1). In this study,

_{t}*T*= 57 years and we choose

*T*

_{0}= 45 years.

### 4.3. Estimation of common stochasticity

We investigate the existence of common stochasticity among all species considered in each ecosystem by extending our model as follows. In a given year *t*, for those species that do not remain at the floor abundance (according to equation (2.3)), we draw a vector of growth fluctuations from a multivariate Laplace distribution [45] characterized by (vector) mean 0, (vector) variance 1, and normalized covariance matrix *Γ*. We set the diagonal elements of *Γ* identically to 1 and the off-diagonal elements identically to *α*^{2}, where *α* is a parameter ranging from 0 to 1. The resulting growth fluctuations for the individual species are distributed according to equation (2.2) but the correlation coefficient is *α*^{2} between and for any pair *i* ≠ *j*. As such, the parameter *α* can be interpreted as the proportion of common stochasticity. The value of *α* that best represents the empirical data is determined in figure 4.

## Funding statement

This study was supported by the U.S. National Oceanic and Atmospheric Administration under grant no. NA09NMF4630406.

## Authors' contributions

All authors participated in the design of the research; J.S. and S.P.C. performed data analysis and simulations; J.S., S.P.C. and A.E.M. wrote the manuscript; all authors commented on the manuscript and gave final approval for publication. J.S. and S.P.C. contributed equally to this work.

## Conflict of interests

The authors declare that they have no competing interests.

- Received March 16, 2015.
- Accepted April 17, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.