# Spatial effects on species persistence and implications for biodiversity

^{a}Laboratory of Ecohydrology, École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland;^{b}Department of Physics, University of Padova, 35131 Padua, Italy;^{c}Department of Civil and Environmental Engineering, Princeton University, Princeton, NJ 08544; and^{d}Dipartimento di Ingegneria Marittima Ambientale e Geotecnica, Universitá di Padova, 35131 Padua, Italy

See allHide authors and affiliations

Edited by Richard T. Durrett, Duke University, Durham, NC, and approved January 27, 2011 (received for review November 18, 2010)

## Abstract

Natural ecosystems are characterized by striking diversity of form and functions and yet exhibit deep symmetries emerging across scales of space, time, and organizational complexity. Species-area relationships and species-abundance distributions are examples of emerging patterns irrespective of the details of the underlying ecosystem functions. Here we present empirical and theoretical evidence for a new macroecological pattern related to the distributions of local species persistence times, defined as the time spans between local colonizations and extinctions in a given geographic region. Empirical distributions pertaining to two different taxa, breeding birds and herbaceous plants, analyzed in a framework that accounts for the finiteness of the observational period exhibit power-law scaling limited by a cutoff determined by the rate of emergence of new species. In spite of the differences between taxa and spatial scales of analysis, the scaling exponents are statistically indistinguishable from each other and significantly different from those predicted by existing models. We theoretically investigate how the scaling features depend on the structure of the spatial interaction network and show that the empirical scaling exponents are reproduced once a two-dimensional isotropic texture is used, regardless of the details of the ecological interactions. The framework developed here also allows to link the cutoff time scale with the spatial scale of analysis, and the persistence-time distribution to the species-area relationship. We conclude that the inherent coherence obtained between spatial and temporal macroecological patterns points at a seemingly general feature of the dynamical evolution of ecosystems.

Understanding local extinction processes has gained urgency as the number of threatened species increases throughout the world because of factors such as habitat destruction or climate change (1–5), but a synthesis of theory and empirical evidence accounting for the relevant ecological dynamics is lacking. In this context, we address here the study of persistence times of trophically equivalent co-occurring species in relation to the spatial scale of observation. The persistence time τ of a species within a geographic region is defined as the time incurred between its emergence and its local extinction (see refs. 6 and 7 and Fig. 1). At a local scale, persistence times are largely controlled by ecological processes operating at short time scales (e.g., population dynamics, dispersal, immigration, and contraction/expansion of species geographic ranges) as local extinctions are dynamically balanced by colonizations (8, 9). At a global scale, originations and extinctions are controlled by mechanisms acting on macroevolutionary time scales.

From a theoretical viewpoint, the simplest baseline model for population dynamics is a random walk without drift, according to which the abundance of a species in a geographic region has the same probability of increasing or decreasing by one individual at every time step. According to this scheme, local extinction is equivalent to a random walker’s first passage to zero, and thus the resulting persistence-time distribution has a power-law decay with exponent 3/2 (10).

A more realistic description can be achieved by accounting for basic ecological processes such as birth, death, migration, and speciation in neutral (11–14) mean-field schemes, as follows. Consider a community of *N* individuals belonging to different species. At every time step a randomly selected individual dies and space or resources are freed up for colonization. With probability *ν* the site is taken by an individual of a species not currently present in the system; *ν* is equivalent to a per-birth diversification rate and it accounts for both speciation and immigration from surrounding communities. With the residual probability 1 - *ν* the individual who died is replaced by one offspring of an individual randomly sampled within the community (15, 16). As such the probability of colonization by a species depends solely on its relative abundance in the community. The asymptotic behavior of the resulting persistence-time distribution [i.e., *p*_{τ}(*t*)] exhibits a power-law scaling limited by an exponential cutoff: [1]with exponent *α* = 2 (7). In formula **1**, time is expressed in generation time units (11); i.e., it has been rescaled such that the birth rate is equal to one. Notably, in the mean-field scheme the probability distribution *p*_{τ}(*t*) depends solely on the diversification rate that accounts for speciation and migration processes and imposes a characteristic time scale 1/*ν* for local extinctions. Although per-birth speciation rates are not expected to vary with the spatial scale of analysis, per-birth immigration rates are argued to decrease as the spatial scale increases. In fact, the possible sources of migration (chiefly dependent on the geometrical properties of the boundary and the nature of dispersal processes) are argued to scale sublinearly with the community size (17), which in turn is typically linearly proportional to geographic area (2, 8). As continental scales are approached, migration processes (almost) vanish and the diversification rate ultimately reflects only the speciation rate.

From an empirical viewpoint, species and genera persistence times deducted from fossil record data have been suggested to follow either power-law [with nontrivial exponents in the range 1.5–2 (18–20)] or exponential distributions (19, 21). It has been argued, however, that data quality, in particular for species, precludes a critical assessment (7). Also, local analyses of species persistence over ecological time scales suggest power-law distributions with nontrivial exponents (6).

In what follows we provide evidence for power-law behavior, either empirically or from a broad spectrum of theoretical derivations. Implications on emerging macroecological patterns are examined, with special attention to possible biogeographical effects.

## Empirical Persistence-Time Distributions

We empirically characterize species persistence-time distributions by analyzing two long-term datasets covering very different spatial scales: (*i*) a 41-y survey of North American breeding birds (22) and (*ii*) a 38-y inventory of herbaceous plants from Kansas prairies (23).

The North American Breeding Bird Survey consists of a record of annual abundance of more than 700 species over the 1966–present period along more than 5,000 observational routes. The spatial location of the routes analyzed is shown in Fig. 1. We consider only routes with a latitude less than 50° because density of routes with a long surveyed period drastically decreases above the 50th parallel. Noting that in many regions the survey started only in 1968, we discard the first two years of observations in order to have simultaneous records for all the regions in the system. The spatial extent of the observational routes allows us to analyze species persistence times at different spatial scales. We consider 20 different scales of analysis with linearly increasing values of the square root of the sampled area starting from *A* = 10,000 km^{2} to *A* = 3.8·10^{6} km^{2}. We also analyze the whole system, which corresponds to an area of *A* = 7.8·10^{6} km^{2}. For every scale of analysis *A* we consider several overlapping square cells of area *A* inside the system. A three-dimensional presence–absence matrix *P* is thus built. Each element *p*_{stc} of the matrix is equal to 1 if species *s* is observed during year *t* in at least one of the observational routes comprised in cell *c*; otherwise *p*_{stc} = 0. For every scale of analysis we discard the cells that (*i*) do not have a continuous record for the whole period (41 y) or (*ii*) have more than 5% of their area falling outside the system. For every cell and every species we measure persistence times from presence–absence time series derived from the second dimension of matrix *P*. Persistence time is defined as the length of a contiguous sequence of ones in the time series. For every scale of analysis we consider all the measured persistence times regardless of the species they belong to and the cell where they are measured. The effect of possible imperfect detection of species (24) on measured persistence times has also been investigated (see *SI Text*).

The herbaceous plant dataset (23) comprises a series of 51 quadrats of 1 m^{2} from mixed Kansas grass prairies where all individual plants were mapped every year from 1932 to 1972. In order to meet the data quality standard required for our analysis as discussed above for the breeding bird data, we discard 10 quadrats and the first three years of observations. Because of the limited number of observational plots in the herbaceous plant dataset, we limit our analysis to quadrat spatial scale *A* = 1 m^{2}. Analogously to the previous case, we reconstruct the matrix *P* from presence–absence data for every species, year, and quadrat.

Note that, when dealing with empirical survey data, the effect of the finiteness of the observational time window on the measured species persistence times must be properly taken into account. To this end, we have developed tools to extend the inference of persistence-time distributions for periods longer than the observational window. In particular, we analytically derive, given the persistence-time probability density function, the distribution of two additional variables that can actually be measured from empirical data: (*i*) the persistence times *τ*^{′} of species that emerge and go locally extinct within the observed time window Δ*T*_{w} and (*ii*) the variable *τ*^{′′} that comprises *τ*^{′} and all the portions of species persistence times that are partially seen inside the observational time window but start or/and end outside (Fig. 2*A* and *Materials and Methods*). The finiteness of the time window imposes a cutoff to *p*_{τ′}(*t*). On the contrary *p*_{τ′′}(*t*) has an atom of probability in *t* = Δ*T*_{w} corresponding to the fraction of species that are always present during the observational time. By matching analytical and observational distributions for *p*_{τ′}(*t*) and *p*_{τ′′}(*t*), it is possible to infer the persistence-time distribution *p*_{τ}(*t*). The scaling exponent and the diversification rate for the herbaceous plant persistence-time distributions have been determined with a simultaneous nonlinear fit of observational and analytical *p*_{τ′}(*t*) and *p*_{τ′′}(*t*). Confidence intervals are equal to the standard error of the fit. For breeding birds, we repeat the nonlinear fit for different spatial scales of analysis. The reported scaling exponent and the confidence interval have been obtained by averaging results across spatial scales.

Remarkably, the persistence times of breeding birds at different spatial scales of analysis and of herbaceous plants prove to be best fitted by a power-law distribution with an exponent *α* = 1.83 ± 0.02 and *α* = 1.78 ± 0.08, respectively (Fig. 2 *B* and *C*). It is important to note that both scaling coefficients derived empirically are significantly different from the predictions of the existing baseline models discussed above (the random-walk persistence time yields an exact exponent *α* = 3/2; the mean-field model yields *α* = 2).

## Theoretical Persistence-Time Distributions

In this section we provide a theoretical rationale for the universality of the scaling behavior of persistence-time distributions with respect to the topology of the interactions allowed by the environmental matrix. In particular, we provide evidence on how nontrivial exponents of the type observed empirically can be reproduced by simple theoretical models once dispersal limitation and the actual network of spatial connections are taken into account. We have implemented the neutral game described above in regular one-, two-, and three-dimensional lattices in which every site represents an individual (15, 16). We have also explored the patterns emerging from the application of the model to dendritic structures mimicking riverine ecosystems where dispersal processes and ecological organization are constrained by the network structure. Indeed, many features of riverine ecosystems have been shown to be affected by the connectivity of river networks (25, 26). In particular, river geometry has been studied in relation to extinction risk (27), migration processes (28), persistence of amphibian populations (29), macroinvertebrate dispersal (30), and freshwater fish biodiversity (14, 31). For general calculations of the topological structure and metric properties relevant to dendritic ecological corridors, we employ the features of optimal channel networks (OCNs) (32). They hold fractal characteristics known to closely conform to the scaling of real networks (33). Among the advantages of the use of OCNs, one recalls the possibility to fit one such construct into any assigned domain (e.g., a square, Fig. 3), thus allowing exactly the same size and number of nodes of a two-dimensional lattice to be endowed with altered directionality of connections. To account for limited dispersal effects, we allow only the offsprings of the nearest neighbors of the individual who died to possibly colonize the empty site. In the networked landscape the neighborhood of a site is defined by the closest upstream and downstream sites. Limited dispersal promotes the clumping in space of species, which enhances their coexistence and survival probability (16, 34). Indeed we find that in all the considered landscapes, persistence-time distributions still follow a power-law behavior characterized by smaller, nontrivial scaling exponents (namely *α* = 1.92 for the 3D, *α* = 1.82 for the 2D, *α* = 1.62 for the OCN, and *α* = 1.50 for the 1D landscape, Fig. 3) limited by an exponential cutoff. Remarkably, the exponent obtained via simulation in a two-dimensional landscape (*α* = 1.82 ± 0.01) is close to those found in both breeding birds and herbaceous plants datasets (*α* = 1.83 ± 0.02 and *α* = 1.78 ± 0.08, respectively).

We also study how persistence-time distributions deducted from the theoretical model change with dispersal broader than nearest neighbors (see *SI Text*). As expected, as long as the mean dispersal distance remains small with respect to the system size, the distribution eventually ends up scaling as the one predicted by the nearest-neighbors dispersal. We also relax the neutral assumption by implementing an individual-based competition/survival trade-off model (16). Specifically, species with higher mortality rates are assumed to hold less competitive ability to colonize empty sites (2, 35). It is important to note that the trade-off model also exhibits power-law persistence-time distribution with exponents indeed close to those shown by the neutral model (see *SI Text*). Our theoretical results are thus robust with respect both to changes in the dispersal range and to relaxations of the neutrality assumptions. This confirms our expectation that a power-law distribution for species persistence times is the result of emergent behaviors independent of fine ecological details, thus supporting the neutral assumption that effective interaction strength among species is weak (36) and does not significantly constrain the dynamics of ecosystems. We also note that our results are not seen as a test for the neutrality hypothesis for breeding birds or herbaceous plants dynamics, but rather as tools to reveal emerging universal and macroscopic patterns (37, 38).

## Discussion and Conclusions

In the previous section we established a hierarchy of scaling exponents ranging from the smallest, proper to 1D interactions, to larger values namely for directional (network-like), 2D, and 3D dispersals. We thus suggest that the coherence of the empirical scalings would stem from the two-dimensional isotropic nature of the environmental matrix available to the ecological processes relevant to both breeding birds and herbaceous plants.

We also suggest that species persistence-time distribution, owing to its robustness and scale-invariant character, is a synthetic descriptor of ecosystem dynamics and of biodiversity. In fact, other key macroecological patterns are intimately related to the persistence-time distribution. A first clear example is the direct link with ecosystem diversity, as explained below. In our framework species emerge as a point Poisson process with rate *λ* = *νN* and last for a persistence time *τ*. The mean number of species *S* in the system at a given time is therefore *S* = *λ*〈*τ*〉 (39) where 〈*τ*〉 is the mean persistence time. Therefore, the smaller exponents found, say, for networked environments with respect to two-dimensional ones, imply longer mean persistence and, in turn, higher diversity. This echoes recent results suggesting a higher diversity of freshwater versus marine ray-finned fishes (40, 41).

Another evidence of the effective way in which species persistence-time distribution can characterize ecosystem diversity is the link with the species-area relationship, which characterizes the increase in the observed number of species with increasing sample area. The spatial extent of the breeding bird dataset and the tools developed for the data analysis allow us to study how the persistence-time distribution depends on the spatial scale of analysis (Fig. 4*A*). As expected, while the scaling exponent remains the same, the diversification rate *ν* decreases with the geographic area *A* and is found to closely follow a scaling relation of the type *ν* ∝ *A*^{-β}, with *β* = 0.84 ± 0.01 (Fig. 4*B*), for a wide range of areas. This scaling form of the cutoff time scale 1/*ν* can be related to the species-area relationship. Assuming that the number of individuals scales isometrically with the sampled geographic area (2, 8), i.e., *N* ∝ *A*, and given that 〈*τ*〉 = ∫*tp*_{τ}(*t*)*dt* ∝ *ν*^{α-2} (see *SI Text*) one gets [2]The observational values *β* = 0.84 ± 0.01 and *α* = 1.83 ± 0.02 give an exponent *z* = 0.30 ± 0.02, which is close to the species-area relation measured directly on the data for the same range of areas (*z* = 0.31 ± 0.02, Fig. 4*C*). Conversely, one could have used the observed species-area exponent to infer the scaling properties of the diversification rate.

Finally, from a conservation perspective, a meaningful assessment of species’ local extinction rates is deemed valuable. We propose the distribution of the times to local extinction *τ*_{e} (Fig. 2*A*) as a tool to quantify the dynamical evolution of the species assembly currently observed within a given geographic area. Mathematically, *τ*_{e} is defined as the time to local extinction of a species randomly sampled from the system, regardless of its current abundance. When formula **1** holds for persistence times, the distribution of the times to local extinction *p*_{τe}(*t*) is shown to scale as *p*_{τe}(*t*) ∝ *t*^{1-α}*e*^{-νt} (see *Materials and Methods*). Therefore, not only do the developed theoretical and operational tools allow to infer the scaling behavior of persistence times, but also of the time to local extinction even from relatively short observational windows. Although these patterns cannot provide information about the behavior of a specific species or of a particular patch inside the ecosystem considered (e.g., a biodiversity hot spot), they can effectively describe the overall dynamical evolution of the ecosystem diversity. In particular, the scaling behavior allows us to extrapolate species persistence-time distributions for wide geographic areas, which are hard to estimate, from measures of persistence on smaller areas, which are, on the contrary, more practical and feasible. We thus conclude that the biogeographical characters of species persistence, stemming from the structure of the spatial interaction networks and from local constraints to species emergence rates, add a previously undescribed ingredient to a rich literature bearing major implications for the inventory of life on Earth.

## Materials and Methods

### Inference of the Persistence-Time Distribution from a Finite Observational Period.

The exact derivation of the probability distribution of the variables *τ*^{′} and *τ*^{′′} (Fig. 2*A*) follows. In this theoretical framework, the probability *νdt* of observing a diversification event in a time step *dt* is assumed to be a constant; thus species emergence in the system due to migration or speciation is seen as a uniform point Poisson process with rate *λ* = *νN* (where *N* is total number of individuals in the system and λ has the dimension of the inverse of a generation time). We term *t*_{0} the emergence time of a species in the system, and *T*_{0} and *T*_{f} = *T*_{0} + Δ*T*_{w} the beginning and the end of the observational time window, respectively. A species emerged at time *t*_{0} will be continuously present in a geographic region for its persistence time τ until its local extinction at time *t*_{0} + *τ*.

We first analyze the distribution of *τ*^{′′}, the most complex case. The variable *τ*^{′′} can be expressed as a function of the random variables *τ* and *t*_{0}, which are probabilistically characterized. We can distinguish four different cases (Fig. 2*A*):

the species emerges and goes locally extinct within the time window;

the species emerges during the observations and it is still present at the end of the time window;

the species emerges before the beginning of the observations and goes locally extinct within the time window;

the species is always present for all the duration of the observations;

or, mathematically, We express the probability of observing *τ*^{′′} conditional on a persistence time of duration *τ* as [3]where the operator 〈·〉 is the ensemble average with respect to the random variable *t*_{0}, *δ*(*x*) and Θ(*x*) are the Dirac delta distribution and the Heaviside function, respectively. is the normalization constant.

When comparing analytical and observational distributions, we assume that the system is at stationarity and unaffected by initial conditions; i.e., *T*_{0} is far from the beginning of the process. Mathematically this is obtained taking the limit *T*_{0},*T*_{f} → +∞ with *T*_{f} - *T*_{0} = Δ*T*_{w}. By solving the ensemble averages and by marginalizing with respect to *τ*, Eq. **3** finally takes the form (see *SI Text* for a step-by-step derivation) [4]where simplifies to with being the exceedance cumulative distribution of the persistence-time probability density function.

The variable *τ*^{′} comprises only the first of the four cases listed in Eq. **3**. Thus the probability distribution *p*_{τ′}(*t*) follows directly from the first term of Eq. **4**,where the normalization constant is equal to which completes the derivation.

### Distribution of Times to Local Extinction.

We term *τ*_{e} the time to local extinction of a species randomly sampled among the observed assembly at a certain time *T* (Fig. 2*A*). Analogously to the derivation described above, we can express *τ*_{e} as We then express the probability distribution of the times to local extinction conditioned to a persistence time *τ* as where the constant ensures proper normalization. Solving the ensemble average operators yields Marginalizing over *τ* and considering the system at stationarity (*T* → +∞), we finally obtain [5]where is simply 〈*τ*〉. Eq. **5** allows us to derive the distribution of the times to local extinction given the persistence-time distribution. Particularizing now to the case of persistence-time distributions of the shape *p*_{τ}(*t*) ∝ *t*^{-α}*e*^{-νt}, Eq. **5** translates into *p*_{τe}(*t*) ∝ *t*^{1-α}*e*^{-νt}.

## Acknowledgments

The authors thank Dr. Stephen Hubbell and two other anonymous referees for their helpful comments and suggestions. We thank Miguel Munoz for discussions on the scaling properties of survival probabilities and Marino Gatto and Renato Casagrandi for useful comments. E.B., S.S., L.M., and A.R. gratefully acknowledge the support by the European Research Council (ERC Advanced Grant RINEC-227612) and by the Swiss National Science Foundation (Project 200021_124930/1). I.R.-I. acknowledges the support of the James S. McDonnell Foundation (Grant 220020138).

## Footnotes

^{1}To whom correspondence should be addressed. E-mail: andrea.rinaldo{at}epfl.ch.

Author contributions: E.B., S.S., L.M., A.M., I.R.-I., and A.R. designed research, performed research, analyzed data, and 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.1017274108/-/DCSupplemental.

## References

- ↵
- Diamond J

- ↵
- Brown JH

- ↵
- ↵
- Svenning J.-C,
- Condit R

- ↵
- May RM

- ↵
- ↵
- Pigolotti S,
- Flammini A,
- Marsili M,
- Maritan A

- ↵
- MacArthur RH,
- Wilson EO

- ↵
- Ricklefs R

- ↵
- ↵
- Hubbell S

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Sneppen K,
- Bak P,
- Flyvbjerg H,
- Jensen M

- ↵
- Solé R,
- Bascompte J

- ↵
- Newman M,
- Sibani P

- ↵
- VanValen L

- ↵
- Department of the Interior Geological Survey, Patuxent Wildlife Research Center, Laurel, MD

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Grant EHC,
- Nichols JD,
- Lowe WH,
- Fagan WF

- ↵
- ↵
- Bertuzzo E,
- et al.

- ↵
- Rodriguez-Iturbe I,
- et al.

- ↵
- ↵
- ↵
- ↵
- Volkov I,
- Banavar JR,
- Hubbell SP,
- Maritan A

- ↵
- ↵
- ↵
- Rodriguez-Iturbe I,
- Cox D,
- Isham V

- ↵
- ↵
- Moyle P,
- Chech J

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Ecology

- Physical Sciences
- Applied Mathematics