Models of range placement

Here we describe four contrasting models of range placement that were used for our theoretical arguments (Fig. 2). The simulated realizations of Models 1–4 (Fig. 2a–d) use rectangular ranges placed on a rectangular grid, which is computationally convenient. The analytical reasoning based on Models 3 and 4 (Fig. 2e–g) assumes circular ranges placed into a circular region, which is analytically tractable28.

In Model 1, we assume random placement of scattered ranges. This model uses a region consisting of 20 × 20 grid cells of the same area, where each grid cell has a probability P i of being occupied by an i-th species, and where P i is constant across grid cells. The presence or absence of i-th species in each cell is then simulated as an outcome of Bernoulli process with probability P i . This model produces no spatial aggregation of ranges, and is similar to the random model of He and Hubbell12. When averaged across many realizations, this model gives identical EAR in and EAR out curves12. Note that He and Hubbell12 demonstrate this only implicitly and use different terminology; they also operate at fine scale using individuals as smallest spatial units, whereas here we consider rectangular grid.

In Model 2, we assume non-random placement of contiguous ranges. Real-world large-scale distributions are usually to some degree contiguous, and often aggregated. Model 2 takes this to the extreme by using only contiguous ranges and placing all of them into the upper-left corner of region consisting of 20 × 20 grid cells. As a consequence, when A in =A out the difference between E in and E out will simply reflect the number of ranges forced to fall completely into A out . This mimics, for example, the presence of strong environmental gradients or barriers constraining (truncating) species distributions, and in such a case, the difference between EAR in and EAR out is trivial.

Models 3 and 4 use random placement of contiguous ranges (Supplementary Fig. 4). Model 3 places each contiguous range entirely and randomly within the region, causing highest richness in the centre (mid-domain effect). Model 4 also places ranges randomly, but allows them to overlap domain edges, effectively eliminating the mid-domain effect (see ref. 26 for details of this algorithm).

Simulations from models of range placement

The simulations use square ranges placed on an artificial rectangular grid (region) of 20 × 20 grid cells (Fig. 2a–d; see refs 16, 26 for similar approach). In each simulation run, we place 711 ranges with sizes drawn from empirical range-size frequency distribution observed in the species of birds that inhabit 2,200 × 2,200 km2 placed in Western Africa (region AF1 in Fig. 4). In each run, we calculate mean extinction curves of the inward and outward scenario (EAR in and EAR out ), as well as scattered random destruction scenarios in which grid cells are lost stepwise, one by one, and with the same probability at each step. Each simulation run was repeated 100 times, and the resulting 100 extinction curves were averaged to produce Fig. 2a–d. The same algorithm was applied to produce the empirical extinction curves (see below).

Vertebrate distributional data

We used expert-drawn range maps for terrestrial amphibians, birds and mammals in all analyses. The range-maps were based on the IUCN assessment (http://www.iucnredlist.org/) for mammals50 and amphibians51. Distributions for birds were compiled from the best available sources for a given geographical region or taxonomic group1. See also ref. 52 for more details on the data.

Grid system and study regions

We used an equal area grid originally derived in cylindrical equal-area projection with∼1° size at the equator and a constant grid cell area A cell of ∼110 × 110 km2. We selected nine square regions (hereafter the regions, Fig. 4 and Supplementary Fig. 1), each of them having a total area (A tot ) of exactly 20 × 20 grid cells. The regions were selected so that they lie completely within major continental landmasses and each grid cell within the region contains at least some land. The regions were also selected to cover broad variety of terrestrial biomes, altitudes and latitudes. The nearly identical geometry (note that regions in higher latitudes are slightly elongated relative to the low latitudes) of the regions enabled us to control for the potentially confounding effects of area and shape, and the extinction curves are thus comparable among the nine regions. We note that because of their large size, the study regions are bound to contain distinct environmental gradients, and species ranges inside the regions may aggregate along such gradients, similar to the situation in Model 2. This has the potential to enhance the inward–outward discrepancy of the EAR curves.

Extinction curves: species richness (E)

We explored three scenarios of habitable area loss: inward, outward and scattered random. We adopted a strictly nested sampling design26 to explore the effects of inward and outward scenario (Fig. 1). We considered sampling windows of sides l and area A in =A cell × l2, where and 1≤l≤20. We started with the smallest sampling window (l=1) and with the largest area outside of the sampling window (A out =A tot −A in ). We moved the sampling window continually across the 20 × 20 regional grid, and for each position, we recorded number of species with ranges exclusively within (E in ) and number of species with ranges exclusively outside of (E out ) the sampling window. We calculated mean and by averaging E in and E out values across all of the possible positions of the sampling window. We then enlarged the side of the sampling window to l=2 and we repeated the procedure described above. We continued enlarging the window until l reached 20 and A in =A tot , calculating mean and for each of the window areas. We then plotted the proportions against A in /A tot (these are the red EAR in curves in Figs 1 and 4 representing the outward destruction) and against A out /A tot (these are the blue EAR out curves in Figs 1 and 4 that represent the inward destruction).

In the scattered random destruction scenario, we selected the grid cells for the sampling one by one. The area within the set of selected cells was A in , and the richness of species endemic to that area was E in . In each step, each grid cell within the remaining area (A tot −A in ) had the same probability of being selected for the destruction. We repeated this procedure 400 times and averaged the resulting extinction curves to get EAR in and EAR out .

Extinction curves: PD

We used the most recent and dated phylogenies on the three vertebrate taxa. For birds, we chose the first tree in the posterior set (distribution) of trees in Jetz et al.53 (see also Jetz et al.23 for discussion of tree dating and robustness). We also selected one mammal super-tree from Kuhn et al.54. For amphibians, we used the super-tree from Isaac et al.55. We used the term PD for the sum of the branch lengths of a phylogenetic tree (or a sub-tree, see below)24. In one extreme case (the rake-shaped phylogeny), PD is equivalent to species richness (Fig. 3). Let us define several terms (see also Fig. 3 and Supplementary Fig. 2): Let PD tot be the total PD of a region, PD in and PD out are PDs of all of the species occurring inside and outside of the sampling window, respectively, PDE in and PDE out are PDs of species that live exclusively (are endemic to) inside and outside of the window, and PDX in and PDX out are the PD lost due to area loss. We note that PDX in does not always equal PDE in as we demonstrate in Supplementary Fig. 2b—PDX in is calculated using PD out , whereas calculation of PDE in does not require PD out . In fact, PDX in =PD tot −PD out (Supplementary Fig. 2), and the same principle applies for PDX out . For our purpose, we calculated and , which is the PDX in and PDX out averaged over all of the positions of the sampling window of a given area. We plotted the proportions against A in /A tot (the solid PDXAR in curve in Fig. 3, representing outward habitat loss) and against A out /A tot (the dashed PDXAR out curve in Fig. 3, representing the inward habitat loss).

Extinction curves: FD

We used species-level trait databases of all birds and mammals56 to calculate dendrogram-based FD metric for all grid cell assemblages57. Specifically, we used Gower distance to calculate species pairwise dissimilarity, weighting each of five functional trait categories (diet, body size, activity time and two measures of foraging niche) equally. For detailed description of the traits, see Wilman et al.56, Additional Information and Supplementary Methods. The species were then clustered by an UPGMA algorithm58 to obtain a functional dendrogram of each taxon in each region, an approach that in simulation studies has been shown to provide strong representation of original dissimilarities59. To calculate the proportional loss of FD (FDX), we applied exactly the same procedure as described above for PD, but using the functional dendrograms instead of the phylogenetic trees. We did not have a comparably comprehensive species-level trait database for amphibians, so we did not construct the FDX curves for this taxon.

Representing regional extinction vulnerability by AUC

We use the AUC (or simply integral) to represent the abovementioned variation of extinction curves (EAR, PDXAR and FDXAR), with steep curves, that is, high extinction vulnerability of a region, characterized by high AUC.

Regional predictors

For each of the three major taxa in each of the nine regions, we put together nine predictors (see Supplementary Methods for technical details): (i) Gamma statistic. Tree ‘stemminess’ represented by Pybus and Harvey’s60 gamma statistic (γ) characterizes the distribution of branching events within the tree. Trees with γ<0 have relatively longer inter-nodal distances towards the tips of the phylogeny (‘tippy’ trees), whereas trees with γ>0 have relatively longer inter-nodal distances towards the root of the phylogeny (‘stemmy’ trees). (ii) Stemminess. This is an alternative measure of how ‘tippy’ (or ‘stemmy’) a tree is in a given region. It is calculated as L actual /L max , where L actual is the sum of all branch lengths in the phylogeny and L max is the distance from the tips to the root of the tree multiplied by the total number of tips. (iii) Colless’ index of imbalance61, which measures the branching symmetry of the phylogenetic tree of a given taxon in a given region. (iv) Mean range size, which is the arithmetic mean (calculated across all species) of the total number of grid cells in a given region in which a species was ‘detected’ according to the expert-drawn range map. (v) Mean of sqrt(range size)/perimeter. By perimeter we mean the total number of grid-cell sides that form edges of a gridded area occupied by a species in a given region. We calculated the arithmetic mean of the sqrt(range size)/perimeter over all species of a given taxon that live in the region. (vi) Mean Moran’s I of the ranges. For each species in each region, we measured the autocorrelation of the 1 (occupied) and 0 (unoccupied) values by global Moran’s I58, which is a measure of contiguity of the ranges, and we took the arithmetic mean (across all species) of the values. (vii) Total richness (or S tot ) of given taxon in given region. (viii) Richness gradient. We created fine-grain (using 1° cells) map of species richness for each taxon in each region. We then calculated an ordinary least squares linear regression of the richness in the cells against latitude and longitude of the cells, and their interaction. R2 of this regression is our index of richness gradient. We interpret it as the magnitude of the large-scale spatial autocorrelation of species richness; it also measures how clumped are species ranges towards the edge of the region. (ix) Moran’s I of richness, which measured the global first-distance class autocorrelation of cell-specific values of species richness. In contrast to the previous measure, the Moran’s I of richness captures short-distance autocorrelation.

Predictors i–iii capture various aspects of departure of phylogenetic trees from the rake-shaped topology, which we hypothesized to be responsible for the PDXAR-EAR discrepancy (Fig. 3 and Supplementary Fig. 2). Predictors iv–vi describe some basic geometrical properties of species geographic ranges in the regions such as range shape, contiguity and complexity of the range edge, which have potential links to scaling patterns of biodiversity26,32. Predictor vii represents the size of the species pool (regional diversity), which was also suggested to affect beta diversity62 and the associated scaling of species richness. Predictors vii and ix describe autocorrelation structure of gradients of species richness, which we expect to broadly correspond with aggregation of species ranges within the domain, as illustrated by our simulations in the inset maps in Fig. 2a–d.

Models explaining AUC and ΔAUC by the regional predictors

We built three sets of statistical models: (i) Models that explain AUC of all of the EAR and PDXAR curves (162 data points). (ii) Models that explain ΔAUC of all pairs of EAR in and EAR out (27 data points). (iii) Models that explain ΔAUC of all pairs of EAR and PDXAR curves (81 data points). By ‘data points’, we mean unique combinations of taxa, regions and types of extinction curves. All models used the predictors described above as predictors and the taxon- and region-specific AUC or ΔAUC as a response. The specific predictors used in each of the three sets of models are listed in Supplementary Tables 2–4. For all models, we used ordinary least squares regression (normal error distribution). We first fitted models with all combinations of predictors (no interactions or nonlinear terms) and we ranked the models by their AIC c (Akaike Information Criterion with small sample size correction). For all predictors, we also calculated standardized regression coefficients (β) by rescaling each predictor to zero mean and variance of 1, and we averaged the betas over all of the models (Supplementary Tables 2–4).

All result from modelling, model selection and model averaging are summarized in Supplementary Tables 2–4. Models presented in Fig. 6 were selected among the models given in Supplementary Tables 2–4. When selecting these models we had in mind (i) their interpretability, (ii) out-of-sample predictive performance measured by AIC c , (iii) number of predictors, (iv) collinearity between the predictors (Supplementary Table 1) and (v) how well they corresponded with the averaged standardized coefficients (β) presented in the second line of Supplementary Tables 2–4. Finally, we note that our 27 combinations of taxa and regions used for the modelling give a relatively small sample size and are not independent realizations of a random process. This makes it problematic to calculate likelihoods and limits the interpretation of AIC c values and model rankings (Supplementary Tables 2–4). Consequently, we do not report P-values or standard errors.

We note that the standard, published procedures for dendrogram-based FD calculations are affected by the number, weighting and types of traits assessed44, and the clustering algorithm40,41. These differences may impact the exact shape of FDXAR curves. We therefore do not analyse these curves or FDXAR and PDXAR differences in detail and instead focus on the qualitative comparison between FDXAR and EAR curves, which is not affected by these methodological effects.

All of the distributional data, phylogenies, functional dendrograms and shapefiles used for the analyses described above are also provided (see Additional Information). Further technical details are in Supplementary Methods.