DOI: 10.5061/dryad.m905qfv8g
title: Data and R code from: Rethinking global hotspots for threatened terrestrial vertebrates

description: Aim We aimed to delimit hotspots for terrestrial threatened vertebrate
 species (HTV) through novel macroecological and statistical approaches.
 Location Global. Time Period Present day (1979–2024). Major Taxa Studied
 Terrestrial threatened vertebrate species (n = 7188). Methods In
 comparison with previous delimitations of hotspots, we: (i) considered
 richness and degree of endemism together through a robust statistical
 framework; (ii) focused on a priority set of species extremely important
 in terms of conservation, based on IUCN threat status; and (iii) used a
 fine spatial scale which allowed us to define key sub-areas within classic
 hotspots. We also assessed the degree of protection and human impact
 within the proposed HTV. Results We propose 20 global hotspots for
 threatened terrestrial vertebrates. In comparison with classic hotspots,
 proposed HTV have a significantly more limited distribution, covering ~27%
 of classic hotspots' area. In addition, a large proportion of HTV
 (~27%) does not match with classic hotspots. The overlap between HTV and
 protected areas (PAs) is low (&lt; 11%), and extremely low when only
 strict protected areas are considered (&lt; 1.5%). Also, a great
 degree of HTV exhibits high to extreme levels of human modification. On
 average, the velocity of climate change within HTV has been low, but
 attention must be given to notable areas presenting medium to high
 velocities. Interestingly, the geographical locations of highly endemic
 and rich areas considerably varied across individual vertebrate taxa. Yet,
 a high proportion of these priority areas for individual taxa are covered
 by the proposed HTV (74%–89%). Main Conclusions Our findings present key
 areas of the world for threatened terrestrial vertebrate species, many of
 these at high risk due to an interplay among low levels of protection,
 extreme levels of human modification and climate change. The proposed HTV
 are highly relevant in terms of decision-making, serving as a guide for
 allocating the limited conservation resources.
description: <strong>Species</strong> We
 obtained digital range maps (extent of occurrence maps) for amphibians,
 birds, mammals and reptiles from the IUCN database (IUCN 2024). Threatened
 species were defined as all terrestrial vertebrate species in the IUCN
 categories vulnerable (VU), endangered (EN) and critically endangered (CR)
 (retrieved from the IUCN webpage using the available filters in June
 2024). We discarded strictly marine species and included only the
 terrestrial portion of the distribution range of marine–terrestrial
 species. Polygons were converted to binary raster layers (0 = absence of
 species; 1 = presence of species) to a spatial resolution of 0.5°. A given
 raster cell was classified with a 1 if a polygon of a species overlapped
 with it to any extent. Next, we excluded the most widely distributed
 species, given the sensitivity of the range-diversity analysis to the
 inclusion of cosmopolitan species. To do this, we calculated the total
 number of cells covered by the distribution of each species, and discarded
 those species whose number of cells was above the 99th percentile
 considering all threatened species together. The subset of species
 resulting from these filters was used for all analysis, both considering
 all species together and for each individual taxon. Overall, 2869
 amphibians, 1294 birds, 1296 mammals and 1729 reptile species were
 considered after applying the mentioned filters, totalling 7188
 vulnerable, endangered or critically endangered species.
   <strong>Range-diversity
 analysis</strong> The range-diversity analysis is
 a statistical method that uses biodiversity presence-absence matrices to
 jointly describe the community composition of every location (i.e., cells
 of a grid or raster layer) within the analysed extent of analysis (Arita
 et al. 2008; Borregaard and Rahbek 2010; Soberón and Ceballos 2011;
 Soberón, Cobos, and Nuñez-Penichet 2021). To do this, the number of
 species and the mean dispersion field are calculated for each cell in the
 grid. The dispersion field is simply the total range size (i.e., number of
 cells) of all the species occupying a given cell and is considered an
 inverse measure of endemism (Soberón, Cobos, and Nuñez-Penichet 2021;
 Osorio-Olvera and Soberón 2023). By applying appropriate statistical
 thresholds to the dispersion field values (outside the bounds of 95% CI of
 random expectations; see Appendix 1 in Supporting Information), which are
 related to the level of endemism in each location, each cell was
 classified as exhibiting low or high levels of endemism (LE or HE,
 respectively) (Soberón, Cobos, and Nuñez-Penichet 2021). Additionally,
 each cell was classified according to the species richness: low richness
 (LR), less than 25th percentile; intermediate richness (IR), between 25th
 and 75th percentile; and high richness (HR): more than the 75th
 percentile. Combinations of levels of endemism and richness determined six
 distinct categories: HE/LR, HE/IR, HE/HR, LE/LR, LE/IR and LE/HR. Also, a
 cell was classified as ‘random’ if its level of endemism was not
 statistically significant according to its dispersion field value. We
 performed rangediversity analyses for all vertebrate species, and
 separately for each evaluated taxonomic group (amphibians, birds, mammals
 and reptiles). Previous works show that patterns of endemism can be
 scale-dependent (Daru et al. 2020; Hurlbert and Jetz 2007). Thus, as an
 additional analysis, we repeated all range-diversity analyses at 1°
 resolution. Overall, results are very similar to those obtained at 0.5°
 resolution (see Figures S1.10–S1.14 and Table S1.4 in Appendix 1 in
 Supporting Information). Highly endemic and rich areas
 (hereafter HE/HR areas) were classified into distinct hotspots for
 threatened vertebrate species (HTV), based on the level of aggregation and
 clumpiness of individual HE/HR areas (i.e., cells). The names for HTV were
 selected based on their geographical location and their overlap with
 classic biodiversity hotspots and ecoregions. Moreover, we calculated the
 representativeness of HE/HR areas (considering all threatened vertebrate
 species together) within biomes, as well as the percentage of the areas of
 each biome that is occupied by these HE/HR areas. Polygons of biomes were
 obtained from Olson et al. (2001). Finally, HTV were compared with the
 areas occupied by the latest biodiversity classic hotspot's
 revisitation (totalising 36 hotspots) (Mittermeier et al. 2011; Noss
 2016). Polygons of hotspots were downloaded from http://www.cepf.net   <strong>Protected area analysis</strong> PAs polygons were downloaded from the World Database of Protected Areas (UNEP-WCMC and IUCN 2022; accessed in November 2022) by selecting all PAs under the status of ‘adopted’, ‘established’ or ‘designated’, including the six categories of PAs, as well as PAs classified as ‘Not Applicable’, ‘Not Assigned’ and ‘Not Reported’. As the information about PAs of some countries is currently publicly missing, the database was completed with other sources. Protected areas from China were obtained from a previous database from the World Database of Protected Areas of 2016, while PAs from India, Turkey and Estonia were obtained from OpenStreetMap services (OpenStreetMap contributors 2024). We classified the PAs as PAs with strict conservation goals (categories I to IV) and PAs with non-strict conservation goals (categories V and VI, as well as those classified as ‘Not Applicable’, ‘Not Assigned’ and ‘Not Reported’). Next, we generated three distinct raster files with a spatial resolution of 0.0042° (~0.5 km), one including all PAs, another including strict PAs and another including non-strict PAs. The distribution of the various categories of endemism/richness (the percentage of the assessed area that is covered by each category) was calculated for the area occupied by PAs and compared with the distribution of the categories of endemism/ richness of all cells (including protected and unprotected cells: global). Complementarily, we calculated the Euclidean distance between the distribution of percentages inside PAs and the distribution of percentages that we obtained for cells of any kind (either protected or not protected). The calculated Euclidean distance serves as a measure of dissimilarity between distributions: the higher the value, the greater the differences between the distribution of values inside PAs (all PAs, strict PAs or non-strict PAs) and the distribution of values for cells of any kind. Lastly, to detect how well protected each category of endemism/richness is, we calculated the percentage of the area classified into each category of endemism/richness that fell inside a PA of any kind, and by a PA with strict and non-strict conservation goals. These analyses were performed for all threatened species together, and for each taxonomic group separately.   <strong>Human Modification Index Analysis</strong> To estimate the degree of human modification of each category of endemism/richness, we used the latest estimation of the global human modification index of terrestrial systems (HMI) (Kennedy et al. 2019). The HMI index provides a cumulative measure of human modification of terrestrial lands based on modelling 13 anthropogenic stressors and their estimated impacts using spatially explicit global datasets at a resolution of 1 km2. The score of HMI ranges from 0 (no human modification) to 1 (maximum human modification) for each 1 km2 pixel. Following Watson et al. (2016), Jones et al. (2018), and Nori et al. (2022), we determined a minimum HMI value of 0.2 as a high level of human modification, which represents the mean value on pasture lands as a reasonable threshold when species are likely to be threatened by habitat conversion; and a minimum value of 0.4 as an extreme level of human modification, which represents the average HMI in the 10 grassland ecoregions with the highest average values for this index. Low levels of human modification were defined as areas with values between 0 and 0.1 (Kennedy et al. 2019), and moderate levels as areas with values between 0.1 and 0.2. We overlapped the raster file with HMI values with the raster containing the areas occupied by each category of endemism/richness obtained from the range-diversity analysis that included all species together. Mean, median and standard deviation of HMI values were calculated for all the cells classified within each category of endemism/richness.   <strong>Velocity of climate change analysis</strong> The velocity of climate change estimates the velocity at which a species would need to move to maintain constant climatic conditions (Loarie et al. 2009). We calculated a gradient-based climate change velocity index, which is calculated as the ratio between the long-term temporal trend in the climate conditions (e.g., temperature) by the spatial gradient in the climate conditions, and its units are in kilometres per unit of time (e.g., years or decades). High values of climate change velocity indicate that a species in a location should move fast and/or far to encounter a new location with an analogous climate condition. To estimate the velocity of climate change, we used temporal trends of mean monthly air temperature at 2 m for the 1979–2019 period, provided by CHELSA (Karger et al. 2017) at a resolution of 0.0083° (~1 km). Next, we calculated a mean annual temperature, which was the variable we used to calculate the gradient-based climate change velocity index. The resulting worldwide raster of velocity of climate change was then overlapped with the raster containing the areas occupied by each category of endemism/richness obtained from the range-diversity analysis that included all species together. Mean, median and standard deviation of climate change velocity values were calculated for all the cells within each category of endemism/richness. All analyses were performed in the R environment (R Core Team 2023) by using base R packages and packages <em>raster </em>(Hijmans 2024a), <em>terra </em>(Hijmans 2024b), <em>sp </em>(Bivand et al. 2008), <em>sf </em>(Pebesma and Bivand 2023) and <em>fasterize </em>(Ross 2020). Rangediversity analyses were performed with R package <em>bamm </em>(Osorio-Olvera and Soberón 2023). The velocity of climate change was calculated with the R package <em>VoCC </em>(García Molinos et al. 2019). All maps were created with QGIS in a Mollweide projection (QGIS Development Team 2021). See Appendix 1 in Supporting Information for more details of the applied methodology. A list of the included species for our analyses is available as a Supporting Information.
description: # Data and R code from: Rethinking global hotspots for threatened
 terrestrial vertebrates ## Description of the data and file structure ###
 R\_scripts\_and\_data.zip Contains all used datasets and raster layers
 used for statistical analyses. It also includes the R code scripts to
 perform these analyses. R scripts: * range-diversity analysis.R: main
 analysis that classifies cells according to the level of endemism and
 species richness. * hotspots analysis.R: analysis that overlaps proposed
 hotspots with classic hotspots. * biome_analysis.R: analysis that overlaps
 proposed hotspots with global biomes. * PAs analysis.R: analysis that
 determinas the level of protection by PAs of the proposed hotspots. * HMI
 analysis.R: analysis that determines the level of human modification
 within proposed hotspots. * velocity of climate change analysis.R:
 analysis that determines the velocity of climate change within proposed
 hotspots. Directories: Note: file names containing the terms
 "ANF",  "AVE", "MAM", and "REP"
 correspond to data on amphibian, bird, mammal, and reptile species,
 respectively. * pams: presence-absences matrices (.csv) for each taxonomic
 group and for all species, at 0.5 degree resolution. These matrices
 depcits the presence (1) or absence (0) of each species in each
 geographical location (raster cell).  * pams_1-degree: presence-absences
 matrices (.csv) for each taxonomic group and for all species, at 1-degree
 resolution. These matrices depcits the presence (1) or absence (0) of each
 species in each geographical location (raster cell). * rdivan_objects: RDS
 files (.RDS) outputs of range-diversity analysis at 0.5 degree resolution.
 * rdivan_objects_1-degree: RDS files outputs (.RDS) of range-diversity
 analysis at 1 degree resolution. * spp_distr_rasters: raster layers (.tif)
 of species distribution at 0.5 degree resolution. *
 spp_distr_rasters_1-degree: raster layers (.tif) of species distribution
 at 1 degree resolution. ### List of species.xlsx Table depicting all
 species considered for the analysis: * Column "Taxon": the taxon
 of each species. * Column "Species": the scientific name of the
 species. ### HTV.zip Contains the vector data for the proposed hotspots
 for terrestrial threatened species. Files and directories: * HTV.gpkg:
 vector file of proposed hotspots. The attribute tables cointains an id
 (first column) and the name of each hotspot (second column). * Directory
 "Shapefile": contains vector files of proposed hotspots. The
 information is identical to "HTV.gpkg" but in shapefile format.
 * world_land.gpkg: vector file of world land. * map_htv_QGIS_PROJECT.qgz:
 Quantum GIS project that relates all mentioned vector files. * README.pdf:
 help file that explains the content of the zip file. # Usage .R files are
 R scripts that can be accessed with R or RStudio. .RDS files are R
 outputs, objects created in the R environment, and can be loaded to an R
 session with the base function readRDS(). The R package "bamm"
 must be installed and loaded in the session to properly view these
 specific objects. .tif, .gpkg and .shp files (and associated files within
 the directory "Shapefile" inside HTV.zip) can be accessed with
 any GIS software, such as Quantum GIS. **Please read the published paper
 for more details and context.**