Three-quarters of species’ ranges have not been covered by protected areas in global borders

Protected area data

Terrestrial protected area data were derived from the February 2023 version of the World Database on Protected Areas (WDPA, accessible on: https://www.protectedplanet.net/en)⁵². In total, we had location data for 253,674 protected areas. We followed previous global studies²⁴ and the WDPA’s recommendations (https://www.protectedplanet.net/c/calculating-protected-areacoverage) on cleaning data and used only the polygon boundary data layer. Point data layers, PAs without designated, inscribed, or established status, and the UNESCO Man and Biosphere Reserves were excluded from our analysis. We assigned PAs to each country based on the ISO3 code reported in the WDPA, which represents the country or territory in which the PA is located. For PAs that have more than one ISO3 code, we assigned the PA to each ISO3 code reported. The PA data were clipped to terrestrial areas, inland lakes and waterways by the polygons of Global Administrative Areas version 3.4 (GADM, 2020; https://www.gadm.org). For a few PAs lacking establishment time information, we randomly selected a year of establishment among PAs in the same country, following the approach used in previous studies²⁷. Finally, we dissolved the overlapping polygons and retained larger areas for polygons with earlier establishment. Although there is variation in number of PAs in the WDPA, we did not include additional data from country-level or regional databases because the WDPA is the only authoritative dataset following globally consistent standards and is regularly validated and updated to maintain the highest data quality^21,53.

Species distribution data

Data on the species distributions and global conservation status of terrestrial mammals (n = 5577) and amphibians (n = 7231) were obtained from the IUCN Red List of Threatened Species⁵⁴. Bird distribution data (n = 10,987) were sourced from BirdLife International and Handbook of the Birds of the World⁵⁵, and reptile data (n = 10,899) were collected from Roll et al.^56,57. These distributional maps represent the most comprehensive spatial databases for different taxonomic groups based on expert reviews, especially for spatial boundaries of the distributions of each species²⁴. Following Titley et al.³, we included only the parts of each species’ ranges where the species was considered to be “Extant” or “Probably Extant” (“Presence” codes 1 or 2), where the species was native (“Origin” code 1), and where the species had breeding or resident ranges (“Seasonality” code 1 or 2). We did not account for subspecies. Species are classified according to their Red List conservation status as Critically Endangered (CR), Endangered (EN), Vulnerable (VU), Near Threatened (NT), Least Concern (LC), or Data Deficient (DD). Species classified as either Extinct (EX) or Extinct in the Wild (EW) were excluded from the analysis. Species classified as VU, EN, or CR are classed as “threatened”.

Border area identification

We used data from the GADM and identified a total of 330 national terrestrial borders worldwide based on unique ISO3 codes, excluding coastlines. In addition, we use data from the Ministry of Natural Resources’ standard mapping service (http://bzdt.ch.mnr.gov.cn/; No. GS (2016) 1666) as a reference in matters related to territorial boundaries. Our analyses did not include the Antarctic areas, Australia and New Zealand as they do not share terrestrial neighboring borders. Following Thornton et al.¹¹, we created three distances namely 25 km, 50 km, and 100 km on both sides to border lines as buffer zones. Next, we overlaid these buffer zones with grid cells of 5 km × 5 km resolution to define candidate transboundary areas.

Protected area coverage across species ranges

We overlaid the terrestrial species polygons (n = 34,694) onto the grid cells identified as transboundary areas, and included species in our analysis if any part of their range overlapped with transboundary areas. Finally, we obtained a dataset consisting of 2943 amphibians, 4906 reptiles, 7710 birds, and 3480 mammals, which accounts for 54.9% of all the 34,694 terrestrial species. The average proportion of their distributional ranges overlapping the transboundary areas was ~30.39% (Supplementary Data 2). Specifically, 3251 (17.1%) species had ≥ 50%, and 1409 (7.4%) species had > 90% of their ranges in transboundary areas. For threatened species, these proportion were 41.8% and 26.2%, respectively. Overall, the threatened species have a greater proportion of distributional ranges (48.1%) than the non-threatened species (24.0%) in transboundary areas (Kruskal-Wallis test, χ² = 590.51, P < 0.001, Supplementary Fig. 5a).

Next, we overlaid the distributional boundaries of each species with the boundaries of PAs to identify proportion of their transboundary ranges protected. All geospatial data processing was carried out in the Mollweide equal-area projection using ESRI ArcGIS Pro version 3.0.1 (https://www.esri.com/).

For each species, we calculated the mean coverage C_i using the following formula (1):

where c_ij is the coverage by PAs of species i in cell j and n is the number of cells that overlap with species i’s ranges.

c_ij is calculated using the following formula (2):

where area_protected is the size of species transboundary range covered by PAs and area_range is the size of species’ transboundary range in the grid cell. For each gird cell, species range coverage by PAs was defined by evaluating the average of the transboundary ranges of all species in that grid that were covered by PAs.

To compare PA coverage of species ranges between transboundary areas and non-transboundary areas, we applied a resampling approach by randomly selecting 50,000 pairs of grid cells from transboundary areas and from non-transboundary areas, respectively, and this process was repeated 999 times⁵⁸. The grid cells were selected either randomly (Random), or with similar land areas (i.e., within a deviation of 10%) covered by PAs (PA availability), sampling intensities (GBIF), mean species richness (Species richness), socioeconomic conditions (Night time light) or topographic complexity (data source and references in Supplementary Table 13) to control for variations in the sampling efforts reflected by different anthropogenic factors.

Sensitivity analyses

We used buffer distances of 25 km and 100 km to test whether our results were sensitive to the buffer distance used to calculate the coverage of species ranges by PAs in transboundary areas, and obtained qualitatively similar results to the analyses using 50 km (Supplementary Fig. 7a–b). We also repeated the analysis using 10 km × 10 km and 25 km × 25 km grid cells and obtained similar results to the analyses conducted using 5 km × 5 km grid cells (Supplementary Fig. 7c–d). As our results were robust to different buffer classes and grid sizes, we show the results using the 5 km × 5 km grid cells with a buffer zone of 50 km in the main text and provide the results of additional analyses in the supplementary materials.

Ecological and socioeconomic variables related with PA coverage of species ranges

To explain spatial variations in PA coverage of species ranges across global borders, we selected 10 variables, namely, topographic complexity, governance, cropland expansion, human development index (HDI, which captures elements of life expectancy, education, and wealth at the country level), human population density, frequency of armed conflicts, bilateral cooperation abilities, size, years of establishment, and IUCN categories of PAs, which are considered potentially important variables that influence the PA coverage of species ranges in transboundary areas (see Supplementary Table 14 for rationales). We do not aim to indicate causal relationships between the outcome and the factors but instead aim to understand how the PA coverage of species ranges might vary across different socioeconomic and environmental contexts.

The difference between the maximum and minimum altitudes was used to measure the topographic complexity, because higher altitude differences can create a wide variety of habitats and environmental niches⁵⁹, and thus more biodiversity protection is expected. The ‘Zonal Statistics’ function in ESRI ArcGIS Pro was applied to extract the highest and lowest points of these altitudes from EarthEnv (https://www.earthenv.org/) at a 1 km resolution⁶⁰ in each grid cell polygon defined as transboundary areas. We included cropland expansion as an indicator of land-use needs using data from Potapov et al.¹⁷ which is a global time series cropland maps from 2000 to 2019 at a spatial resolution of 30 m (available at https://glad.umd.edu/dataset/croplands). Human development index (HDI) values were derived from the 2022 Human Development Report published by the United Nations of Development Programme (https://hdr.undp.org/content/human-development-report-2021-22), which is measured annually from 1990 to 2022. The mean governance score for each country was calculated based on the Worldwide Governance Indicators (WGI, 2023; https://www.worldbank.org/en/publication/worldwide-governance-indicators), which includes voice and accountability, political stability and absence of violence/terrorism, government effectiveness, regulatory quality, rule of law, and control of corruption⁶¹. Following previous studies^3,62, we first calculated the mean score for each indicator for each country and then calculated the mean governance score for each country between 1996 and 2021. Human population density data were from the GPWv4 database⁶³ (https://doi.org/10.7927/H4F47M2C), which included UN-adjusted estimates of human population density data for the years 2000, 2005, 2010, 2015 and 2020. We obtained data of frequency of armed conflicts from the UCDP Georeferenced Event Dataset (GED) Global version 23.1 (UCDP, 2023; https://ucdp.uu.se/downloads/), which was developed by the UCDP at the PRIO to provide data to a wide range of users on armed conflicts that have occurred since 1946^64,65. To quantify bilateral cooperation abilities between each neighboring country pair, we used a parameter of “Goldstein score”, to determine the extent to which neighbors have been cooperative in the past years⁴. The Goldstein score was quantified as ranging from -10 to 10, with larger values indicating a higher bilateral cooperation ability. Physical fencing distribution is a unique and crucial conservation threat in the transboundary areas⁶. We collected this variable from previous global reports with either fully fenced or under constructions. As physical fencings were mainly located in Eurasia, our supplementary analyses incorporating this variable were only conducted in Eurasia (Supplementary Fig. 15). For bilateral cooperation abilities and physical fencing distribution data, we used the ‘spatial join’ function in ESRI ArcGIS Pro to group the grid cells by neighboring country pair, and assigned the grid cells with the same values of the nearest neighboring country pair. Sizes, years of establishment, and IUCN categories of PAs were extracted from the WDPA database. Previous studies showed that large reserves were more effective than smaller ones in attracting investment in conservation interventions⁶⁶ and protect larger species⁶⁷. We thus included PA size in our model, by assigning grid cells a value of the “GIS_AREA” of their intersected PA. For the protection level of PAs, following the previous approach¹³, we categorized PAs into two general levels: “strict” PAs (IUCN categories I–IV), which are well suited for protecting biodiversity and restricting human activities, and “multiple-use” PAs (IUCN categories V–VI and uncategorized), which permit multiple land uses and resource extraction.

One potential issue with the predictors we used is that there might be differences in the establishment time of PAs and the period of the predictors calculated (i.e., HDI, cropland expansion, human population density and bilateral cooperation ability). Therefore, we extracted the time-related predictor values for each 5 km × 5 km grid for the same period of time when a PA was established. For country-level HDI and governance score, we first identified the grid cells overlapping with the PAs and then extracted and calculated the mean values of HDI and governance score with the same period when the PA was established for each grid cell. For human population density data, we obtained data of the nearest year to the PA establishment for each grid cell.

All continuous predictor variables were normalized (scaled to mean = 0, s.d. = 1) to allow a direct comparison of effect sizes^68,69. One potential issue in our data analysis is that PA coverage might be dependent on the variation in species richness among regions. To evaluate PA coverage by correcting for the effect of species richness, we first fitted a linear model of PA coverage as a function of species richness. Species richness indeed had a positive relationship with PA coverage (Linear Model, estimate ± SEM, 6.7 × 10⁻⁴ ± 2.5 × 10⁻⁶, F_1,363214 = 71,160, P < 0.001). We then extracted the residuals from the model above as a response variable for the further analysis. Species richness was measured as the number of species in each grid cell and was determined by overlaying IUCN range maps of all species recorded with grid cells. Each grid cell was assigned a country identity and neighboring country pair using ‘spatial join’ function in ESRI ArcGIS Pro. Autocorrelation between independent variables would lead to overfitting of models. Therefore, all continuous predictor variables were tested for multicollinearity (for all pairwise correlations, −0.7 < r < 0.7; see Supplementary Fig. 16). We applied linear mixed-effects models as implemented in the R package ‘lme4’⁷⁰ to identify the correlation of ecological and socioeconomic factors with PA coverage for species ranges. Species taxonomic class, country pairing, and country identity were included as random effects to account for the sample pseudoreplications. Specific to the taxonomic class, we calculated the proportion of species ranges protected in each grid cell for different taxa separately, and thus the species taxa information was identified into one of the four with an independent class ID. We then conducted a model averaging procedure to select the best model that explains species geographic ranges protected. To do this, we applied model selection analyses by generating all combinations of the initial variables using the dredge function of ‘MuMIn’ package⁷¹ and ranked them according to the information theory (Akaike’s Information Criterion, AIC; see Supplementary Table 4 for the top models that were within 2 AIC units, i.e., ΔAIC < 2). Our analyses were conducted at both the global scale and continental scale considering the fact that there might be regional heterogeneity in factors influencing the proportion of species ranges covered by PAs. All analyses were conducted in R version 4.2.0 (http://www.r-project.org).

As our analysis is based on grid analysis, we also conducted a sensitivity analysis to account for the potential problem of spatial autocorrelation on the estimation of regression coefficients by applying a generalized additive model using the R package ‘mgcv’⁷² with the same combination of predictor variables as in the LMM plus a tensor product of the longitude and latitude of each grid cell. Furthermore, as relationships between predictor and response variables might be nonlinear, we conducted additional analyses by applying a random forest model using the R package ‘randomForest’⁷³ with the same combination of predictor variables as in the LMM. All analyses were conducted in R version 4.2.0 (http://www.r-project.org).

Global change threats to PAs at global borders

Threats from global changes can cause a large number of species to be threatened¹⁴. Here, we used the spatial overlap between the species ranges covered by PAs and threats from global changes to identify areas with high levels of human pressure but low PA coverage. We used global change maps from Bowler et al. 2020²⁵, which are cumulative indices of spatially explicit datasets of three major anthropogenic stressors: climate change, land use change, and alien species invasion. Climate change, as a recent and accelerating direct driver of biodiversity loss, has attracted attention in transboundary areas^3,74,75. Other drivers, including land use change¹⁶ and alien species invasion²¹, also cause widespread biodiversity loss. These data were consisted of 11 associated variables²⁵ and were contained in numerical raster files at a spatial resolution of 100 km. All the raster files were resampled to a resolution of 5 km × 5 km to match the grid size used in the PA coverage analyses using a bilinear interpolation function that is widely regarded to be more realistic than the simpler nearest-neighbor method⁷⁶. We calculated the mean global change threats for each grid cell by overlaying the grid cells with the anthropogenic stressors raster file. Grid cells were classed as impacted by the threat when they had a mean raster pixel value > 0. To compare the extent of global change threats to PAs between transboundary areas and non-transboundary areas, we applied the resampling approach as mentioned above again⁵⁸. Then, we calculated the proportion of grid cells impacted by the global change threats.

Reporting summary

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