Quantification of habitat fragmentation reveals extinction risk in terrestrial mammals

As the basis of our fragmentation models, we used habitat-suitability models developed by Rondinini et al. (10) for the world’s extant terrestrial mammals (n = 5,027 species). The completion of the Global Mammal Assessment by the IUCN Species Survival Commission (IUCN/SSC) (2) provided an unprecedented opportunity to develop global-habitat models for all terrestrial mammals. The Global Mammal Assessment dataset, based on the taxonomy published by Wilson and Reeder (45), received the input of thousands of mammal experts belonging to more than 30 Specialist Groups of the IUCN/SSC. The dataset, which contains the known geographic range (i.e., extent of occurrence as defined by IUCN), species-specific qualitative textual descriptions of habitat preferences, and conservation status of each mammal species, is available through the IUCN Red List of Threatened Species (1).

From this database, the habitat-suitability models were developed at a 300-m resolution and limited to occur within the geographic range of the species (10). Three environmental variables formed the basis of the habitat models: elevation, type of land cover, and hydrological features. Elevational range where a species is found, when known and recorded in the IUCN Red List, was incorporated into the habitat models. Expert information on other habitat affinities, including preferred types of land cover, tolerance to human impact, and relationship to water bodies, were extracted from textual description within the IUCN database and input as quantitative data into the habitat models. Habitat models ranked areas with a three-level gradient of habitat suitability: (i) high, corresponding to primary habitat (i.e., preferred habitat where the species can persist); (ii) medium, corresponding to habitat where a species can occur, but not persist without nearby primary (i.e., high-suitability) habitat; and (iii) unsuitable, where a species is expected to seldom or never be found (10). A subset of models and their associated habitat-suitability ranks were validated against available points of known species occurrences. Habitat was further classified as to if it occurred within protected areas, using IUCN categories I–IV from the World Database of Protected Areas (46, 47). Full details regarding the development of these habitat models are available elsewhere (10), and data are available upon request from the model developers (https://globalmammal.org/activities/research/distribution-modelling/).

We measured fragmentation of high-suitability habitat, because such habitat is defined as essential for species persistence. When delineating high-suitability patches in the geographic range, we eliminated small patches (<four adjacent cells at 300-m resolution) potentially created by artifacts contained in underlying land-use and cover maps. This approach reduced the influence of spurious, isolated patches and improved computational efficiency of our fragmentation analyses. We then defined high-suitability habitat patches to be formed as clusters of suitable cells that were adjacent in any of the eight-neighborhood cells.

For each species we quantified the degree of fragmentation of high-suitability habitat by calculating the average Euclidean distance of all cells within high-suitability habitat from the nearest edge [that is, “GISfrag” (11, 12)]; edges demarcated the boundary between high-suitability and medium-suitability or unsuitable habitat, distributed either in the matrix external to habitat patches or as internal perforations within a patch. Low values of the average Euclidean distance into habitat for each species represent more highly fragmented habitat, whereas high values represent more core habitat and less-fragmented habitat. We considered a variety of other landscape metrics (including FRAGSTATS) and selected the GISfrag metric because it does not require an arbitrary distance threshold of what constitutes an “edge,” accounts for different shapes of patches and landscapes patterns and arrangements, is a robust measure that accounts for the distribution of patch area (48), is comparable across landscapes of different extents, and provides stable, readily interpretable information (12). Furthermore, Euclidean distance-to-edge was found to be singularly valuable in quantifying global forest fragmentation (8). We did not calculate distance between patches through the intervening landscape “matrix,” so our metric does not report on patch isolation or landscape configuration. For graphical display (Figs. 1 and 2 and Fig. S1), we inverse-coded the GISfrag metric to facilitate a more intuitive interpretation, where high GISfrag values represented high degrees of fragmentation.

Our analysis of fragmented habitat is also robust to any specific scaling assumption. Although information on dispersal ability can help inform the scale of fragmentation effects (49, 50), a comprehensive database of dispersal distance for all mammal species does not exist. Prior studies have attempted to approximate dispersal using allometric relationships based on body mass alone (51, 52). However, dispersal distances are highly heterogeneous even within the same species and are influenced by a variety of life-history traits and ecological factors, such as diet, sociality, and home-range area; thus, body size can be an inaccurate predictor of dispersal (50). Additionally, our inclusion of body mass as a predictor of extinction risk (see next section) further accounts for the potential latent influence of scale-dependent effects on species-specific responses to fragmentation, and inclusion of dispersal estimates would potentially lead to circularity and redundancy with our body size covariate. Rather than speculate about the effects of spatial scale on dispersal distance or other species-specific responses to fragmentation, our primary assumption is that further distance into the interior of a habitat patch is beneficial, given that it signifies larger patches and fewer edge and fragmentation effects [see also Haddad et al. (8) for a similar assumption to assess global forest fragmentation].

We used phylogenetic comparative methods to examine the relationship between the fragmentation metric and various species characteristics, including IUCN Red List status, range size, body mass, and the proportion of high-suitability habitat within the geographic range. We first used Pagel’s λ (53) to determine the strength of the phylogenetic signal in these variables. We found that the phylogenetic signal in the fragmentation metric was moderate (λ = 0.60) and similar to the values we obtained for range size (λ = 0.56), the proportion of high-suitability habitat within the range (λ = 0.45), and IUCN Red List status (λ = 0.67). Pagel’s λ value, however, was much higher for body size (λ = 0.99), a variable for which a strong phylogenetic signal was expected given that it is an evolved trait intrinsic to organisms. We therefore corrected for phylogenetic signal in our analyses. Because of nonnormal distributions, we ln-transformed the fragmentation metric, range size, and body mass, and arcsine square-root–transformed the proportion of high-suitability habitat in all statistical analyses.

We used the extant mammalian supertree phylogeny developed by Bininda-Emonds et al. (54, 55) that used the best estimates of divergence times (n = 4,510). To link this phylogeny to the taxonomy used by Rondinini et al. (10), we edited the two datasets to reconcile synonym species using Wilson and Reeder (45), and deleted species that were not common to both datasets. We did not incorporate new species into either dataset, and pruned all nonterrestrial mammals (cetaceans, seals, and sirenians) and Homo sapiens from the phylogeny, leaving a total of 4,018 species from 26 taxonomic orders. All phylogenetic analyses were conducted in the R statistical package (v3.1.3, R Development Core Team, 2015) using the ape (56), Geiger (57), and phylolm (58) packages.

Data on body mass were extracted from the PanTHERIA database (59). Where no body mass was listed, we first checked if it was present in the MOM database (60). If it was not, we followed the same methods as in Crooks et al. (12) and used the midpoint of the body mass range presented in Nowak (61) for that species. If not listed in Nowak (61), we used data from the closest member of the same species group, sister species of the same genus, or closest relative as indicated by Wilson and Reeder (45). If no other data were available, an average of all species in the genus was used. In other cases where PanTHERIA data were missing, information was copied and pasted from a synonym or sister species.

We used a phylogenetic generalized linear model (58) to evaluate the relative capacity of our fragmentation metric and two widely used macroecological variables, body mass and range size, to explain extinction risk, as indexed by IUCN status categories (13). Data Deficient species (n = 388) were initially excluded from these analyses given they have inadequate information to assess extinction risk (13), so we evaluated five status categories, including Least Concern (n = 2,633 for phylogenetic analyses), Near Threatened (n = 266), Vulnerable (n = 345), Endangered (n = 315), and Critically Endangered (n = 71). A species is classified by the IUCN as Least Concern when it is widespread and abundant and does not qualify for listing in other categories, and as Near Threatened when it does not qualify for a threatened category now but is close to qualifying for a threatened category in the near future (13). A species is classified as Vulnerable, Endangered, or Critically Endangered when the best available evidence indicates it faces a high, very high, or extremely high risk of extinction in the wild, respectively. Given these criteria, we reclassified the IUCN categories into a binary threat using two classification schemes. In the first analysis, we assumed species classified by the IUCN as Least Concern and Near Threatened faced no immediate extinction risk, whereas species classified as “threatened” (i.e., Vulnerable, Endangered, and Critically Endangered) faced some degree of extinction risk. This is a standard binary classification scheme in prior extinction-risk studies (16). We then conducted a second analysis, again categorizing Least Concern species as having no immediate extinction risk, but assuming Near Threatened species face some extinction risk, given that the IUCN identifies such species as close to being threatened; we thus view this comparison as a conservative and precautionary approach regarding extinction risk of Near Threatened species (13). We used AIC (62, 63) to compare models with all combinations of fragmentation, body mass, and range size as predictors of extinction risk; AIC analyses included calculations of model probabilities (the likelihood of the model given the data and model set) and evidence ratios (the weight of evidence of each model relative to the best model). Finally, to further explore patterns of fragmentation in species classified as nonthreatened by the IUCN, we also used phylogenetic generalized linear models to compare species listed as Least Concern to those listed as Near Threatened or Data Deficient.

The IUCN Red List process is expert-driven, and the IUCN uses five criteria to assign species to categories of extinction risk (13). One of those criteria is criterion B, used to list species that have restricted geographic ranges, with subcriteria for species that are severely fragmented or exist in few locations (subcriteria B1a/B2a), or are undergoing continuing decline (B1b/B2b) or extreme fluctuations in population size or distribution (B1c/B2c). To evaluate the reliability of criteria used by the IUCN to assess if a taxon belongs in the threatened category, in particular those subcriteria (B1a/B2a) relating to habitat fragmentation, we used a phylogenetic generalized linear model to compare fragmentation levels for threatened species classified under subcriteria B1a/B2a to fragmentation levels for threatened species that did not meet the B1a/B2a subcriteria. As expected, threatened species classified under subcriteria B1a/B2a had significantly higher degrees of fragmentation (Results and Discussion).

These results, however, might imply that our use of habitat fragmentation as a predictor of extinction risk may lead to potential circularity with the IUCN Red List criteria. To explore this issue, we followed prior approaches (25–27) and conducted additional analyses after first excluding threatened species listed under relevant listing criteria. Red List assessment guidelines stipulate that assessors should evaluate each species against all listing criteria, and all relevant criteria should be applied (13). Furthermore, the Red List guidelines encourage assessors to adopt a precautionary attitude to uncertainty when applying the criteria, so if range restriction or fragmentation was deemed a plausible threat by assessors, criterion B should be applied, even in the presence of uncertainty. Because of the relationship between range size and fragmentation (Results and Discussion, Fig. 2, and Fig. S1), we decided to adopt a conservative approach to account for potential circularity and exclude all species categorized because of range restriction (i.e., under any combination of subcriteria of criterion B, representing 46.8% of threatened species). We then used phylogenetic generalized linear models and AIC to compare models with all combinations of fragmentation, body mass, and range size as predictors of binary extinction risk. As above, we used two classification schemes to model extinction risk, assuming Near Threatened species (i) face no immediate extinction risk or (ii) face some extinction risk. For each model, we compared Least Concern/Near Threatened species to those threatened species not listed under criterion B.

We also evaluated how our fragmentation metric related to body mass and range size. We first conducted simple correlation analyses, corrected for phylogeny, relating fragmentation to body mass and also to range size. We then evaluated these relationships further with a phylogenetic generalized least-squares regression (64) that assessed the relative influence of body mass and range size on the fragmentation metric; we used AIC to compare models with all combinations of body mass, range size, and their interaction as predictors of habitat fragmentation. For our generalized least-squares regressions, we assumed that the continuous covariates evolved randomly after a speciation event (65), so we used a Brownian-motion model to define the structure of the correlation among species.

In addition to the macroecological variables, we also conducted two additional phylogenetic generalized least-squares regressions relating fragmentation to proportion of high-suitability habitat within the range and to proportion of protected high-suitability habitat within the range. Prior work has demonstrated that prevalence of remaining suitable habitat within the geographic range, or conversely the proportion of unsuitable human-affected land cover, is predictive of threat status in mammals (10, 14, 16) and other taxa (66, 67). We also conducted separate phylogenetic generalized linear models to evaluate the relative importance of fragmentation versus proportion of suitable habitat as sole predictors of extinction risk. These tests thus helped distinguish between the effects of habitat loss per se (as assessed by proportion of suitable habitat remaining in the range) from how remaining suitable habitat was distributed into intact, core habitat patches (as assessed by our fragmentation metric).

We generated global maps of fragmentation by summing the fragmentation metric at each 300 × 300-m cell for all terrestrial mammal species worldwide. Because we summed metrics for all species within a given area (Fig. 3), degree of fragmentation in part reflects global patterns of species distribution and richness. Thus, we also generated global maps displaying an average fragmentation index that was corrected for species richness (for details, see Supporting Information).

Our fragmentation analyses are based on the extent and geometry of high-suitability habitat and do not distinguish as to whether fragmentation is caused by natural factors (e.g., elevational or hydrological barriers) or anthropogenic disturbances (e.g., human development). To further identify the sources of fragmentation, we generated another map that highlighted anthropogenic fragmentation specifically by weighting the global fragmentation map with a recently constructed human modification (HM) model, which combines the effects of multiple stressors (e.g., urban and agricultural land cover, energy production, nighttime lights, and roads) into an overall score of HM globally [approach and technical details in Theobald (29) and Supporting Information]. When weighting our original fragmentation map with the HM model, the resulting map identifies regions that have been fragmented by human development specifically, and de-emphasizes regions that are naturally fragmented such as high-elevation areas and landscapes with water bodies interspersed (e.g., relictual glaciated areas).

Our methodology also allowed us to develop spatial models predicting the degree of habitat fragmentation globally for the world’s mammals. Degree of fragmentation was spatially quantified by summing the fragmentation metric at each 300 × 300-m cell for all terrestrial mammal species worldwide. When generating global maps, we adopted the finest grain possible (300 m), given that the majority (∼58%) of all terrestrial mammals are nonvolant species with a body size below 1 kg, so likely to respond to fragmentation at that grain.

Because we summed metrics for all species within a given area (Fig. 3), degree of fragmentation in part reflects global patterns of species distribution and richness. To further evaluate the potential influence of mammalian diversity on our spatial fragmentation models, we also generated a species-richness map (Fig. S2A) by summing the number of terrestrial mammal species with high-suitability habitat (as depicted with our GISfrag metric) for each 300 × 300-m cell globally. We note that these estimates of species richness are based not on the total extent of the geographic ranges of mammals, but rather on the extent of suitable habitat within their ranges; these distribution data predict known species occurrences better than geographic ranges (10). We then divided our fragmentation metric within each cell globally by the number of species with suitable habitat within that cell. Dividing by the number of species generated an average fragmentation index that was corrected for species richness. For comparative purposes, we then mapped this corrected fragmentation index globally (Fig. S2B).

To further investigate the sources of fragmentation, we generated a spatial model that highlighted anthropogenic fragmentation specifically by weighting the global fragmentation map with a recently constructed HM model (Fig. S3) [approach and technical details in Theobald (29)]. The HM model provides a comprehensive but parsimonious approach that combines the effects of multiple stressors into an overall score of HM globally. Specifically, the HM model compiles data on urban land cover from the Global Human Settlement Layer [68] (300 m resolution) and on urban and agricultural land cover from European Space Agency’s L4 LCCS dataset for 2010 (69) (300-m resolution) using intensity values from Theobald (29). For energy production, mining, and other miscellaneous human land activities, the HM model compiles data on “nighttime lights” from the 2013 National Oceanic and Atmospheric Agency (NOAA) VIIRS dataset (70) (400-m resolution, maximum normalized). Finally, the HM model maps roads from the 2010 gRoads dataset compiled by Center for International Earth Science Information Network (CIESIN; converted to 300-m resolution reflecting the proportion of a cell covered by the physical footprint of the road) (71). The HM model assumes the contribution of a given threat decreases as values from other threats overlap. Locations with multiple threats have a higher HM value than locations with just a single threat, but the cumulative HM score converges to 1.0 as multiple human-impact data layers are added. Individual factors were combined across multiple data layers using an increasive function (29), technically called a fuzzy sum (72), to create a 300-m resolution map. When weighting our original fragmentation map (Fig. S3A) with the HM model (Fig. S3B), the resulting map identifies regions (Fig. S3C) that have been fragmented by human development specifically, and de-emphasizes regions that are naturally fragmented such as high-elevation areas and landscapes with water bodies interspersed (e.g., relictual glaciated areas).