Extinction debt in local habitats: quantifying the roles of random drift, immigration and emigration

We developed a time-dependent stochastic neutral model for predicting diverse temporal trajectories of biodiversity change in response to ecological disturbance (i.e. habitat destruction) and dispersal dynamic (i.e. emigration and immigration). The model is general and predicts how transition behaviours of extinction may accumulate according to a different combination of random drift, immigration rate, emigration rate and the degree of habitat destruction. We show that immigration, emigration, the areal size of the destroyed habitat and initial species abundance distribution (SAD) can impact the total biodiversity loss in an intact local area. Among these, the SAD plays the most deterministic role, as it directly determines the initial species richness in the local target area. By contrast, immigration was found to slow down total biodiversity loss and can drive the emergence of species credits (i.e. a gain of species) over time. However, the emigration process would increase the extinction risk of species and accelerate biodiversity loss. Finally but notably, we found that a shift in the emigration rate after a habitat destruction event may be a new mechanism to generate species credits.

species loss is a time-delayed process that is usually described as extinction debts [1][2][3][4], many profound ecological problems are still poorly known and remain open for exploration, including: (i) How fast will species go extinct due to abrupt habitat destruction? (ii) To which factors will the velocity of species extinction be correlated? (iii) Can we accurately predict the temporal trajectory of species loss? (iv) How will community patterns (e.g. species abundance distributions) be altered if species loss is inevitable?
Many previous studies attempted to demystify extinction debts using a variety of approaches. For example, several previous studies [5,6] used 'static' (i.e. time-irrelevant) approaches to model delayed extinction patterns caused by habitat destruction. However, such methods are unable to characterize the temporal trajectories of species extinctions, which are clearly time-dependent. In comparison, inspired by the theory of island biogeography, some other papers developed dynamic models (mostly related to species-area relationships) to predict the delayed loss of species richness [7][8][9][10]. However, these methods typically require an estimation of a relaxation parameter, which is usually difficult to fit [7] and has few empirical values to which one can refer [11]. Moreover, other than species richness, community-level delayed consequences (e.g. compositional change: how abundant and rare species will change) cannot be clearly demonstrated using these dynamic models.
Neutral theory and models can be powerful for quantifying delayed consequences of habitat loss on biodiversity. By assuming the functional equivalence of different species subject to stochastic birth, death and dispersal processes, neutral models are able to quantify how individual species and the entire ecological community will respond to habitat destruction [12][13][14]. However, no previous studies have systematically employed and developed neutral models to address the above-mentioned questions.
In this study, we developed a general neutral model that can predict the transitional behaviours of species extinctions by incorporating multiple essential ecological mechanisms: emigration/immigration, random drift and habitat destruction. More importantly, compared to previous studies, our model can explicitly project how community patterns will shift from both spatial and temporal perspectives.

A stochastic dispersal-birth-death model
Here, we use the following stochastic dispersal-birth-death (SBD) model to describe the temporal dynamics of a species' population under random drift, immigration and emigration processes dp 0 ðtÞ dt ¼ p 1 ðtÞ À vp 0 ðtÞ dp 1 ðtÞ dt ¼ vp 0 ðtÞ þ 2p 2 ðtÞ À ð2 À uÞp 1 ðtÞ . . . dp j ðtÞ dt ¼ ðj þ 1Þp jþ1 ðtÞ þ ðj À 1Þð1 À uÞp jÀ1 ðtÞ À jð2 À uÞp j ðtÞ: Note that p j (t) in equation (2.1) is the probability that a given species has abundance j at time t. We treat v as the immigration rate and u as the emigration rate. The immigration rate is necessary for the model; otherwise, no species will exist in the study area, and community dynamics are not possible if the area originally had no organisms in existence, and no individuals could disperse into the area from the outside. A flowchart for showing the transition dynamic of neighbouring abundance states from equation (2.1) is visually illustrated in electronic supplementary material, figure S1A of the appendix. This model only allows the flow of immigration when there are no individuals inside the target area (vp 0 (t)). This assumption seems a bit strong, but there are a variety of reasons for assuming so (details are introduced in the Discussion section). To further demonstrate the usefulness of the proposed model (equation (2.1) as to the evaluation of stochastic dispersal process in population dynamic and extinction, we also introduce and compare a full immigration model [15,16], detail of which is shown in the Additional Methods of the appendix (electronic supplementary material, figure S1B and equation S1). In that full immigration model, immigration can take place as long as the population of species changes in an increasing direction (electronic supplementary material, figure S1B). By defining the probability generating function (PGF) as When there is no immigration effect (i.e. v = 0), the above solution can be recovered to a previously described one [17], and there is an analytical solution for the PGF H(z, t), which has a form as follows [17]: for the specific initial condition p n (0) = δ(n − j ) ( j ≥ 0), where δ(w) = 1 if w = 0; and δ(w) = 0 if w ≠ 0. The probability of observing n individuals at time t can be derived by taking the n-ordered derivatives of the PGF as However, it becomes challenging to analytically solve the PDE (equation (2.2a)) using the traditional method of characteristics when the immigration rate is not zero (i.e. v ≠ 0), since the rightmost term of equation (2.2a) still involves H(0, t). To this end, we use a matrix exponential method [18,19] when conducting numerical analyses. Specifically, we use the following matrix to record the coefficients of the ordinary differential equation (ODE) system in equation (2.1) as

ð2:3Þ
and a column vector to record the time-dependent probability at each abundance state as X m (t) ¼ ( p 0 (t),p 1 (t),p 2 (t), . . . ,p m (t)): For numerical computing purpose, we convert the infinite ODE system equation (2.1) into a finite ODE system with m + 1 states (i.e. from 0 to a maximal abundance value denoted by m; integer m > 1). Therefore, equation (2.1) can be re-formulated as the following finite matrix form Accordingly, we can numerically solve the time-dependent probabilities by using matrix exponential method as With the result from equation (2.5), the time-dependent probability of extinction can be easily extracted from X m (t), i.e. its first element ( p 0 (t)). As a remark, because the original stochastic system in equation (2.1) is infinite, it is necessary to check the convergence of the derived extinction probability p 0 (t) for different matrix sizes (m) under the finite setting.

Total biodiversity loss at the metacommunity level
The probability that a species can still survive at time t given that its initial abundance, j, under the above stochastic SBD process (equation (2.1)) is given by Therefore, without both emigration and immigration effects (u = v = 0), equation (2.6) characterizes the effect of random shift on extinction and survival dynamic [12], which can be derived using equation (2.6) as p n.0 ðtÞ ¼ 1 À (t=ð1 þ tÞ) j . If the emigration rate of a species is not zero, but immigration is absent (i.e. u > 0 and v = 0), equation (2.6) returns to p n.0 ðtÞ ¼ 1 À ð1 À e Àut Þ=ð 1 À ð1 À uÞe Àut Þ ð Þ j .
Finally, if both emigration and immigration rates are not zero, we employ the matrix exponential method (equation (2.5)) to numerically evaluate its temporal patterns. At the metacommunity level, suppose that there are S(0|A) species in the region with size A at the initial time t = 0. Moreover, let S j (0|A) denote the number of species with abundance j at the initial time; then with no emigration or immigration effects (u = 0 and v = 0; under the pure random drift scenario), the expected species number at any time t is given by Therefore, the expected number of individuals at time t is and represents the number of species with abundance j in the region with size A at time t. The expected total biodiversity loss at time t at the metacommunity level, under either random drift or emigration processes, can be estimated as Eðt j AÞ ¼ Sð0 j AÞ À Sðt j AÞ: ð2:8Þ When the immigration rate does not exist (i.e. v = 0), we can express E(t|A) by This time-dependent solution for the total species loss at the metacommunity level is a concave function with respect to time t when u = 0; however, the curvilinear pattern of E(t|A) for the case u > 0 is indeterministic; one can refer to Additional Methods of the appendix in electronic supplementary material for a proof. By contrast, the species richness dynamic model presented in equation (2.4) of Sgardeli et al.'s [20] paper is a convex function in most cases when studying extinction debt [1,7,10,21]. A proof is also presented in electronic supplementary material, appendix. Accordingly, the derived biodiversity loss model in Sgardeli et al.'s [20] paper would become a concave function again when applying equation (2.8). Moreover, by properly assigning model parameters, both equation (2.8) and the model in Sgardeli et al.'s [20] paper will have similar prediction of temporal trajectory of species loss pattern (electronic supplementary material, figure S4). However, we have to mention here, when the immigration rate is not zero, the expected total biodiversity loss in our model (equation (2.8)) has to be solved out numerically.

Total biodiversity loss at the local community level
Suppose in the initial stage, there is a local area of size a that is part of the region (a , A), and the local species abundance distribution (SAD) is given by S j (0|a), which represents the number of species with abundance j at the initial time in local area a. Accordingly, the expected total biodiversity loss at time t at the local level, under random drift, immigration and emigration processes, can be estimated as Eðt j aÞ ¼ Sð0 j aÞ À Sðt j aÞ, ð2:10Þ where S(0 j a) ¼ P 1 j¼1 S j ð0 j aÞ; S(t|a) represents the expected species number of the local area at any time t and is estimated as  To better incorporate instant habitat destruction on the transition behaviour of the species richness with Fisher's explicit statistical background, and accordingly, estimate the expected species loss over time, we assume there are Sð0 j AÞ ¼ 1000 species at the initial time point t = 0 at the regional scale; and the SAD at either the regional or local scale follows an area-based logseries abundance model [22]. To be specific, the area-based model is given as follows.
Suppose a truncated negative binomial model is best fitted to the abundance distribution of all species necessarily present in target region A, and its probability mass function (PMF) is ; N A is a random variate to depict the abundance of each species in the entire region A; k and ω are two positive parameters of the PMF. Using this model with k → 0, the limiting distribution can be derived by 13Þ This model can be further simplified to contain only one unknown parameter, α A , resulting in the following form: Suppose we only sample a local area of size a from the entire region A. To do this, we define the number of individuals of each species observed in the local sampled area a as N a ; then the PMF of N a is given by Detailed derivation of equation (2.15) can be found in our previous paper [22] and is thus omitted here. Therefore, the expected species richness with abundance j in the local ecological community a in the entire region A, given the parameter value α A , is given by As a result, the time-dependent and area-dependent total biodiversity loss of local sample a, before habitat destruction point τ, can be computed as E a ðt j t , tÞ ¼ S a ð0Þ À S a ðtÞ ¼ X j!1 S j ð0 j a,a a Þp 0 ðtÞ: ð2:16Þ Suppose at some time point, τ, a local area of size b of region A is instantly destroyed; then the number of extinct endemic species found in local area b can be estimated by Accordingly, the number of species found in the remaining habitat (A − b) is given by Then, the expected species richness with abundance j in an intact local area, given the instantaneous destruction of another local area b at time τ, is given by j a a ,a A ): Therefore, the total biodiversity loss at local sample a, after the habitat destruction point τ, becomes  figure S2 of the appendix). Moreover, the immigration rate would retard the risk of species extinction, while the emigration rate would accelerate the extinction risk, as demonstrated from the diverse curve shape patterns of different combinations of immigration and emigration rates (electronic supplementary material, figure S3 of the appendix).
In the absence of habitat destruction, the total biodiversity loss dynamic tends to present a concave trend over time ( figure 1). When an immediate habitat destruction event takes place, because of imminent extinction (i.e. loss of endemic species that are only found in the destroyed habitat), the curve is discontinuous (i.e. has an abrupt change) at the time point of habitat destruction ( figure 1).  After that, biodiversity loss continually accumulates and becomes smooth again, but the debt magnitude is always smaller than the original dynamic without habitat destruction (figure 1). The curve shape patterns were actually very similar between the proposed model (equation (2.1)) and the full immigration model (electronic supplementary material, equation S1). High imminent extinction is expected when the areal size of the destroyed local habitat is great ( figure 1). However, when the areal size of the destroyed habitat is higher, the expected total biodiversity loss after habitat destruction is lower ( figure 1). This is because imminent extinction is paid off at the point in time when habitat destruction occurs.
The shape of the curve of the original SAD influences the shape of the curve of the total biodiversity loss dynamic (figure 2). When parameter ω in Fisher's logseries model is lower, the expected number of species in the local target area will be larger (electronic supplementary material, figure S5). In such a situation, total biodiversity loss over time will tend to be higher, and it takes a longer time to reach asymptotic values (figure 2). By contrast, when parameter ω in Fisher's logseries model is higher, the corresponding total biodiversity loss quickly becomes asymptotically stable because of low species richness in the target area ( figure 2 and electronic supplementary material, figure S5). Note that the curvilinear trends were very similar between the proposed model (equation (2.1)) and the full immigration model (electronic supplementary material, equation S1).
When the emigration rate of a species is high, the expected total biodiversity loss is also high ( figure 3). Moreover, for large emigration rates, 0.1 and 0.2, the time to equilibrium is very short as their expected total biodiversity losses keep unchanged after habitat destruction. However, the expected total biodiversity loss can increase first then decrease for the small emigration rate cases (u = 0 and 0.01). By contrast, increasing the immigration rate (v) can remarkably retard biodiversity loss in both the proposed (equation (2.1)) and the full stochastic immigration (electronic supplementary material, equation S1) models (figure 4). More interestingly, when looking at the temporal trajectory of species loss curve in figure 4, potential immigration credits can be observed (blue and green solid curves when v = 0.1 and 0.2, respectively). Analogous immigration credits were also observed for the full immigration model (blue and green dashed curves when v = 0.1 and 0.2, respectively, in figure 4).
Finally, when the intact area of interest has a larger areal size, the expected total biodiversity loss tends to be higher after habitat destruction ( figure 5). However, this positive relationship is nonlinear. When the local intact area has a sufficiently large areal size, the total biodiversity loss tends to reach the maximal asymptotic value and only slightly increases when the size of the local target area is

Discussion
Modelling biodiversity change, because it involves the coupling of spatial and temporal ecological factors, is an intriguing but challenging topic for global ecologists. Recent studies suggested that both  ecological and evolutionary processes can take place on the same time scale [23,24]. To this end, our present model serves as a bridge to disentangle how ecological and evolutionary processes jointly explain the temporal complexity of the extinction of species. In our study, we quantified temporal biodiversity loss pattern, which differs from previous studies that examined extinction debts. However, both concepts are closely related and become complementary to each other when discussing total extinction. That is, the sum of total biodiversity loss and extinction debt is equal to the total expected extinction when species richness reaches equilibrium after a sufficient time. As a consequence, our study actually already quantified temporal behaviours of species loss or the probability of extinction over time (electronic supplementary material, figure S3).
The parameter characterizing the immigration effect (v) can mitigate the potential biodiversity loss, while, by contrast, the parameter characterizing the emigration effect (m) would accelerate biodiversity loss (figures 3 versus 4). To this end, bidirectional dispersal processes (emigration versus immigration) could have great but contradictory impacts on ecosystem stability by accelerating versus reducing species extinctions, respectively (figures 3 and 4 and electronic supplementary material, figure S3). From a perspective of conservation biology, it is expected that the immigration effect represents a form of a rescue effect to maintain biodiversity in the studied local area by providing individuals of extinct species from the outside to the local community of interest. In our model (equation (2.1)), we assume that the arrival of a new immigrant is only possible when there are no individuals of the same species in the target community. This means, we assume that strict competitive exclusion effect exists between native individuals and immigrants, and the immigration event is possible when there are no native conspecifics in the community. However, native conspecifics (i.e. those individuals that are created by birth events) can coexist without exclusion.
This assumption seems a bit strong but realistic and a bit similar to the rescue effect in ecological literature [25][26][27][28][29][30]. In our paper, the rescue effect takes place through immigration from other areas only when the population size of a target species in a target area is sufficiently small. That is, the rescue effect (v) is magnified at most when there are no individuals of the target species in the target area. By contrast, the rescue effect (v) is ignored in our model (equation (2.1)) when there are some individuals of the target species in the focused area. This is because the dynamic of population is dominated by density-dependent growth and death rates (both rates = 1, see table 1 for details). The density-independent immigration rescue effect is therefore negligible, as v is usually small (table 1). Previous studies [27][28][29][30] working on the rescue effect also emphasized the importance of rescue effect for species that are extinction-prone and in small population sizes, akin to our model structure and assumption. However, the difference between theirs and ours is that we assume the rescue effect or immigration only takes place when there are no individuals of the species in the target area. Additionally, the model structure and assumption implemented in our paper are also a bit similar to the one-migrant-per-generation rule in conservation biology and genetics [33][34][35]. In our study, one migrant (vp 0 ) is only allowed for initiating the dynamic of a target population of a species in the studied area.
Finally but most importantly, even a more realistic immigration model is used (electronic supplementary material, equation S1 and figure S1B of the appendix), one can see that its resultant curve shape patterns were similar to those of our proposed model (figures 1, 2 and 4). Therefore, the proposed model (equation (2.1)) is sufficient for us to explore the potential roles of immigration versus emigration on biodiversity change over time. As summarized below, the present model reports the temporal trajectory of species loss, showing the emergence of immigration credits after habitat loss (figure 4) for the first time.
Previous studies extensively showed how spatial distributional patterns ( particularly aggregation) can affect species extinctions at both the local and regional scales; when the spatial distributions of species vary from random to aggregated, either net losses or immigration credits are possible [6]. However, predictions of local or regional species losses or gains in those studies were debatable given that a clear statistical sampling framework in those papers was lacking. To resolve this statistical issue, the present study goes a step further to provide an area-based Fisher's model in which the localregional sampling process is properly developed. Our numerical results showed that time-dependent immigration credits could occur when there is a mass of migrants from neighbouring regions (i.e. outside the targeted area; figure 4).
Furthermore, if species abundance distributions can be altered after habitat destruction, species gain might possibly be observed, as found in a previous study [6]. For the time-dependent model studied here, if the emigration rate of species is allowed to drop after habitat destruction, species gain is also possible (electronic supplementary material, figure S6). This phenomenon has never been observed in previous studies. To this end, we argue that a new process 'emigration credit', which differs from previous studies reiterating 'immigration credit' [3,4,6], is useful in explaining species gains due to historical habitat loss. That is, two independent dispersal processes, immigration and emigration, can contribute to the process of species gain. Our model is the first one to show, because of abrupt habitat destruction, how emigration can influence species loss and gain, and how immigration credit can emerge in local habitats over time, further extending previous static findings on the emergence of immigration credits in local habitats [6]. density dependence yes, density dependence is reflected by the fact that the changing rate of probability at abundance j is dependent on the probability at abundance j − 1, j + 1 and j; and the associated magnitude ( j + 1), ( j − 1)(1 − u) and j(2 − u). Such a kind of modelling density dependence in neutral models is also seen in previous studies [31,32] interaction with other species no, the model assumes the dynamic of a single species irrespective of the dynamic of other species immigration rate v, which is a small value in comparison to the intrinsic growth or death rate (which is equal to 1) in our model emigration rate u, which is a small value in comparison to the intrinsic growth or death rate (which is equal to 1) in our model royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 191039 As discussed above, our conceptual framework is flexible, and the trajectory of biodiversity changes over time is directly relevant to the SAD (figure 2). This means that in addition to Fisher's logseries model studied in the present work, other ecologically interesting models of SAD (e.g. the Poisson lognormal model) can be incorporated and compared to evaluate how biodiversity will change over time. Furthermore, spatial aggregation patterns of species distributions might be incorporated as in previous studies [5,6,[36][37][38][39][40]. Because our central goals were to integrate both ecological and evolutionary processes in modelling species extinction, we did not explicitly investigate how spatial distributional aggregation would influence temporal biodiversity changes. However, spatial aggregation is indeed an essential factor affecting the magnitude of extinction debt, because it is closely related to habitat fragmentation. Through numerical simulations, previous studies [37][38][39][40] have assessed the impacts of spatial percolation and scales on area-based estimation of species diversity and extinction. Given that traditional species-area relationships usually fail to capture the difference on the short-versus long-term species loss and the degree of habitat fragmentation [40,41], for future research, it becomes interesting to assess how complex spatial structures influence the temporal dynamic of species extinction.
Data accessibility. R code for generating figures 1 and 2 of main text is available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.n710qb3 [42].