Species-area relationships and biodiversity loss in fragmented landscapes.

To estimate species loss from habitat destruction, ecologists typically use species-area relationships, but this approach neglects the spatial pattern of habitat fragmentation. Here, we provide new, easily applied, analytical methods that place upper and lower bounds on immediate species loss at any spatial scale and for any spatial pattern of habitat loss. Our formulas are expressed in terms of what we name the 'Preston function', which describes triphasic species-area relationships for contiguous regions. We apply our method to case studies of deforestation and tropical tree species loss at three different scales: a 50 ha forest plot in Panama, the tropical city-state of Singapore and the Brazilian Amazon. Our results show that immediate species loss is somewhat insensitive to fragmentation pattern at small scales but highly sensitive at larger scales: predicted species loss in the Amazon varies by a factor of 16 across different spatial structures of habitat loss.

Introduction

To estimate the number of species lost as a result of habitat destruction, ecologists typically turn to species–area relationships (SARs), which take habitat area as an input and give species richness as an output. This approach has several limitations (Lewis 2006). First, it ignores the dependence of species loss on the time elapsed since habitat destruction occurred. Here, we study the shortest timescale possible – immediate species loss – which is a prerequisite to understanding more complex long‐term problems where losses from extinction debt repayment are added to immediate losses. Second, while traditional SARs assume that species persist only in remnant habitat areas, in practice some species may persist in non‐habitat areas too. Previous studies have tackled this second problem (Pereira & Daily 2006). Here, we retain the assumption that species live only inside habitat areas in order to focus on the third major limitation of the traditional species–area approach: the assumption that total habitat area is the only spatial parameter of importance. In practice, fragmentation and spatial structure – of both habitat loss and species distributions – are also critical.

If one has complete information on the spatial distribution of species, including their abundances, one can compute the number of species and the number of endemics in an arbitrary region of habitat; the number of endemics gives the number of species immediately lost when that habitat is destroyed (Kinzig & Harte 2000; He & Legendre 2002; Green & Ostling 2003; He & Hubbell 2011, 2013; Pereira et al. 2012). In practice, however, such complete information is usually lacking. The traditional fallback approach has been to estimate species loss from an SAR by comparing the species richness calculated from an initial area to that from a reduced area. But this is only accurate if one of two restrictive conditions is satisfied: either species are distributed randomly in space, or the SAR is constructed such that the spatial pattern of habitat considered is exactly opposite (complementary) to the spatial structure of destroyed habitat in the target system (He & Hubbell 2013). These considerations are routinely ignored in practice, with the typical approach being to parameterise an SAR with data from contiguous patches of habitat and apply it to a fragmented target system.

The above concerns apply to any SAR, that is, any function that expresses species richness solely as a function of habitat area. The most common specific approach is to use a power‐law function of area (Arrhenius 1921) with exponent in the range 0.2–0.3. But these typical exponents come from fits to data from contiguous habitat patches (e.g. Storch et al. 2012), which are not complementary to the complex and usually fragmented patterns of real world habitat destruction and thus not suitable for predicting species loss. Further problems with the power‐law SAR are that it is phenomenological, rather than mechanistic, and that it ignores scale‐dependent variation in the power‐law exponent observed in empirical data (Harte et al. 2009).

Beyond the power‐law SAR, more sophisticated mechanistic SAR models have been developed (Durrett & Levin 1996; Rosindell & Cornell 2007; O'Dwyer & Green 2010; Grilli et al. 2012). Such models accurately characterise the contiguous SARs seen in nature, with scale‐dependent power‐law exponents and an overall triphasic shape (Preston 1960; Plotkin et al. 2000; Storch et al. 2012), but share the one key limitation with the power‐law SAR when it comes to predicting species loss from habitat destruction: they implicitly assume contiguous habitat.

Several studies have pushed the envelope and studied fragmented SARs, which allow one to predict how species loss is affected by the spatial pattern of habitat fragmentation, not just total habitat area loss. These studies have typically used phenomenological extensions of the power‐law model (Kinzig & Harte 2000; He & Hubbell 2011; Pereira et al. 2012) or spatial simulations (Campos et al. 2012; Hanski et al. 2013; Keil et al. 2015), and have mostly found that fragmented habitat destruction causes less immediate species loss than does contiguous habitat destruction, because fragmented habitat arrangements preserve more of a landscape's beta diversity. Indeed, simulations suggest that immediate species loss is lowest when the landscape is cleared randomly (Claudino et al. 2015). Simulations also suggest that fragmented SARs exhibit a qualitatively similar triphasic shape to contiguous SARs (Campos et al. 2013). But to establish the generality of these conclusions and to facilitate species loss calculations at spatial scales where simulations are infeasible, we require analytical results from mechanistic models, which have thus far been elusive. The only such result we know of is for when habitat destruction is spatially random and each species’ spatial distribution follows a finite version of a negative binomial distribution (Zillio & He 2010; He & Hubbell 2011).

Here, we make substantial progress towards the problem of analytically calculating immediate species loss in fragmented landscapes. Our analyses are based on spatially explicit communities generated by a neutral model (Hubbell 2001; O'Dwyer & Chisholm 2013). Neutral theory is not the focus of our study; we apply it here because it generates artificial communities that look similar to natural ones while being analytically tractable across multiple scales. Our results place upper and lower bounds on immediate species loss at any spatial scale and for any spatial pattern of habitat destruction.

Methods

Clearing patterns

Empirical habitat destruction typically exhibits complex fractal‐like patterns (Gardner et al. 1987; Sugihara & May 1990; With 1997; Hargrove et al. 2002; Ewers & Laurance 2006). In theoretical studies, patterns used include fractal structures (With 1997), percolation maps (Campos et al. 2012) and distinct clumps of habitat or non‐habitat (He & Hubbell 2011). Some results we derive here apply to all clearing patterns but we focus on three cases: ‘contiguous clearing’, where remaining habitat is a contiguous block (Fig. 1a inset); ‘clumped clearing’, where remaining habitat comprises several clumps that are randomly distributed throughout a landscape (Fig. 1b,c insets); and ‘random clearing’, where the clearing pattern is random at the scale of individual organisms, with no spatial autocorrelation (Fig. 1d inset). Random and contiguous clearing provide upper and lower bounds on species richness, and hence lower and upper bounds on immediate species loss (Claudino et al. 2015). Heuristically, this is because remaining habitat under random clearing is sampled uniformly across the landscape, maximising beta diversity, whereas contiguous clearing maximises the spatial autocorrelation among remaining habitat, minimising beta diversity.

Figure 1

Insets: Types of habitat clearing considered in this study. Green pixels indicate habitat; grey pixels indicate cleared areas. Panel a shows classic contiguous clearing (n = 1); panels b and c show clumped clearing with n = 16 and n = 256 clumps respectively; panel d shows random clearing (n = A). Main panels: SARs resulting from the four types of habitat clearing shown in insets, for a range of values of landscape extent A and total habitat area A = A /16 (thus A and A vary proportionally on the horizontal axes). The ecological communities to which the clearing maps are applied are spatial neutral communities with speciation rate ν = 0.01 and a range of values for the dispersal standard deviation σ (colours; see legend on panel a).

To characterise our clearing patterns, we use three parameters: the area of remaining habitat A; the total area of the landscape of interest A  ≥ A; and the number of fragments n, with n = 1 corresponding to contiguous clearing, n = A to random clearing, and 1 < n < A to clumped clearing. We constructed the clumped clearing maps by randomly choosing n clump centroids on a gridded landscape of area A , and then marking as habitat the A grid cells nearest to centroids.

Spatial neutral model

To create artificial and generic communities on which to study the immediate effects of habitat loss, we required a spatially explicit model that produces plausible species distributions and is tractable for large landscape sizes. We chose the spatial neutral model on an infinite two‐dimensional landscape, which comes in two versions: zero‐sum and non‐zero‐sum. In the zero‐sum version (e.g. Rosindell & Cornell 2007), exactly one individual occupies each grid cell. Individuals rain propagules down in a pattern governed by a radially symmetric dispersal kernel with variance σ2. In most of the following we assume that the kernel is a bivariate normal distribution, but many of our results generalise to any dispersal kernel with finite variance (σ2 < ∞). Individual death follows a Poisson process, and each cell vacated by a death is immediately replaced by a new recruit drawn at random from the propagule rain onto the cell. There is a small but non‐zero probability ν that the new recruit transforms into a new species, that is, ν is the per‐capita probability of point speciation. The non‐zero‐sum version of the model (e.g. O'Dwyer & Cornell 2017) is similar, except that deaths and births are independent processes, positions of individuals are not constrained to a grid, and new species enter the community via an independent Poisson process whose rate parameter ν is constant in time and space and measured per unit area.

Of the two versions of the spatial neutral model, the zero‐sum is generally easier to simulate, and non‐zero‐sum is more tractable analytically. But they yield almost identical results except in some limiting cases, for example, σ2 → 0. In particular, they reproduce the triphasic contiguous SAR characteristic of empirical data (Rosindell & Cornell 2007), with the first phase occurring at scales A < σ2, the second phase at σ2 < A < σ2/ν and the third phase at A > σ2/ν (O'Dwyer & Green 2010). Here, we draw on coalescence theory to simulate the zero‐sum model (Rosindell et al. 2008) and analytical solutions for the non‐zero‐sum model (O'Dwyer & Cornell 2017; also see Appendix 1).

Because we are here concerned with immediate species loss, we solve the spatially explicit neutral model at the dynamic equilibrium on a contiguous landscape, and then subsequently fragment the community and measure immediate species loss (if we were instead investigating long‐term species loss, we would need to find the dynamic equilibrium on the fragmented landscape).

Our results from the neutral model can be applied to non‐neutral systems provided that the constituent species’ spatial distributions are statistically similar to those of neutral species. In such applications, the parameters σ and ν can be chosen to match data on species’ spatial distributions, with σ controlling spatial autocorrelation and ν controlling overall landscape diversity (Condit et al. 2002) (see Appendix 2 for details).

Rescaling approach

We now define functions that relate species richness immediately after habitat clearing to the parameters of the clearing map and the community model. Except where otherwise noted, the sampling area A, landscape area A and dispersal parameter σ are measured in units such that there is one individual per unit area. For the contiguous‐clearing case, we define , which is a contiguous SAR and depends only on the area A of a contiguous habitat region and the two underlying community model parameters. For the random‐clearing case, we define , which additionally depends on the underlying landscape area A . The clumped‐clearing case depends also on n and is described by .

The following rescaling relationship for has been derived previously (Rosindell & Cornell 2007) for σ2 ≫ 1:for some function Ψ, or, equivalently,for any scalar k. Formula (1) reduces the number of effective parameters so that instead of having to understand the three‐parameter function S , we need to only understand the simpler two‐parameter function Ψ.

In view of the central importance of the neutral contiguous‐SAR function in ecological theory, we henceforth refer to Ψ as the ‘Preston function’, after Frank Preston who first observed the triphasic SAR in empirical data and understood the basic mechanisms behind it (Preston 1960). Specifically, we define the Preston function to be the SAR for contiguous circular areas in the non‐zero‐sum spatial neutral model with a bivariate normal dispersal kernel (O'Dwyer & Green 2010; O'Dwyer & Cornell 2017). We choose this definition because the non‐zero‐sum model is more amenable to analytical treatment and, anyway, produces almost identical SARs to the zero‐sum model. Abstracting the Preston function away from the main focus of our work enables us to derive rescaling formulas independently of the mathematical study of the Preston function itself (O'Dwyer & Green 2010; Grilli et al. 2012; O'Dwyer & Cornell 2017). The combination of the new scaling results we present below and any method to evaluate the Preston function yields analytical approximations of S , S and S that predict immediate species loss under a range of scenarios at arbitrary spatial scales.

Case studies

After deriving our rescaling formulas, we applied them to three real species‐loss problems spanning a range of spatial scales. We focused our applications on tropical tree communities, for several reasons: (1) Beta diversity of tropical trees over distances of 0.5–50 km is statistically similar to beta diversity under neutral theory (Condit et al. 2002), which suggests that tree species’ spatial distributions at these scales are also statistically similar to those under neutral theory (notwithstanding a current paucity of individual‐mapped data at these scales). (2) Trees are sessile and this simplifies the question of immediate species loss because individuals cannot relocate to escape habitat destruction. (3) Immediate species loss is most relevant to communities in which extinction debt is small or repaid slowly, and in tree communities the latter condition is likely to be satisfied because trees are sedentary and long‐lived (the same is not true for, e.g. vertebrates). (4) There are substantial data available for tropical tree communities.

Our definition of ‘tree’ is a vascular plant that can grow to at least 10 cm diameter‐at‐breast‐height (DBH). For all case studies, we used estimates of tree density and dispersal distance from the well‐studied tropical forests of central Panama: our estimate of tree density was ρ = 42 000 km−2 (Hubbell et al. 2005); our estimate of dispersal standard deviation was σ = 40.2 m (Condit et al. 2002). Although in reality ρ and σ vary across tropical forests, these ballpark estimates serve the purpose of demonstrating our methods and producing approximate estimates of species loss.

For each case study, we first determined from empirical data, the initial area A , post‐clearing area A, and initial species richness S . The values of A and S were used as constraints to estimate the speciation rate parameter ν separately for each case study (Appendix 2). Initial species richness, initial area, final area and dispersal distance were thus the only data used to fully parameterise our model: our predictions were emergent from these data and not fitted to the empirical species loss numbers. Next, we used our analytical formulas to estimate the post‐clearing species richness S and S under contiguous‐ and random‐clearing scenarios, respectively, and thus the number of endemic species in the cleared areas, and , under the same two scenarios. Species loss under each scenario can then be equated to the number of species that were endemic to the cleared region (He & Hubbell 2013): and . These provide upper and lower bounds on immediate species loss under arbitrary clearing scenarios (Claudino et al. 2015).

The first case study was at the local scale of an A  = 50 ha forest plot on Barro Colorado Island in the Panama Canal (Hubbell et al. 2005). The number of tree species in the plot in 2010 was S  = 222. Our clearing scenario here was hypothetical, with area losses of 10–90% explored. We compared these estimates to those from simulated random and contiguous clearing on the actual 50 ha plot spatial census data.

Our second case study was at the regional scale and corresponded to the entire tropical island city‐state of Singapore (Corlett 1992). Prior to 1819, Singapore was almost entirely forested. From 1819 to 1992, Singapore's forested area decreased from approximately to (Corlett 1992), representing 95.7% forest loss. From herbarium records, we estimated Singapore's tree species richness in 1819 as S  = 840 and current tree species richness as S = 728, implying ΔS = 112 extinctions. In addition, we estimated the number of undetected extinctions (species that went extinct before they could be recorded) from a nonparametric statistical method (Chisholm et al. 2016) that requires as input only the time series of herbarium records. The method estimated 89 undetected extinctions and a corresponding upper bound on species loss of ΔS = 112 + 89 = 201.

Our third case study was at the continental scale and covered the Brazilian Amazon. The true magnitudes of tree species diversity and extinction in the Amazon are highly uncertain. Accordingly, this case study is not an exercise in model validation, but a demonstration of our formulas’ power to make predictions easily and quickly at arbitrarily large spatial scales. We used input parameter values from the optimistic scenario of Hubbell et al. (2008): original area A  = 4 652 400 km2, current area A = 0.633A and original tree species richness S  = 11 210. We also compared our model's predictions to those of Hubbell et al. (2008), who used a more complex spatial simulation approach.

In Appendix 3, we provide R code for reproducing our case study results and for applying our methods more broadly. As input, the R code requires five parameter values: initial habitat area A , final habitat area A, initial species richness S , density of individuals per unit area ρ and standard deviation of dispersal distance σ. The parameter σ can be estimated either directly from dispersal data (e.g. from seed trap data for plants (Muller‐Landau et al. 2008) or telemetry data for animals) or indirectly from data on the spatial distribution of organisms, as done for our case studies here (Condit et al. 2002) (Appendix 2). The R code returns estimated species loss under contiguous and random clearing (running time is on the order of milliseconds).

Results

Changing spatial scale

Our first theoretical result (Appendix 4) is that rescaling relationships similar to those shown in formula (2) also apply to arbitrary clearing patterns, including clumped clearing and random clearing:

These formulas can be understood intuitively in terms of coalescence processes (Appendix 4). We confirmed them numerically for the clumped clearing case (see Fig. 1b,c and dashed and dash‐dotted curves in Figs 2 and 3), the random clearing case (see Fig. 1d and dotted curves in Figs 2 and 3) and for arbitrary clearing patterns (Fig. 4).

Figure 2

The SARs from Fig. 1 rescaled by 1/σ2 on both axes. Solid curves correspond to contiguous clearing (Fig. 1a); dashed and dash‐dotted curves correspond to clumped clearing with n = 16 and n = 256 clumps respectively (Fig. 1b,c); dotted curves correspond to random clearing (Fig. 1d). As predicted by formulas (2), (3), (4), there is a scaling collapse within each of the four clearing types and the collapse is better for large σ (in this case, for σ ≥ 2). The shaded grey region indicates the second phase of the SAR: σ2 < A < σ2/ν. (See Fig. 3a for a zoomed‐in version of the second and third SAR phases of this graph, without the σ = 1 curve).

Figure 3

Rescaled SARs as in Fig. 2 but showing only the second and third SAR phases, and also omitting the σ = 1 case (which does not exhibit a good scaling collapse; see Fig. 2). Panel a is for the same parameter values as Fig. 2 (habitat‐to‐landscape area ratio A/A  = 1/16 and speciation rate ν = 0.01); Panel b has a lower speciation rate ν = 0.001; Panels c and d repeat panels a and b with a lower habitat‐to‐landscape area ratio A/A  = 1/256.

Figure 4

Arbitrary clearing patterns (left) and a test of the hypothesis that the species richness after clearing (S) divided by the variance of the dispersal kernel (σ2) should be invariant under the spatial rescaling of both the clearing map and σ, provided that σ is not too small (the clearing map is ‘stretched’ equally in all directions proportionately with increasing dispersal distance σ; see Appendix 4 for details). Clearing patterns (from top‐left, moving clockwise) are ‘smiley’, ‘illusion’, ‘fractal’, and ‘panda’. A ‘clumped’ clearing pattern with 16 clumps (Fig. 1b inset) is also included. Each data point shows the mean species richness for the corresponding clearing pattern (colours) applied to 10 replicate communities (speciation rate ν = 0.01, σ values on horizontal axis).

Relating clumped clearing to other clearing patterns

We also derived theoretical rescaling relationships between qualitatively different types of clearing. Three different results emerged relating a clumped‐clearing pattern (see Fig. 1b,c insets) to a contiguous‐clearing pattern (see Fig. 1a inset). The first result occurs when the mean species range size is small so that any given species is unlikely to appear in more than one clump: A /n ≫ σ2/ν. This decouples the clumps from one another so that total species richness is simply the sum of the species richness in each independent clump:

Our second result for clumped clearing occurs under the stricter condition that A/n ≫ σ2/ν so that the clumps themselves become larger than the average species range size. Consequently, there is negligible overlap of species composition between clumps, even if the clumps are repositioned right next to one another. Merging the clumps to form a contiguous‐clearing pattern therefore does not affect species loss and clumped clearing becomes equivalent to contiguous clearing:

This result is demonstrated in Fig. 3, where the dashed and dash‐dotted curves (clumped clearing) begin to converge to the solid curves (contiguous clearing) as A/σ2 becomes large, with faster convergence for smaller n and larger ν.

Our third result for clumped clearing occurs at the opposite extreme where the average species range size (σ2/ν) is large enough to overlap many patches, that is, A /n ≪ σ2/ν. In this scenario, species are well mixed at the scale of a few habitat clumps and hence species richness is affected similarly by random and clumped clearing:

We numerically confirmed formula (7) for simulated clumped clearing with n = 16, and n = 256. In Fig. 3, for instance, we see that clumped clearing converges to random clearing as A/σ2 (and hence A /σ2) becomes small, with faster convergence for larger n and smaller ν (see Fig. S1 for a broader range of n). Our rescaling relationships (5)–(7) assume equal‐sized clumps (Figs 1b,c insets), but we expect they will be robust to moderate perturbations of this idealised case, for example, unevenness in clump shape and size.

Relating random clearing to contiguous clearing

We derived a further scaling result relating random‐clearing to contiguous‐clearing scenarios provided that σ ≫ 1:

This, our flagship theoretical result, is highly non‐trivial. To put it in words, if a fraction A/A of the landscape remains as habitat after a spatially random clearing event, this has a similar effect on species richness to a scenario of contiguous clearing leaving the same habitat area A, but on a more species‐rich landscape constructed with a larger speciation rate. This result can again be understood in terms of coalescence processes (Appendix 5).

Formula (8) gives a rescaling from random clearing to contiguous clearing. Formula (2) gives a rescaling between contiguous‐clearing scenarios with different sets of parameters. We can combine formulas (2) and (8) to get a rescaling between random‐clearing scenarios (Appendix 5):

We verified numerically that this result holds true except for σ ≈ 1 (Fig. 5).

Figure 5

Panel a: SARs resulting from random clearing with parameters ν and σ related by , where ν0 = 0.01, σ2 = 212 k and k ∈ {2−12, 2−10, … 20}. Panel b: The same SARs rescaled by 1/σ2 on both axes. These graphs demonstrate a different scaling collapse from that in Figs 2, 3, 4, because here the landscape extent is fixed (A  = 222) and only habitat A varies along the horizontal axis, and because the speciation rate ν now varies with dispersal standard deviation σ. As predicted by formula (9), a scaling collapse occurs in panel b and the collapse is better for large σ.

Relating the results from all clearing patterns to the Preston function

The function is already related to the two‐parameter Preston function by formula (1). The equivalent result for comes from a combination of formulas (1), (2) and (8):and is broadly valid for σ ≫ 1 (Fig. S2). This shows that the original four‐parameter random‐clearing function collapses on to the well‐studied two‐parameter Preston function and allows us to transfer knowledge about contiguous SARs to random‐clearing SARs. In the limits of large area and small speciation rate, our formulas for S and S reduce to known SAR and endemics–area results for the log‐series distribution (He & Legendre 2002; Green & Ostling 2003; Appendix 6).

For clumped clearing there are two scenarios depending on the degree of dispersal limitation. If interpatch distance is large compared to the average species range A /n ≫ σ2/ν, and σ2 ≫ 1, then we can combine formulas (1) and (5) to give

In the opposite case A /n ≪ σ2/ν, we can combine formulas (7) and (10) to give

Figure S3 summarises all the rescaling relationships derived in this study. The lower bound on species richness after habitat loss, and hence the upper bound on species loss, comes directly from the classic contiguous‐SAR approach by calculating . The upper bound on species richness, and hence the lower bound on species loss, comes from random‐clearing SARs by calculating . Importantly, both can now be expressed analytically in terms of the well‐studied Preston function.

Local‐scale findings and case study: Barro Colorado Island, Panama

We observed that on local scales, random and contiguous clearing give qualitatively similar SARs (Fig. 6a); thus fragmented SARs will be similar too, because random and contiguous clearing are the bounding cases. Despite the qualitative similarity of the SARs, the quantitative differences in species loss estimates between clearing scenarios can be fairly large (e.g. by a factor of ≈ 2.5 for 90% forest loss in the BCI plot; Fig. 6b). When we compared the predictions from the formulas to the results of simulated clearing on the spatial census data, we found that the simulated results were within the bounds from the formulas, but showed less variation between the random and contiguous scenarios (Figs 6b; S4, Appendix 8). This discrepancy is expected, given the known limitations of current spatial neutral models at small scales, that is, their failure to accurately reproduce small‐scale beta diversity patterns (Condit et al. 2002) and species abundance distributions (O'Dwyer & Cornell 2017).

Figure 6

Results from our random‐clearing and contiguous‐clearing SAR formulas with realistic parameter values (see formula (8), and formula (S1) in Appendix 1). Panels a, c and e show illustrative random‐clearing SARs (dashed curves) and contiguous‐clearing SARs (solid curves) with landscape area A  = 104, 108 and 1012, respectively, and σ = 8 and ν = 10−8. Panels b, d and f show estimated tree species loss in our case studies at the local scale (hypothetical forest loss of 90% in Barro Colorado Island (BCI) 50 ha plot, Panama), regional scale (observed forest loss in Singapore, 1819–present) and continental scale (observed forest loss in the Brazilian Amazon since pre‐1970), with estimates from our formulas (dark coloured bars) compared to estimates from other methods (light coloured bars; see main text for parameterisations and other details of the case studies).

Regional‐scale findings and case study: Singapore

When clearing occurs at a regional scale, qualitative differences between the random and contiguous scenarios emerge (Fig. 6c), and the quantitative differences become larger. In our application to Singapore, estimated species loss in the contiguous‐clearing scenario is five times that in the random‐clearing scenario, and empirical tree species loss is closer to the random‐clearing lower bound (Fig. 6d).

Continental‐scale findings and case study: the Brazilian Amazon

At continental scales, where the third phase of the SAR dominates, major qualitative differences arise in the SAR across fragmentation types: while the contiguous model predicts that the proportion of species lost equals the proportion of area lost in the initial stages of habitat destruction, the random‐clearing model predicts very little species loss until the vast majority of habitat has been destroyed (Fig. 6e). The stark differences between the predictions of the random and contiguous models leave a huge envelope of possible intermediate outcomes for fragmented‐clearing scenarios (Fig. 6f).

Discussion

How many species are lost when an area of grassland goes under the plough or a forest falls by the axe? Ecologists have pondered the species–area question for almost a century (Arrhenius 1921), but have only recently come to grapple with the spatial nature of the problem, with the help of advanced computer simulations and mathematics (Durrett & Levin 1996; Rosindell & Cornell 2007; O'Dwyer & Green 2010; Grilli et al. 2012). Here, we have compartmentalised the problem by first designating the ‘Preston function’ to be the well‐studied triphasic contiguous SAR that emerges from a spatial neutral model, and then deriving new formulas that relate this existing work to a range of new fragmented SAR scenarios (Fig. S3). Our results are derived from neutral models, but their applicability extends to any situation where the spatial patterns of species distributions are statistically similar to those of neutral species – the actual mechanisms generating the patterns need not be neutral.

The spatial pattern of habitat loss can have a major influence on species loss, and traditional contiguous SARs may overestimate this loss (Kinzig & Harte 2000), but most studies on fragmented SARs have been restricted to phenomenological approaches or numerical simulations. Although useful, simulations come with computational limitations that restrict the spatial scale and resolution at which questions can be asked and the generality of conclusions drawn. While our local‐scale case studies (Fig. 6a,b) could easily be investigated numerically, our regional‐scale case studies (Fig. 6c,d) would be more difficult to simulate at the resolution of individuals, and our continental‐scale case studies (Fig. 6e,f) would be impossible with present technology: we estimate that even efficient neutral coalescence simulations of our Amazon case study would take millennia on a modern desktop machine, whereas our analytical formulas can be evaluated in milliseconds. Our method allows for rapid estimation of upper and lower bounds on species loss from fragmented habitat clearing, with the upper bound given by a contiguous SAR and the lower bound given by a random‐clearing SAR (see Appendix 7 for a technical discussion).

Our case studies reveal a strong scale‐dependence to the fragmented SAR problem: while the pattern of clearing makes only a moderate difference to species loss on local scales (Fig. 6a,b), the effect on regional and, especially, continental scales can be enormous (Fig. 6c,f). The scale of observation also affects whether empirical species loss from fragmented clearing is closer to the theoretical lower or upper bound. In the idealised case where species are distributed randomly at the scale of observation, species loss is predicted exactly from the lower bound, corresponding to the random SAR, under any spatial pattern of habitat destruction (He & Hubbell 2011). In more realistic cases, where species distributions have some characteristic scale, the random SAR formula will continue to work well providing that this characteristic scale is substantially greater than the scale of clumping in the habitat destruction spatial pattern (formula (7)). This explains why empirical tree species loss in Singapore is close to the theoretical random‐clearing lower bound (Fig. 6c,d): the average tree species’ range is much larger than the distance between forest clumps in Singapore.

At continental scales, the range of possible outcomes from fragmented clearing becomes enormous (Fig. 6e) and, lacking empirical data on species loss, we cannot be certain whether our contiguous‐clearing upper bound or the random‐clearing lower bound is more accurate. This is because the condition for the clumped‐to‐random scaling collapse (formula (7)) may fail at large scales: the typical species ranges may cover just one or a few clumps. It is also at continental scales that the power and utility of our methods become fully apparent. In the Amazon, biologists lack certainty about the number of extant tree species, let alone the extinct ones. In such cases, if we want to estimate tree diversity and extinction, we must leverage theoretical models, despite their inevitable caveats (Hubbell et al. 2008).

Two important general lessons emerge from the application of our new methods to the three case studies at different scales. First, the relationship between fragmentation and immediate species loss is fundamentally scale dependent: fragmentation has only moderate effects independent of area at small scales (Fig. 6a), but enormous effects at large scales (Fig. 6e). Second, the 16‐fold difference between our lower bound random‐clearing and upper bound contiguous‐clearing predictions at the Amazon scale (Fig. 6f) illustrates starkly that we can have no hope of even approximately characterising continental‐scale species loss unless we come to grips with the mathematics of fragmented SARs.

We focused on immediate species loss from habitat destruction here, but there is also a need to study long‐term loss (e.g. Claudino et al. 2015). Immediate species loss may be more relevant for communities of long‐lived sedentary organisms, such as trees, or small‐bodied organisms that can persist in tiny patches of habitat. But long‐term loss is important for communities of organisms that are mobile or turn over relatively quickly, for example, birds and mammals. Critically, the effects of fragmentation on long‐term loss can be opposite to those on immediate loss, with diversity negatively related to fragmentation (Hanski et al. 2013). Long‐term loss is a more complex technical problem, because it necessarily includes immediate loss as one component, but also involves a richer cocktail of biological processes, including edge effects and population structure. Our methods can be extended to study cases where long‐term relaxation occurs via neutral drift, but a comprehensive evaluation of long‐term loss will likely involve a broader suite of modelling techniques.

In summary, we have introduced a suite of formulas that allow for analytical, mechanistic estimation of lower and upper bounds on immediate species loss at arbitrary spatial scales. The resulting insights about fragmented SARs help provide a foundation for future work on this topic. Our formulas predict, consistent with other models and with intuition, that immediate species loss is lower under random clearing than under contiguous clearing, and in between for other cases of fragmented clearing. We have shown that true species loss will, all else being equal, be closer to the lower bound (random clearing) than the upper bound (contiguous clearing) for (1) well‐dispersed taxa, (2) highly fragmented landscapes and (3) small spatial scales. Our rescaling formulas advance research on fragmented SARs by providing an analytical way forward on a topic where analytical solutions have been scarce.

Authorship

RAC and JR designed the study. RAC, JR and FL performed the analyses. YSY, WWS and RC provided and collated the data. RAC and JR wrote the paper. All authors revised the paper.

Click here for additional data file.