On the decline of biodiversity due to area loss

Predictions of how different facets of biodiversity decline with habitat loss are broadly needed, yet challenging. Here we provide theory and a global empirical evaluation to address this challenge. We show that extinction estimates based on endemics–area and backward species–area relationships are complementary, and the crucial difference comprises the geometry of area loss. Across three taxa on four continents, the relative loss of species, and of phylogenetic and functional diversity, is highest when habitable area disappears inward from the edge of a region, lower when it disappears from the centre outwards, and lowest when area is lost at random. In inward destruction, species loss is almost proportional to area loss, although the decline in phylogenetic and functional diversity is less severe. These trends are explained by the geometry of species ranges and the shape of phylogenetic and functional trees, which may allow baseline predictions of biodiversity decline for underexplored taxa.

Beware: PD that goes extinct ( ) with loss of any area is not There is a hypothetical terrestrial species living in the southern part of the island (light green area); in the past there used to be no water and the species also used to live south from the island (dark green area), and it will be able to re-occupy the potential range if the water disappears. Panel (b) shows how the island coast truncates the potential range, while in panel (c) the artificial plot boundary truncates the realized range. When estimating extinction, one may be interested in extinction of the species only within the plot (c), within the whole island (b), or within the whole world. In any of these cases the species range, realized or potential, can be truncated by natural or artificial barriers, making Model 4 a useful null model of range placement.  Figure 2g. Here we show that, in Model 4, P out ≥ P in even when r r > r tot . Note that, in contrast to Figure 2g, the x-axis here has log 10 scale; hence, valueg of −2 means that r r is 100× smaller than r tot , 0 means that r r = r tot , and value of 2 means that r r is 100× larger than r tot .

Supplementary Tables
Supplementary Table 1 | Pearson's correlations of the regional characteristics that were used to predict AU C and ∆AU C (see the Methods for details). Cells marked by grey fill contain values higher than 0.5 or lower than -0.5 (an arbitrary threshold).  For framing magnitudes of extinction with habitable area loss (using EAR) two aspects of spatial scale are important: extent and grain. The critical spatial extent for assessing extinctions is global, since a global loss is irreversible. Grain of habitable area loss can be defined as the size of the smallest contiguous blocks of habitable area that are lost; it is the smallest spatial unit of the loss which is represented by the area axis (x-axis) of EAR, and in such definition small (local) grains are again irrelevant for extinction estimates. An alternative view is that grain is the spatial resolution at which complex shapes of immediate loss (e.g. networks of logging roads) emerge, and in such view the fine-grain configuration of loss can be critical. However, examining such configurations of loss is only possible at small extents as it requires comparably fine-grain information on species distributions.
At global extent and for the best known taxa the distributional data are only available at large grains above roughly 110 × 110 km (1 degree), and hence any study of empirical EARs at continental to global extents can only examine patterns of loss at grains of 1 degree and coarser.

Supplementary Note 2 | The difference between Models 3 and 4, with additional derivations, details and illustrations
Here we show, in richer detail than in the main text, that the difference between the inward and outward extinction curves produced by Models 3 and 4 is inevitable. Moreover, the two models give inverse relative positions of the inward and outward curves. Specifically, we show that Model 3 (with mid-domain effect) always gives a steeper extinction curve for the outward area loss relatively to the inward loss. In contrast, Model 4 (without mid-domain effect) always leads to a steeper or equally steep inward extinction curve, relatively to the outward curve.

Notation
In this section we use the following abbreviations and terms, which are also illustrated in Figure  2e in the main text.
R -The entire circular region that consists of the inner and outer domain. r r -In Model 3 this is the radius of a circular species' range that is placed entirely within the boundaries of R. In Model 4 this is the radius of a circular potential species' range that is placed into the region R, and possibly truncated by its boundary, so that a realized range emerges.  A in -Area of the inner circular domain, which is A in = πr 2 in .
A out -Area of the outer circular domain (annulus), which is A out = πr 2 tot − πr 2 in .
P in -In Model 3 this is the probability that the range with radius r r will fall entirely within the inner circular domain. In Model 4 this is the probability that the potential range with radius r r will produce a realized range falling entirely within the inner circular domain.
P out -In Model 3 this is the probability that the range with radius r r will fall entirely within the outer circular domain (annulus). In Model 4 this is the probability that the potential range with radius r r will produce a realized range falling entirely within the outer circular domain (annulus).
P overlap -In Model 3 this is the probability that the range with radius r r will overlap the boundary between the inner and the outer domain. In Model 4 this is the probability that the potential range with radius r r will produce a realized range overlapping the boundary between the inner and the outer domain.

Model 3 with mid-domain effect
In this model we are randomly placing a contiguous circular range entirely within region R, so that it never overlaps its outer boundary (Supplementary Figure 4a). This model inevitably leads to what is known as the mid-domain effect -a higher concentration of range centroids in the centre of the region than around its boundary. Here we provide formulas for the following probabilities: The probability P out that a range with radius r r will fall entirely within the outer domain is: X out is the ratio of the outer area A out to the total area A tot , both adjusted by r r , which eliminates the possibility of the range overlapping the region R's outer boundary, and which is a sufficient condition to invoke the mid-domain effect.
The probability P in that the range with radius r r will fall entirely within the inner domain is: where X in is the ratio of the inner area to the total area, both adjusted by r r to invoke mid-domain effect.
The probability P overlap that the range will overlap the boundary between the inner and outer domain is:

Model 4 without mid-domain effect
In this model the circular range (which we call potential range) is placed randomly into region R, it is allowed to overlap the R's outer boundary, and the part of the range that ends up outside of the boundary is eliminated (Supplementary Figure 4b). Hence, the resulting realized range is only the cropped part of the original range that lies inside R. The probability P out that the realized range lies entirely within the outer domain is: X out is the ratio of the outer area to the total area, both adjusted by r r of the potential range to eliminate the mid-domain effect. Note that here the r r is added, rather then subtracted from r tot as in the 'mid-domain' case of Model 3. Also note that in Model 4 the centres of the potential circular ranges may lie close to, or even beyond, the boundary of R (Supplementary Figure 4b).
The probability P in that the realized range lies entirely within the inner domain is: X in is the ratio of the inner area to the total area, both adjusted by the r r of the potential range to eliminate the mid-domain effect.
The probability P overlap that the range will overlap the boundary between the inner and outer domain is: The proof: P in and P out as functions of range radius In order to have the the effect of inner (A in ) and outer (A out ) areas fully under control, and hence to focus fully on the difference between the inward and outward direction of destruction, we set A in = A out . Under this condition it follows that r in = r tot / √ 2, which allows us to calculate P in , P out (equations 1, 3, 6, 8) and P overlap for any r r .
For the total number of S tot species that live in the region R, and are indexed by i, it follows that E in = Stot i=1 P in i and E out = Stot i=1 P out i . In other words, to get the mean number of species E in that go extinct with the loss of A in we need to sum up the probabilities P in i for all species that live in the region R. Thus: • For Model 3 it always holds that P in ≥ P out , from which follows that E in ≥ E out . This holds for any range size and for any imaginable frequency distribution of range sizes, as the curves in Figure 2f involve all possible range sizes (represented by r r ).
• Model 4 is somewhat trickier since the potential circular range can be truncated by R's boundary, and hence the r r of the potential range is not directly proportional to the area of the realized range inside R (Supplementary Figure 4b). Figure 2g shows that, for potential ranges with 0 < r r ≤ r tot , P in < P out and hence E in < E out . Supplementary Figure 6 extends this reasoning for potential ranges with r r > r tot , showing that E in < E out holds for potential ranges with r r roughly up to 100× the size of r tot ; above that threshold we see that E in ≈ E out (Supplementary Figure 6). This means that any imaginable frequency distribution of potential range sizes in which at least one potential range satisfies 0 < r r < r tot × 100 will lead to E in < E out , and E in ≈ E out otherwise 1 .

Probability that the range is lost as a function of lost area
To further illustrate the discrepancies between Models 3 and 4 we can use equations 1, 3, 6 and 8 to plot the probability that a randomly placed range will be lost (P in or P out ) as a function of lost area (A in or A out ) -this contrasts with the previous section where A in = A out . We also set A tot = 1, leading to r tot = 1/π.
• In case of the outward area loss we eliminate the inner circular domain and calculate the probability that the species range occurs entirely within the inner domain; hence, for a given range perimeter r r , we plot P in (from equations 3 or 8) against A in .
• In case of the inward area loss we eliminate the outer domain (annulus) and calculate the probability that the species range will occur entirely within the outer domain; for a given range perimeter r r we plot P out (from equations 1 or 6) against A out .
We have chosen to plot the curves for r r of 0.01, 0.1 and 0.5 (Supplementary Figure 7), and in all three cases the results are consistent with the simulations (Fig. 2c-d): In Model 3 the outward area loss leads to higher probability that a species is lost, while the opposite holds for Model 4.

Supplementary Methods
Functional traits used to calculate functional diversity The trait categories were composed of single (body size, activity time) or multiple variables (diet, foraging niche).
Diet was characterized as proportional use of each of seven dietary categories for mammals (seeds, fleshy fruits, nectar and pollen, other plant material, invertebrates, fish, vertebrates) and eight dietary categories for birds (seeds, fleshy fruits, nectar and pollen, other plant material, invertebrates, fish, carrion, other vertebrates).
To represent variability in daytime activity patterns we matched bird foraging height data to that of mammals in the form of four ordinal categories: (1) ground level, (2) scansorial or low vegetation or understory, (3) fully arboreal or canopy, (4) aerial.
Finally, we classified the species by the environment in which they forage to: (1) aquatic, (2) semi-aquatic, (3) terrestrial or non-aquatic.
For additional details see Supplementary References 1 and 2.
Details of the calculation of the regional predictors • Gamma statistic γ. We used function gamStat in R package laser to calculate the γ.
• Colless index. We calculated the index by using the colless function in R package apTreeshape.
• Mean Moran's I of the ranges. We measured the global Moran's I 3 using function Moran in R package raster (using the default 3 × 3 'King-style' neighborhood). We note that this measure will always be inevitably correlated with average range size (Supplementary Table 1), and hence a region with high mean Moran's I of the ranges will also have large mean range area.
• Richness gradient. We created map of species richness for each taxon in each region. We then took geographic coordinates (latitude and longitude) of each grid cell and calculated an ordinary least squares regression (function lm in R) of the richness at the cells against latitude and longitude of the cells, and their interaction. In the R-style formula specification the model would be written as: lm(richness ∼ latitude + longitude + latitude:longitude). We then calculated R 2 of this regression, which is our index of richness gradient.
• Moran's I of richness. We measured the global 1st distance class autocorrelation of cell-specific values of species richness by the global Moran's I (using the 3x3 'King-style' neighbourhood). We used the function Moran in R package raster for that.