Using the IUCN Red List to map threats to terrestrial vertebrates at global scale

Species-level data

Species range maps were derived from BirdLife International and NatureServe⁵⁰ and the IUCN⁵¹. The threat data were from the IUCN Threats Classification Scheme (Version 3.2), which contains 11 primary threat classes and almost 50 subclasses⁵². In Red List assessments, assessors assign those threats that impact the species. For birds, the scope of the impact is also recorded categorically as the percentage of the species population that the threat impacts (unknown, negligible, <50%, 50–90% or >90%) and the severity, describing the scale of the impact on population declines: unknown, no decline, negligible declines, fluctuations, slow but significant declines (<20% over ten years or three generations), rapid declines (20–30%) or very rapid declines (>30%).

Model development approach

We designed our analytical framework with three considerations in mind. First, the threat location information is limited: for each species, the data only describe whether a species is threatened by a given activity anywhere within its range (data on the timing, scope and severity of threats are available only for birds and are not spatially explicit). Second, we wanted to compare the spatial patterns of threat against independent data on spatial distributions of human activities. Third, for many activities, the relationship between human activity (for example, hunting or invasive species and diseases) and biodiversity response is poorly understood. We therefore chose not to incorporate known patterns of human activity as explanatory variables in our models.

In the absence of global datasets on the spatial patterns of the impact probability of each threat, we used a simulation approach to develop our models and assess the ability of different model parameterizations to reproduce our simulated threat. This process had four steps (Extended Data Fig. 1).

Simulated threat intensity maps

First, we simulated a continuous synthetic threat across sub-Saharan Africa. The concept behind this is that a credible model should be able to reproduce a ‘true’, synthetic threat pattern on the basis of information comparable to that available in the Red List. To test this, we generated a set of synthetic, continuous surfaces of threat intensity with different levels of spatial autocorrelation and random variation (Supplementary Fig. 1). This was achieved by taking a grid of 50 km × 50 km (2,500 km²) pixels across the Afrotropic biogeographic realm (i.e., sub-Saharan Africa). Threat intensity was modelled as a vector of random variables, Z, one for each pixel i, generated with a correlation structure given by the distance matrix between points weighted by a scalar value, r, indicating the degree of correlation (equations (1–3)). Four values of r were used: 1 × 10⁻⁶, which yields very strong autocorrelation; 1 × 10⁻⁴, which yields strong autocorrelation; 0.05, which yields moderate autocorrelation; and 0.3, which produces a low-correlation, localized pattern (Supplementary Fig. 1). The model included the following equations:

where r is a scalar determining the degree of spatial autocorrelation (as r decreases, the autocorrelation increases), D is the Euclidean distance matrix between each pair of pixels, W is the matrix of weights for the threat intensity, U and U^* are the upper triangular factors of the Choleski decomposition of W and its conjugate transpose, U^T is the transpose of U and n is the number of pixels.

We chose the Afrotropic biogeographic realm (sub-Saharan Africa) as our geography within which to develop the modelling approach because it permitted more rapid iterations than a global-scale simulation while also retaining characteristics of importance for the model evaluation such as strong environmental gradients and heterogeneity in species richness. However, for the simulation, no information from the geography or overlapping species ranges was used, except the spatial configuration of the polygons. Thus, the use of the Afrotropic realm was purely to avoid generating thousands of complex geometries for the purpose of the simulation. Using a real geography and actual species ranges ensures that our simulation contains conditions that are observed in reality (for example, areas of high and low species richness also observed in the real world). We took the simulated threat maps generated through this process to be our ‘true’ likelihood of a randomly drawn species that occurs in that location being impacted by the synthetic threat (Supplementary Fig. 1).

Simulating the red-listing process

Second, we wanted to simulate the red-listing process whereby experts evaluate whether a threat is impacting a species on the basis of the overall threat intensity within its range. For this, we used the range maps for all mammal species in Africa and assigned a binary threat classification (that is, affected or not affected) to each species on the basis of the values of the synthetic threat within each species’ range. We assumed that the binary assessment of threat for a species is based on whether the level of impact across a proportion of its range is judged as significant. This step was intended to replicate the real red-listing process, where assessors define threats that impact the species on the basis of an assessment of the information available on threatening mechanisms and species responses. In practice, this was done by overlaying the real range maps for mammals over the four simulated threat surfaces and assessing the intensity of synthetic threat within each species range map. We wanted to assign species impacts considering that species will be more likely to be impacted if a greater part of their range has a high threat intensity. Understanding how to set a threshold for what intensity would constitute sufficient threat to be assessed as affected is a complicated exercise. We thus tested three thresholds to capture different assumptions. These thresholds were chosen after discussion with leading experts on the red-listing process. More specifically, we calculated the 25th, 50th and 75th percentiles of threat intensity across pixels within the species range. We then used a stochastic test to convert these quantiles to binary threat class, C. For each species, we produced a set of ten draws from a uniform distribution bounded by 0 and 1. If over half of the draws were lower than the threat intensity quantile, the species was classified as threatened for that percentile.

The above simulation assumes perfect knowledge of the threat intensities across the species range, which might not always be the case in the actual red-listing process. In real life, certain areas within species ranges are less well known for a suite of different reasons. To incorporate some uncertainty about the knowledge of the red-listing experts about the ‘true’ threat intensity, we constructed a layer to describe the spatial data uncertainty associated with the Red List. This aspect was intended to simulate the imperfect knowledge of the simulated ‘Red List assessors’. This layer was calculated as the proportion of species present in a given location that are categorized as Data Deficient—in other words, there is insufficient information known about the species to assess its extinction risk using the IUCN Red List Criteria (Extended Data Fig. 7). Then, when calculating the 25th, 50th and 75th percentiles of threat intensity across each range, we weighted this calculation by one minus the proportion of Data Deficient species, so that more uncertain places (those with a greater proportion of Data Deficient species) contributed less to the calculation than locations where knowledge was more certain. These were then converted to a binary threat class accounting for uncertainty in expert knowledge among the simulated ‘assessors’, C_Uncertain, using the same stochastic process described above for the calculation of C.

This step produced, for each species, a threat classification analogous to the threat classification assigned by experts as part of the IUCN Red List process. Six sets of threat classifications were produced for each synthetic threat surface, on the basis of the 25th, 50th and 75th percentiles with perfect (C_0.25, C_0.5 and C_0.75) or uncertain (C_{Uncertain-0.25}, C_{Uncertain-0.5} and C_{Uncertain-0.75}) spatial knowledge.

Model formulation and selection

Third, using all species polygons with assigned threat assessments from step 2 (that is, affected or not affected), we fitted nine candidate models and predicted the estimated probability of impact for each grid cell. Then, in a fourth step, we compared the predicted probabilities of impact produced in step 3 with the original synthetic threat maps created in step 1 to test the predictive ability of our models.

The Red List threat assessment does not contain information on where in the range the impact occurs. Therefore, a species with a very small range provides higher spatial precision about the location of the impact, whereas a species with a large range may be impacted anywhere within a wide region. To address this lack of precision in the impact location, we took the area of each species range to serve as a proxy for the spatial certainty of the impact information. The certainty that a species was impacted or not impacted in a given cell depended on its range size, R. The models we evaluated therefore incorporated R in different ways (Supplementary Table 1).

The models were fitted as a binomial regression with a logit link function. For each pixel, the model predicts the probability of impact, P_Th—in other words, the probability that if you sampled a species at random from those that occur in that pixel, the species would be impacted by the activity being considered. To account for uncertainties in the simulation of the threat assessment process (thresholds for impact and perfect or imperfect knowledge), models were fitted to the six sets of threat codes (C_0.25, C_0.5, C_0.75, C_{Uncertain-0.25}, C_{Uncertain-0.5} and C_{Uncertain-0.75}), and the root mean squared error (RMSE) was calculated between P_Th and the simulated threat intensity, Z(r), for each value of r. For each simulation, we ranked the different models according to their model fit as measured by the RMSE. We assessed these ranks across all simulations and sets of threat codes. We evaluated the models on the basis of the ranks of RMSE, across the threat code sets and threat intensity maps. Rank distributions for each model are shown in Extended Data Fig. 2, and the results from these models are shown in Supplementary Tables 1 and 2.

All models were correlated (Pearson’s r² > 0.5), albeit with some variation between model types and across the simulation parameters (Supplementary Fig. 2). However, some models had greater predictive accuracy when evaluated using the RMSE. The top four ranking models were, in order of decreasing summed rank, (1) inverse of cube root of range size as a weight, (2) inverse 2.5 root of range size as a weight, (3) inverse square root of range size as a weight and (4) inverse natural logarithm of range size as a weight. The fact that these four models showed good model fit suggests that the best model structure had a measure of range size as a weight but that the model was not particularly sensitive to the transformation of range size.

The best-fitting model across the range of simulation parameters was an intercept-only logistic regression where the response variable was the binary threat code (1 = threatened, 0 = not threatened) for each species in the pixel and where the inverse cube root of the range size of each species was used as a weight. The model was concordant across the set of simulated datasets with a relationship that was predominantly linear with r² between 0.47 and 0.7, depending on simulation parameters for Z(r) in 0.05, 10⁻⁴ and 10⁻⁶, centred around unity and with the RMSE ranging between 0.129 and 0.337 depending on simulation parameters (Supplementary Figs. 2 and 3). The choice of the inverse cube root range size weight was based on the performance of this against eight other model types (Supplementary Fig. 4 and Supplementary Table 1).

We conducted a decomposition of variance in model performance using a binomial regression model, with RMSE as the dependent variable and model type, knowledge level and autocorrelation structure as the independent factorial variables. This showed that knowledge about the threats underlying each species range and how that threat information is used in the assessment explained the vast majority (94.7%) of the variation in RMSE outcomes (Supplementary Fig. 4).

For birds, further information on the scope of the threat was available as an ordinal variable describing the fraction of range that the threat covers. We explored the use of scope in our models but concluded that to avoid arbitrary decisions about the scope of non-threatened species (where they are either not threatened anywhere or threatened in only a small part of their range), and for consistency with other taxonomic groups, we would model birds using the same model structure as used for mammals and amphibians (see the Supplementary Methods for further details).

Mapping probability of impact

Once the best-performing model was identified using the simulated data, we then used this model on the actual Red List threat and range data to develop threat maps. This model produced threat maps for each taxonomic group (amphibians, birds and mammals) of the probability of impact, P_Th, for each individual threat. For a given pixel, threat and taxonomic group, this estimates the probability that a randomly sampled species with a range overlapping with that pixel is being impacted by the threat, while taking into account spatial imprecision in the Red List data.

Threat maps were generated using range map data and threat assessments from the IUCN Red List¹⁸. We intersected range maps for 22,898 extant terrestrial amphibians (n = 6,458), birds (n = 10,928; excluding the spatial areas within the range that are associated with ‘Passage’—where the species is known or thought very likely to occur regularly during relatively short periods of the year on migration) and mammals (n = 5,512; including those with uncertain ranges) with a global 50 km × 50 km (2,500 km²) resolution, equal-area grid for the terrestrial world. This provided, for each 50 km × 50 km pixel, a list of the species whose range overlapped it, along with the associated range size of each species. For each pixel and taxonomic group (amphibians, birds and mammals) independently, we then modelled the probability of impact, P_Th,Activity (for example, P_Th,Logging for logging, P_{Th,Agriculture} for agriculture or P_Th,Pollution for pollution), for each of the six threats: agriculture, hunting and trapping, logging, pollution, invasive species and diseases, and climate change. We focused on these as the six main threats as defined by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services⁴, but our methodological framework is flexible and could be expanded to other threats in the IUCN classification¹⁹. We used only taxonomic groups with a sufficiently high total number of species and where they have been comprehensively assessed so that potential biases associated with the groups of species prioritized by experts are avoided.

Calculating uncertainties for the threat probability

We estimated a measure of uncertainty associated with our impact probability predictions using maps of the proportions of Data Deficient species in each cell within each taxonomic class (amphibians, birds or mammals) as a measure of knowledge certainty in that cell. The rationale for this approach is that places with more Data Deficient species with unknown threatened status should have greater uncertainty in the probability of impact. We therefore created greater variation in the data where there were more Data Deficient species. We used the knowledge-certainty map to probabilistically draw a sample of 100 threat codes for each species, on the basis of the median Data Deficiency across the species range. The random sample changed the species threat code with a probability related to the proportion of Data Deficient species within its range. If the median proportion of Data Deficient species was zero, then we assumed that there was a small probability (0.005) that the species could have been incorrectly coded. Where the median proportion was greater than zero, the probability increased linearly. So, for a species with 5% Data Deficient species within its range, the sample changed the species threat code with a probability close to 5%; if the median proportion was equal to 0.5, then the probability of the species being incorrectly assigned was equal to 0.5. We then fitted the impact probability model with each of the 100 species threat codes and generated a distribution of predicted threat probabilities in each grid cell, from which we took the 95% confidence intervals as the uncertainty estimates (Extended Data Figs. 8–10).

Evaluating modelled threat patterns

We evaluated the spatial patterns of threat on the basis of the real Red List threat assessment data against empirical data in two independent ways. First, we compared the probability of impact from logging and agriculture combined within forested biomes (that is, corresponding to remotely detected forest loss, which we refer to as the probability of impact from forest loss, P_{Th,Forest-loss}) with data on forest cover change¹⁰. Forest cover change was aggregated from their native 30 m × 30 m (900 m²) resolution pixels to our 50 km × 50 km resolution pixels using Google Earth Engine. For each 50 km × 50 km pixel, we calculated the total area lost between 2000 and 2013 and the area lost as a proportion of the area in 2000. We restricted our analysis to forested biomes: (1) tropical and subtropical moist broadleaf forests, (2) tropical and subtropical dry broadleaf forests, (3) tropical and subtropical coniferous forests, (4) temperate broadleaf and mixed forests, (5) temperate coniferous forests and (6) boreal forests/taiga, following the World Wildlife Fund’s ecoregions classification⁵³. The relationship between forest loss and the probability of impact from forest loss as captured by agriculture and logging overall showed a significant positive correlation: P_{Th,Forest-loss} increased with increasing forest cover loss (P < 1 × 10⁻⁵, Supplementary Fig. 5) but also showed some nuances.

Second, we evaluated threat levels against threat for about 6,000 KBAs assessed by specialists through a rigorously tested and standardized approach developed by Bird Life International²⁰. For a given activity, we calculated the median of predicted impact probabilities within each KBA and then grouped these KBA estimates by KBA severity class. On average, P_Th was higher in KBAs assessed as having a high severity of threat from an activity than in KBAs classed as having low threat. Significant relationships (P < 0.05, Wilcoxon test) were found in one or more taxonomic groups for logging, agriculture, hunting and climate change. For both evaluations, we conclude that the modelled spatial patterns of threat were consistent with expectations from the empirical data (Supplementary Figs. 5–11). We subsequently shared threat maps with taxon-specific experts from the Red List assessment groups, who qualitatively reviewed the patterns and assessed them as consistent with expert knowledge. Further details on the evaluation methods can be found in the Supplementary Methods.

We suggest that the broad concordance with empirical data in two independent evaluations demonstrates that the maps of impact probability are sufficiently plausible to underpin the findings of this study. We provide a framework that can easily be updated with future versions of the IUCN data, and we also stress that our approach should be viewed as an initial attempt to map patterns of threat impacts, which should be used iteratively to guide the collection of new data and improvement of the mapping approaches used.

Comparing threat occurrence likelihoods and the HFI

HFI data for the year 2009 were taken from Venter et al.¹¹. The native resolution of the index was 1 km², so we calculated the mean HFI in each 50 km × 50 km pixel used in our analysis. The HFI was standardized to lie within the bounds 0 and 1 by dividing by the maximum HFI value (50).

We compared the standardized, averaged HFI values to the predicted likelihood of threat occurring from land use change, hunting and trapping, and climate change. For land use change, we combined agriculture and logging by calculating the unweighted mean threat occurrence likelihood across taxa for these two threats. For hunting and trapping, we took the mean threat occurrence likelihood across taxa, while for climate change we used the predicted threat occurrence likelihood for birds.

To plot the spatial relationship between HFI and mean threat occurrence likelihood, we used a two-dimensional kernel density estimator (MASS package⁵⁴) to estimate the variation in the density of pixels for a given HFI and mean threat occurrence likelihood.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com