Using the IUCN Red List to map threats to terrestrial vertebrates at global scale
Species-level dataSpecies range maps were derived from BirdLife International and NatureServe50 and the IUCN51. The threat data were from the IUCN Threats Classification Scheme (Version 3.2), which contains 11 primary threat classes and almost 50 subclasses52. 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, 90%) and the severity, describing the scale of the impact on population declines: unknown, no decline, negligible declines, fluctuations, slow but significant declines (30%).Model development approachWe 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 mapsFirst, 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 km2) 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−6, which yields very strong autocorrelation; 1 × 10−4, 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:$${mathbf{Z}}(r) = U^{mathrm{T}}{mathrm{Norm}}left( {n,0,1} right)$$
(1)
$$W = UU^{ast}$$
(2)
$$W = {mathrm{e}}^{left( { – rD} right)}$$
(3)
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, UT 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 processSecond, 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’, CUncertain, 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 (C0.25, C0.5 and C0.75) or uncertain (CUncertain-0.25, CUncertain-0.5 and CUncertain-0.75) spatial knowledge.Model formulation and selectionThird, 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, PTh—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 (C0.25, C0.5, C0.75, CUncertain-0.25, CUncertain-0.5 and CUncertain-0.75), and the root mean squared error (RMSE) was calculated between PTh 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 r2 > 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 r2 between 0.47 and 0.7, depending on simulation parameters for Z(r) in 0.05, 10−4 and 10−6, 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 impactOnce 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, PTh, 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 List18. 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 km2) 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, PTh,Activity (for example, PTh,Logging for logging, PTh,Agriculture for agriculture or PTh,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 Services4, but our methodological framework is flexible and could be expanded to other threats in the IUCN classification19. 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 probabilityWe 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 patternsWe 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, PTh,Forest-loss) with data on forest cover change10. Forest cover change was aggregated from their native 30 m × 30 m (900 m2) 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 classification53. 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: PTh,Forest-loss increased with increasing forest cover loss (P More
