Species across the planet are shifting their ranges to track suitable climate conditions in response to climate change. Given that protected areas have higher quality habitat and often harbor higher levels of biodiversity compared to unprotected lands, it is often assumed that protected areas can serve as steppingstones for species undergoing climate-induced range shifts. However, there are several factors that may impede successful range shifts among protected areas, including the distance that must be travelled, unfavorable human land uses and climate conditions along potential movement routes, and lack of analogous climates. Through a species-agnostic lens, we evaluate these factors across the global terrestrial protected area network as measures of climate connectivity, which is defined as the ability of a landscape to facilitate or impede climate-induced movement. We found that over half of protected land areas and two-thirds of the number of protected units across the globe are at risk of climate connectivity failure, casting doubt on whether many species can successfully undergo climate-induced range shifts among protected areas. Consequently, protected areas are unlikely to serve as steppingstones for a large number of species under a warming climate. As species disappear from protected areas without commensurate immigration of species suited to the emerging climate (due to climate connectivity failure), many protected areas may be left with a depauperate suite of species under climate change. Our findings are highly relevant given recent pledges to conserve 30% of the planet by 2030 (30x30), underscore the need for innovative land management strategies that allow for species range shifts, and suggest that assisted colonization may be necessary to promote species that are adapted to the emerging climate.

Identifying climate analogs

We followed the methods of Abatzoglou et al. (2020) and Parks et al. (2022) to characterize climate and identify backward and forward climate analogs. The specific climate variables we used were average minimum temperature of the coldest month (T_min), average maximum temperature of the warmest month (T_max), annual actual evapotranspiration (AET), and annual climate water deficit (CWD). AET and CWD concurrently account for evaporative demand and availability of water (N. L. Stephenson, 1990). These four variables provide complementary information pertinent to ecological systems and collectively capture the major climatic constraints on species distributions and ecological processes across a range of taxa (Dobrowski et al., 2021; Lutz et al., 2010; Parker & Abatzoglou, 2016; N. Stephenson, 1998; C. M. Williams et al., 2015). Monthly data acquired from TerraClimate (Abatzoglou et al., 2018) were used to produce these annual summaries from 1961-1990 (resolution = ~4km), which were then averaged over the same time period to represent reference period climate normals. The reference time period (1961–1990) is meant to represent climate conditions and climate niches prior to the bulk of recent warming. Future climate conditions were also computed from TerraClimate (available from www.climatologylab.org/terraclimate.html) and correspond to a 2°C increase above pre-industrial levels that are likely to manifest by mid-21^st century without immediate and massive changes in global climate policies (Friedlingstein et al., 2014). As with the reference period climate, we summarized the four +2°C climate metrics annually and over a 30-year time period to represent future climate normals. All analyses in this study were conducted in the R statistical platform (R Core Team, 2020).

We identified backwards and forwards analogs by estimating the climatic dissimilarity between each protected focal pixel (resolution = ~4km to match gridded climate data) and all protected pixels within a 500-km radius using a standardized Mahalanobis distance (Mahony et al., 2017). We chose the 500-km search radius as it encompasses an upper range of dispersal for some terrestrial animals and plants (Chen et al., 2011) when assuming 2°C warming by the mid-21^st century; this search radius has also been used in previous studies (Bellard et al., 2014; Parks et al., 2022; J. W. Williams et al., 2007). The Mahalanobis distance metric synthesized the four climate variables (i.e. T_min, T_max, AET, and CWD; fig. 2a) by measuring distance in multivariate space away from a centroid using principal components analysis of standardized anomalies. Mahalanobis distance scales multivariate mean climate conditions between a pixel and those within the search radius by the focal pixel’s covariance and magnitude of interannual climate variability (ICV) across the four metrics. For backwards analogs, we characterized +2°C ICV and reference period climate normals to calculate climatic dissimilarity; for forward analogs, we used reference period ICV and +2°C climatic normals to calculate climatic dissimilarity. We standardized Mahalanobis distance to account for data dimensionality by calculating a multivariate z-score (σ_d) based on a Chi distribution (Mahony et al., 2017). σ_d represents the climate similarity between each focal pixel and its candidate backward and forward analogs (i.e. all other protected terrestrial pixels within 500 km), and we considered any protected pixels with σ_d ≤ 0.5 as climate analogs (fig. 2b) (following Parks et al., 2022). We were unable to calculate Mahalanobis distance when there was no ICV for any one of the four variables, and as a consequence, these areas are omitted from all analyses; this affects, for example, a relatively small tropical area in South America (CWD=0 each year) and areas perennially covered by snow (CWD=0 each year; e.g. most of Greenland).

We focused our analyses on protected areas as defined by the World Database on Protected Areas (WDPA) (IUCN & UNEP-WCMC, 2019) and included protected areas classified as IUCN (International Union of Conservation for Nature) Management Categories I-VI, except those identified as ‘proposed’, ‘marine’, or otherwise aquatic (e.g. wetland, riverine, endorheic). A large number of protected areas, however, were not assigned an IUCN category in the WDPA (identified as ‘Not Reported’, ‘Not Assigned’, or ‘Not Applicable’) but are likely to have reasonably high levels of protection (e.g. Kruger National Park in South Africa). We included these additional protected areas if the level of human modification was similar or less than that observed within IUCN category I-VI protected areas. To do so, we measured mean land-use intensity within each IUCN category I-VI protected area using the Human Modification Gradient (HMG) raster dataset (Kennedy et al., 2019) and calculated the 80^th percentile of the resulting distribution. Any unassigned protected areas with a mean HMG less than or equal to this identified threshold were included in our study (following Dobrowski et al., 2021). We then converted this vector-based polygon dataset to raster format (resolution = ~4km to match gridded climate data; n=1,063,748 pixels). It is well-recognized that the WDPA contains a large number of duplicate and overlapping polygons (Palfrey et al., 2022; Vimal et al., 2021). Although this does not affect summaries across the globe or for individual countries (described below), it provides a challenge when trying to summarize by individual protected areas (due to double-counting). Consequently, we ‘cleaned’ the WDPA prior to summarizing the climate connectivity metrics for individual protected areas by removing polygons that exhibited ≥ 90% overlap with another; this resulted in 29,752 individual protected areas (available in the Electronic Supplemental Material).

Least-cost path modelling

Following Dobrowski and Parks (2016) and Carroll et al. (2018), we used least-cost path modelling (Adriaensen et al. 2003) to build potential climate-induced movement routes between each protected focal pixel and its backward and forward analogs. The least-cost models were parameterized with resistance surfaces based on climate dissimilarity and the human modification gradient (HMG) (Kennedy et al., 2019). For backward analog modelling, we characterized climatic dissimilarity (i.e. climatic resistance) using two intermediate surfaces, the first being the Mahalanobis distance between each focal pixel (using +2°C ICV) and all other pixels using reference period climate normals (fig. 2c) and the second being the Mahalanobis distance (using +2°C ICV) and all other pixels using +2°C climate normals (fig. 2d). These two surfaces provide a proxy for climate similarity designed to capture transient changes between the reference period and +2°C climate; these were then averaged to characterize the overall climatic resistance across time and space (fig. 2d). For forward analog modelling, the process is similar except we used reference period ICV when characterizing climatic resistance (fig. 2a-2d). We then multiplied the climatic resistance (fig. 2d) by HMG (fig. 2e) to create the final resistance surface for least-cost path modeling (cf. Parks et al., 2020). Prior to this step, we rescaled HMG from its native range (0–1) to 1–25 to correspond with the range of Mahalanobis distance values and thereby grant comparable weights to climatic resistance and HMG resistance (~95% of all Mahalanobis distance values are below 25 within a 500km radius). Open water was given a resistance=25 so that paths would avoid water when possible. Least-cost path modelling was achieved using the gdistance package (van Etten, 2017); paths represent the least accumulated cost across the final resistance surface (fig. 2f) between each focal pixel and analog (fig. 2g). Because paths were rarely straight lines, some were longer than the 500km that we established as a search radius. We removed these longer paths to abide by the biologically informed upper dispersal constraint.

Calculating climate connectivity metrics and climate connectivity failure

We calculated the length (i.e. dispersal exposure), land-use modification (i.e. human exposure), and climatic resistance (i.e. climate exposure) for each path, remembering that each focal pixel may have many analogs and resultant paths. Human exposure represents cumulative HMG (fig. 2e) across all pixels in a path and climate exposure represents cumulative climate resistance (fig. 2d) along a path. Human exposure and climate exposure were calculated by multiplying the mean HMG (unscaled; fig. 2f) and mean climate resistance (fig. 2d) along each path by the length of each path, respectively. Each path’s climate connectivity metric (dispersal, human, and climate exposure) was converted to a percentile (range = 0–100) to facilitate easier interpretation and comparison among metrics; relative to other protected pixels, small percentiles represent low exposure and large percentiles represent elevated exposure. We summarized (i.e. averaged the percentiles) dispersal exposure, human exposure, and climate exposure across each protected focal pixel (again, remembering that each pixel may have multiple analogs and resultant paths). Our fourth climate connectivity metric, analog exposure, can’t be summarized on a per-path basis, because by definition, there is no least-cost path when there are no protected climate analogs. Instead, protected pixels either do or do not have protected climate analogs.

Focal pixels were identified as exhibiting climate connectivity failure when they exceeded the 75^th percentile for dispersal or climate exposure, exceeded the 90^th percentile for human exposure, or had no protected climate analog. We assumed that focal pixels exceeding these percentiles are located in landscapes that hinder successful range shifts among protected areas (i.e. climate connectivity failure) for a non-negligible proportion of extant species, considering that the biodiversity at a given site comprises mammals, birds, insects, mollusks, amphibians, reptiles, fish, crustaceans, annelids, vascular plants (e.g. trees grasses, shrubs), and non-vascular plants (e.g. fungi, mosses, lichens). The numerous and diverse species at a given site have a wide range of dispersal abilities, sensitivities to human land uses, and climatic tolerances. We used a higher threshold (90^th percentile) for describing climate connectivity failure due to human exposure because large, remote protected areas in the network skew human exposure towards lower values from a global perspective. These percentile thresholds are likely conservative when considering the large number and diversity of species at a given site. In terms of dispersal, for example, many species have maximum dispersal capabilities on the range of 1 km/year or less (Jenkins et al., 2007; McLachlan et al., 2005; Schwartz et al., 2001). This represents dispersal of 75 km under 2°C warming in the 75 years covering the midpoint of the reference period (1975) to mid-21^st century. In our study, the 75^th percentile path length, corresponding to dispersal exposure, is ~385 km, well above such dispersal limits, supporting our assertion that the 75^th percentile is conservative for estimating climate connectivity failure. Furthermore, the mean HMG value for a 100km path at the 90^th percentile threshold is 0.22, which is well above the 0.1 threshold that Brennen et al. (2022) used to identify areas moderately to highly impacted by human land-uses. Lastly, the mean climatic distance for a 100km path at the 75^th percentile is well over two standard deviations different, on average, from the focal pixel and analog.

We report the percent of protected pixels across the globe and within each country that exhibits climate connectivity failure. We also assessed the potential for each of the 29,752 individual protected areas (e.g. Yellowstone National Park, Serengeti National Park) to undergo climate connectivity failure using a slightly different method. To do so, we calculated the mean percentile among pixels within each protected area for each of dispersal exposure, human exposure, and climate exposure (each metric was averaged across a protected area; the metrics themselves were not averaged with each other). We then calculated the percent of each protected area that did not have a protected climate analog (analog exposure). Although a binary approach (has or does not have an analog) is appropriate when evaluating individual focal pixels, a percent-based valuation is most appropriate and informative when evaluating individual protected areas with up to thousands of pixels. Individual protected areas exhibited climate connectivity failure if the mean dispersal exposure or climate exposure exceeded the 75th percentile, mean human exposure exceeded the 90^th percentile, or the analog exposure exceeded 75%.

References

• Abatzoglou, J. T., Dobrowski, S. Z., & Parks, S. A. (2020). Multivariate climate departures have outpaced univariate changes across global lands. Scientific Reports, 10(1), Article 1. https://doi.org/10.1038/s41598-020-60270-5
• Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A., & Hegewisch, K. C. (2018). TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Scientific Data, 5(1), Article 1. https://doi.org/10.1038/sdata.2017.191
• Bellard, C., Leclerc, C., Leroy, B., Bakkenes, M., Veloz, S., Thuiller, W., & Courchamp, F. (2014). Vulnerability of biodiversity hotspots to global change. Global Ecology and Biogeography, 23(12), 1376–1386. https://doi.org/10.1111/geb.12228
• Brennan, A., Naidoo, R., Greenstreet, L., Mehrabi, Z., Ramankutty, N., & Kremen, C. (2022). Functional connectivity of the world’s protected areas. Science, 376(6597), 1101–1104. https://doi.org/10.1126/science.abl8974
• Carroll, C., Parks, S. A., Dobrowski, S. Z., & Roberts, D. R. (2018). Climatic, topographic, and anthropogenic factors determine connectivity between current and future climate analogs in North America. Global Change Biology, 24(11), 5318–5331. https://doi.org/10.1111/gcb.14373
• Chen, I.-C., Hill, J. K., Ohlemüller, R., Roy, D. B., & Thomas, C. D. (2011). Rapid Range Shifts of Species Associated with High Levels of Climate Warming. Science, 333(6045), 1024–1026. https://doi.org/10.1126/science.1206432
• Dobrowski, S. Z., Littlefield, C. E., Lyons, D. S., Hollenberg, C., Carroll, C., Parks, S. A., Abatzoglou, J. T., Hegewisch, K., & Gage, J. (2021). Protected-area targets could be undermined by climate change-driven shifts in ecoregions and biomes. Communications Earth & Environment, 2(1), Article 1. https://doi.org/10.1038/s43247-021-00270-z
• Dobrowski, S. Z., & Parks, S. A. (2016). Climate change velocity underestimates climate change exposure in mountainous regions. Nature Communications, 7(1), Article 1. https://doi.org/10.1038/ncomms12349
• Friedlingstein, P., Andrew, R. M., Rogelj, J., Peters, G. P., Canadell, J. G., Knutti, R., Luderer, G., Raupach, M. R., Schaeffer, M., van Vuuren, D. P., & Le Quéré, C. (2014). Persistent growth of CO2 emissions and implications for reaching climate targets. Nature Geoscience, 7(10), Article 10. https://doi.org/10.1038/ngeo2248
• IUCN & UNEP-WCMC. (2019). Protected Planet: World Database on Protected Areas (WDPA). Accessed September 2019. Available at www.protectedplanet.net. (Accessed September 2019) [Map]. www.protected.planet.net
• Jenkins, D. G., Brescacin, C. R., Duxbury, C. V., Elliott, J. A., Evans, J. A., Grablow, K. R., Hillegass, M., Lyon, B. N., Metzger, G. A., Olandese, M. L., Pepe, D., Silvers, G. A., Suresch, H. N., Thompson, T. N., Trexler, C. M., Williams, G. E., Williams, N. C., & Williams, S. E. (2007). Does size matter for dispersal distance? Global Ecology and Biogeography, 16(4), 415–425. https://doi.org/10.1111/j.1466-8238.2007.00312.x
• Kennedy, C. M., Oakleaf, J. R., Theobald, D. M., Baruch-Mordo, S., & Kiesecker, J. (2019). Managing the middle: A shift in conservation priorities based on the global human modification gradient. Global Change Biology, 25(3), 811–826. https://doi.org/10.1111/gcb.14549
• Lutz, J. A., van Wagtendonk, J. W., & Franklin, J. F. (2010). Climatic water deficit, tree species ranges, and climate change in Yosemite National Park. Journal of Biogeography, 37(5), 936–950. https://doi.org/10.1111/j.1365-2699.2009.02268.x
• Mahony, C. R., Cannon, A. J., Wang, T., & Aitken, S. N. (2017). A closer look at novel climates: New methods and insights at continental to landscape scales. Global Change Biology, 23(9), 3934–3955. https://doi.org/10.1111/gcb.13645
• McLachlan, J. S., Clark, J. S., & Manos, P. S. (2005). Molecular indicators of tree migration capacity under rapid climate change. Ecology, 86(8), 2088–2098. https://doi.org/10.1890/04-1036
• Palfrey, R., Oldekop, J. A., & Holmes, G. (2022). Privately protected areas increase global protected area coverage and connectivity. Nature Ecology & Evolution, 6(6), Article 6. https://doi.org/10.1038/s41559-022-01715-0
• Parker, L. E., & Abatzoglou, J. T. (2016). Projected changes in cold hardiness zones and suitable overwinter ranges of perennial crops over the United States. Environmental Research Letters, 11(3), 034001. https://doi.org/10.1088/1748-9326/11/3/034001
• Parks, S. A., Carroll, C., Dobrowski, S. Z., & Allred, B. W. (2020). Human land uses reduce climate connectivity across North America. Global Change Biology, 26(5), 2944–2955. https://doi.org/10.1111/gcb.15009
• Parks, S. A., Holsinger, L. M., Littlefield, C. E., Dobrowski, S. Z., Zeller, K. A., Abatzoglou, J. T., Besancon, C., Nordgren, B. L., & Lawler, J. J. (2022). Efficacy of the global protected area network is threatened by disappearing climates and potential transboundary range shifts. Environmental Research Letters, 17(5), 054016. https://doi.org/10.1088/1748-9326/ac6436
• R Core Team. (2020). R: A language and environment for statistical computing.
• Schwartz, M. W., Iverson, L. R., & Prasad, A. M. (2001). Predicting the potential future distribution of four tree species in Ohio using current habitat availability and climatic forcing. Ecosystems, 4(6), 568–581. https://doi.org/10.1007/s10021-001-0030-3
• Stephenson, N. (1998). Actual evapotranspiration and deficit: Biologically meaningful correlates of vegetation distribution across spatial scales. Journal of Biogeography, 25(5), 855–870. https://doi.org/10.1046/j.1365-2699.1998.00233.x
• Stephenson, N. L. (1990). Climatic Control of Vegetation Distribution: The Role of the Water Balance. The American Naturalist, 135(5), 649–670. https://doi.org/10.1086/285067
• van Etten, J. (2017). R Package gdistance: Distances and Routes on Geographical Grids. Journal of Statistical Software, 76, 1–21. https://doi.org/10.18637/jss.v076.i13
• Vimal, R., Navarro, L. M., Jones, Y., Wolf, F., Le Moguédec, G., & Réjou-Méchain, M. (2021). The global distribution of protected areas management strategies and their complementarity for biodiversity conservation. Biological Conservation, 256, 109014. https://doi.org/10.1016/j.biocon.2021.109014
• Williams, C. M., Henry, H. A. L., & Sinclair, B. J. (2015). Cold truths: How winter drives responses of terrestrial organisms to climate change. Biological Reviews, 90(1), 214–235. https://doi.org/10.1111/brv.12105
• Williams, J. W., Jackson, S. T., & Kutzbach, J. E. (2007). Projected distributions of novel and disappearing climates by 2100 AD. Proceedings of the National Academy of Sciences, 104(14), 5738–5742. https://doi.org/10.1073/pnas.0606292104

There are three files in this repository:

1)      backward.analogs - master.table.xlsx – results for backward analogs:

·         Each climate connectivity metric (dispersal exposure, human exposure, climate exposure, and analog exposure) is summarized by country; percent protected lands in each country that exhibit climate connectivity failure is also indicated.  

·         Each climate connectivity metric (dispersal exposure, human exposure, climate exposure, and analog exposure) is summarized by protected area. Values represent the mean pixel-based percentile. Also included is a binary (0, 1) indicator of whether the protected area exhibits climate connectivity failure.

2)      forward.analogs - master.table.xlsx – results for forward analogs:

·         Each climate connectivity metric (dispersal exposure, human exposure, climate exposure, and analog exposure) is summarized by country; percent protected lands in each country that exhibit climate connectivity failure is also indicated.  

·         Each climate connectivity metric (dispersal exposure, human exposure, climate exposure, and analog exposure) is summarized by protected area. Values represent the mean pixel-based percentile. Also included is a binary (0, 1) indicator of whether the protected area exhibits climate connectivity failure.

3)      PA_shapefile - cleaned.zip: This is the ‘cleaned’ (see Methods) protected area shapefile we used as a way to summarize dispersal exposure, human exposure, climate exposure, and analog exposure for each protected area. 

Note that two of these files are Microsoft Excel; they should be accessible via LibreOffice and R and potentially other open-source alternatives.