Protected-area targets could be undermined by climate change-driven shifts in ecoregions and biomes

Expanding the global protected area network is critical for addressing biodiversity declines and the climate crisis. However, how climate change will affect ecosystem representation within the protected area network remains unclear. Here we use spatial climate analogs to examine potential climate-driven shifts in terrestrial ecoregions and biomes under a +2 °C warming scenario and associated implications for achieving 30% area-based protection targets. We find that roughly half of land area will experience climate conditions that correspond with different ecoregions and nearly a quarter will experience climates from a different biome. Of the area projected to remain climatically stable, 46% is currently intact (low human modification). The area required to achieve protection targets in 87% of ecoregions exceeds the area that is intact, not protected, and projected to remain climatically stable within those ecoregions. Therefore, we propose that prioritization schemes will need to explicitly consider climate-driven changes in patterns of biodiversity. Preservation of biodiversity by expanding protected areas is complicated by expected climate change that could dissociate ecoregion or biome from climate conditions in over half of land area, according to analyses of spatial climate analogs

C limate change and biodiversity loss are among the greatest challenges our society faces today and in the coming century 1 . Habitat loss and degradation and their interaction with climate change imperil the planet's biodiversity and the systems that support civilization 2,3 . Protected areas (PAs) are the foundation of modern-day conservation approaches. Their importance in protecting biodiversity has been demonstrated across the globe 4,5 . A clear scientific consensus has emerged that expanding the PA network is critical for maintaining and restoring intact natural ecosystems 6 , for protecting biodiversity, supporting ecosystem services, and for achieving scalable natural climate solutions 7,8 .
As the intactness of the worlds ecosystems declines, there are growing calls to focus PA network expansion on ecosystem-level targets 9,10 . For example, Dinerstein et al. 11 and Chauvenet et al. 12 proposed that global targets for expansion of the protected area system be based on ecoregions as surrogates for ecosystem-level biodiversity. Dinerstein et al. 7 assert that "ecoregions effectively represent similar clusters of not only habitat types but also species", and thus their use ensures the "creation of networks of protected areas that represent the widest array of habitats and, by extension, conserve the widest range of species and their unique adaptations to their environments." These authors thereby call for a "Global Deal for Nature" with area-based targets of 30% for each ecoregion and an additional 20% of the planet designated as climate stabilization areas: in sum, half the planet protected 13 . In such proposals, ecoregions serve not only as ecologically based planning units, but also as a basis for assessing representation of conservation targets within global prioritization processes. The use of ecosystem area-based targets has been questioned 14,15 . Nonetheless, ecosystem-level targets are actively being proposed for the 2030 follow on to the 2010-2020 Convention on Biological Diversity (CBD) 6 .
Although ecosystem-level PA prioritization schemes strive to capture a breadth of natural processes, habitats, and species, climate change may undermine the efficacy of the PA network-even an expanded one 16 . If the past is prologue, species assemblages and ecosystems will change over time, driven by species range shifts, shifts in abundance, and extirpation or extinction [17][18][19] . A growing body of research suggests that terrestrial lands within the protected area estate 16,20 and globally are experiencing heightened climate change exposure 21 and relatively novel climatic conditions 22 . These changes and those projected for the future are expected to result in changes in the extent of biomes [23][24][25] and habitat transformation in terrestrial ecoregions globally 26 .
Collectively, the impermanence of species assemblages, communities, and ecosystems pose a challenge to conservation frameworks that rely on protected areas with static boundaries. Conservation plans based on current geographic patterns of biodiversity may be insufficient to support future biota and natural processes 27 and may fail to afford species access to suitable climates as the Earth warms 28,29 . These challenges raise questions about the efficacy of the existing PA network and how to expand its coverage under a warming climate. With renegotiations of global conservation targets on the horizon and persistent failures to meet existing targets (e.g., 2020 Aichi Biodiversity Targets 30 ), the global conservation community must devise a strategic and flexible PA prioritization approach that expressly addresses climate change and anticipated shifts in biodiversity.
Prioritization schemes for expanding the network tend to fall on a spectrum with respect to how or if they explicitly consider the challenges posed by climate change. On one end are static place-based approaches that focus on existing patterns of biodiversity, are characterized by discrete PA boundaries, and do not explicitly consider climate change. Examples include key biodiversity areas (KBAs) 15 and threatened ecosystems (e.g., IUCN Red List of Ecosystems) 9 . Regardless of the particular focus, these types of place-based strategies treat planning units and the features they represent as static over conservation planning horizons and thus are agnostic to climatic changes per se. At the other end of the spectrum are climate space approaches that focus on representation of climate types under current and future conditions [31][32][33] . These approaches explicitly address climate changes under the premise that protecting a diversity of climate types will in turn protect biodiversity. This assumption deserves scrutiny: biodiversity is not evenly distributed in climate space and though climate is an important filter, other factors affect the distribution of biota including biogeographic and evolutionary history, disturbance, geologic and edaphic factors, dispersal constraints, and biotic interactions [34][35][36] . Moreover, there are multiple approaches that fall between these end-members and borrow from both dynamic and place-based strategies 37,38 .
Here we examine how potential climate-driven changes in the distribution of terrestrial ecoregions and biomes may alter representation within the global PA network, and we consider the implications for area-based conservation targets that call for 30% of Earth to be formally protected 7 . To assess this, we use spatial climate analogs-an approach that borrows from strengths of both place-based and climate-space prioritization schemes. Spatial climate analogs use space-for-time substitution to identify present day locations that share similar climates to those projected for a focal location in the future (Fig. 1). We examine climate projected for mid-century with global mean temperatures that are +2°C above pre-industrial levels. We identify analogs within a 2000 km search radius for all global terrestrial surfaces using four climate variables that represent temperature and climatic water balance, and are strong predictors of species distributions across taxa 39,40 . We quantify climatic similarity between focal sites and candidate analogs based on a multivariate measure of climate distance that is scaled by the interannual climate variability at a focal site.
We project changes to the distribution and extent of biomes and the terrestrial ecoregions they encompass, including the emergence of regionally novel climate conditions that may drive the development of novel ecosystems 41 . We then assess the proportion of each ecoregion that is projected to: (1) remain climatically stable with conditions similar to the current ecoregion (and by default, biome), (2) transition to conditions representative of other ecoregions or biomes, and (3) remain stable and intact (low human modification), as these may represent strategic opportunities for expanding the PA estate. We make a database of climate analogs publicly available for use in conservation and climate adaptation planning and for communicating climate change impacts to the broader public. To facilitate the latter, we provide a web-based tool (https://plus2c.org) that allows users to visualize potential future impacts of climate change (including ecoregions shifts) at user-selected sites.

Results
We found that by mid-century, 53.9% of global land area will experience climatic conditions that correspond to different ecoregions (Fig. 2). A larger proportion (58.2%) of land area within PAs is projected to transition between ecoregions (p = 0.01 based on Wilcoxon signed rank test). At the broader characterization of biomes, 21.8% of land area is projected to transition between biomes, which is consistent with biome transitions projected within PAs (20.6%) (p = 0.58). The prevalence of transitions varies by latitude, with peaks occurring in northern mid-latitude regions (+30°C to +60°C), and southern tropical to sub-tropical regions (0°C to −30°C; Fig. 2). These latitudinal zones roughly correspond with low protected area coverage (Fig. 2). Regionally, we found areas with low protection and high transition potential (for both ecoregions and biomes) in the southeastern U.S., northern India, central Asia, central Africa, and Boreal forests (Fig. 2). On the other hand, we found areas of high protection and high transition potential for ecoregions in the Amazon Basin and central Australia (Fig. 2a). Lastly, 3.6% of total land area and 4.6% of area within PAs are projected to have regionally non-analogous climates (i.e., no current climates within 2000 km are analogous to those projected for the focal site in the future) including extensive areas within the Amazon, Southeastern U.S., Tropical and Subtropical regions of Africa, the Arabian Peninsula, India, and Northern Australia (Fig. 2). An interactive version of Fig. 2 can be found at https://umlandscapee cology.shinyapps.io/PA_shiny_map/.
We examined states and transitions that occurred under climate change between biomes across all land outside of and within PAs (Fig. 3   Modeling climate change impacts on conservation planning units using spatial climate analogs. Under current conditions (a), reverse analogs (sites that have climates that are similar to those projected in the future for a focal cell) were identified within a fixed search radius (dotted line) (adapted from ref. 83 ). Reverse analogs were then used in conjunction with a map of terrestrial ecoregions (b) to project potential changes in ecoregions and biomes. In this study we assess compositional changes at focal sites based on a plurality vote from 100 analog sites (c). Protection status. When averaged across all ecoregions, the change in the percent of protection (ΔPP) between the reference period and mid-century is minimal (current = 15.4%; +2°C = 16.1%). This consistency belies larger transitions that occur at the biome and ecoregion scale (Fig. 4). Based on a bootstrapping procedure (see methods), we anticipate significant increases in protection for Tropical Dry Broadleaf Forests, Temperate Grasslands and Shrublands, and Mediterranean Forests and Scrub (Fig. 4). We anticipate declines in protection for Tundra and Montane Grassland and Shrublands (n.b. a bootstrapped 90% confidence interval of ΔPP for these biomes overlaps zero in both cases).
Nearly seventeen percent (16.8%; n = 128) of ecoregions currently meet or exceed a 30% protection target. Under a +2°C scenario, these 128 ecoregions lose 16% of their protected area on average and roughly a ¼ of these will fall below a 30% protection target. Additionally, for individual ecoregions, large transitions in protection status are apparent: 8.6% of ecoregions are projected to have gains or losses in ΔPP that exceed 20%. Additionally, 56 ecoregions 'disappeared' under the +2°C scenario (were not found to be a reverse analog for any location), with nearly half of these in the Tropical & Subtropical Moist Broadleaf Forests Biome (n = 24). This biome contains the most ecoregions (n = 216) and has a large area of novel climate conditions.
Area-based targets and stable and intact areas. The availability of stable, unprotected, and intact area for expanding the PA estate varies regionally and by biome. Averaged across ecoregions, 40.6% of current ecoregion area is projected to remain within the same ecoregion (stable) under a + 2°C scenario. However, this area is inadequate to meet area-based targets for a majority of ecoregions. For example, 55.9% of ecoregions have protected area targets that exceed the area projected to remain stable for that ecoregion. If we consider only areas that are currently unprotected, intact (low human modification), and are expected to remain climatically stable, this value increases to 86.8%. Unprotected, intact, and stable areas principally lie in high latitude (Boreal Forests, Tundra), mountainous (Temperate Conifer Forests), and desert regions (Deserts and Xeric Shrub) where human modification is limited (Fig. 5). By contrast, many ecoregions within Tropical Biomes, Temperate Grasslands, and Temperate Broadleaf Forests have limited opportunities for expanding the PA estate either due to limited stable area, high human modification, or both (Fig. 5).
Uncertainty in Projections. We use two approaches to characterizing the uncertainty of our projections. The first relies on the dissimilarity in climate (σ d ) between projected future conditions for focal locations and their candidate analogs (see methods). For instance, sites with low σ d values have analog locations that are climatically similar to conditions projected for a focal location whereas high σ d values suggest dissimilar conditions. We found a latitudinal gradient in σ d was apparent. Climate dissimilarity (σ d ) values were higher in tropical and sub-tropical ecoregions associated with areas of low interannual climate variability and warm climates and they declined toward the poles (Fig. S1). The second approach to characterizing uncertainty relies on the plurality of votes used to assign ecoregions to a given focal cell (see methods). We assume greater confidence in a projected ecoregion shift if the plurality is large (i.e., if a large proportion of the 100 reverse analogs agree on the same projected change). Consistent with patterns in σ d , we found that the proportion of votes that comprised the plurality was smallest in For chord diagrams, the relative magnitude of biome-level fluxes is shown by the width of arrows between biomes. The net change in biome area (percent) is shown in parenthesis. The proportion of the circle's circumference associated with a given biome indicates the total amount of land area that is or will be in that biome. Current biome land area is the circumference proportion from which arrows of that biome's color originate; future biome land area is the circumference proportion of that biome's color to which arrows of any color point. Arrows that start and end at the same color represent areas that do not undergo a biome-level transition.
tropical latitudes (suggesting lower confidence) and showed regional variation with lower values in regions with complex topography (Fig. S2). Nonetheless, there was strong agreement amongst the prospective analogs; for most areas the percent assigned to the plurality ranged between 60 and 90%.

Discussion
Climate change has the potential to substantially alter the distribution of the planet's terrestrial ecoregions and biomes. Our results suggest that by mid-century, 54% and 22% of global land area will experience climatic conditions that correspond to different ecoregions and biomes, respectively. Further, 3.6% of land area will experience climatic conditions that are regionally novel, which is consistent with previous research that shows recent climate novelty is most prevalent in the tropics 22,42 (Fig. 2). Our results are also broadly consistent with previous studies which examine climate change exposure of protected areas globally. For example, Hoffmann and Beierkuhnlein 20 showed that climate change exposure in PAs was projected to be most pronounced in tropical and subtropical regions. Our findings suggest these areas, along with northern mid-latitude regions, which have the least PA coverage and greatest human footprint, are poised to transition to other ecoregions and biomes (Fig. 2).
If realized, climate-driven transitions will affect the degree to which individual ecoregions are protected and thus representation within the PA network as a whole. For example, more than ten percent of ecoregions within the Mediterranean and Tundra biomes are projected to experience gains or losses in protection status that exceed 20%. Similarly, for 128 ecoregions (16.9% globally) that currently have at least 30% area protected, we project an average loss of 16% in protected area. Roughly a quarter of these ecoregions would thereby fall below the 30% target under the +2°C scenario. More broadly, we anticipate losses in protection within the Tundra biome and gains in protection in the Tropical Dry Broadleaf Forest, Mediterranean Forest and Scrub, and Temperate Grass and Shrubland biomes (Fig. 4). Moreover, our projections suggest that 56 ecoregions (~6%) will disappear under the +2°C scenario, which if manifest, would represent a conservation loss that is not accounted for in current targets.
Shifting representation within the PA network will challenge our ability to define conservation targets (e.g., area-based protection targets) and develop accounting strategies and indicators for monitoring progress towards those targets at the global and national scale. The 2011-2020 Aichi Biodiversity Target 11 called for protection of at least 17% of terrestrial areas by 2020. The post-2020 CBD framework proposes expanding this target to 30%. A key indicator under Target 11 is PA coverage within ecoregions because this metric is considered "Useful for assessing PA network representativeness, and can be applied at a range of scales (national and global)" 43 . The CBD's forthcoming meeting, termed the Global Biodiversity Framework (GBF), will include Opportunity for expanding the protected area estate. For each ecoregion, the area needed to achieve a 30% area-based protection target (after accounting for existing PAs) was compared against the remaining area within the ecoregion that is currently unprotected, intact (low human modification), and projected to remain climatically stable. A value of 100% means that the remaining unprotected, intact, and stable area is equal to the area needed to achieve the 30% target. Values below 100% imply that the target area exceeds the area that meets these criteria whereas values above imply more area is available than the target. An ecoregion level map of 'opportunity' for expanding the PA estate is presented in a and summarized by biome in b. Ecoregions in brown currently meet the 30% protection target. Areas in black represent no-analog conditions.
Boreal Forests (26) Deserts & Xeric Shrub. (100) (40) Mont. Grass & Shrublands (46) Temp. Broadlf Forests (82) Temp. Conifer Forests (47) Temp. Grass & Shrublands (44) Trop. Conifer Forests (14) Trop. Dry Broadlf Forests (56) Trop. Grass & Shrublands (54) Trop. Moist Broadlf Forests (216) Tundra ( targets and indicators designed to reflect the status of global biodiversity. For example, proposed targets related to biodiversity status and trend include limits on the degradation and loss of intact ecosystems 6 . Our results provide a set of climate-informed null expectations that scientists and policy makers can consider when evaluating indicators for the post 2020 GBF for assessing future losses and gains in the extent of ecosystems. The potential for novel geographic patterns of biodiversity will also challenge how we strategically expand the protected area network 31,44 . The need to anticipate and plan for changes in biodiversity patterns is becoming more acute as we increasingly fall behind global emissions targets set by the Paris Accord. Assuming an RCP 8.5 emissions trajectory, climate models, including those used here, show the globe approaching +2°C around 2050. Dinerstein et al. 7 and others highlight the value of formally protecting 30% of each ecoregion and an additional 20% of terrestrial land area for promoting natural climate solutions to stave off 1.5°C of warming, a value that is considered by some to be a climatic "tipping point" 45 . Unfortunately, current emission trajectories almost certainly lock us into exceeding this threshold and put us on a path of broad-scale disruptions to the planet's ecosystems 3 .

Mediterr. Forests & Scrub
Conservation initiatives should consider emphasizing prioritization schemes that enhance the resilience of the PA network to the emergence of novel geographic patterns in the planet's biodiversity. A logical starting point for this is to identify areas that currently fall short of area-based targets, have minimal human modification, and are projected to remain climatically stable over specific planning horizons. At the global scale, human modification acts as a dominant constraint 46 with opportunities for PA expansion found principally in high latitude ecoregions, mountainous regions, and in Deserts and Scrublands globally (Fig. 5). Moreover, we find that prioritizing PA network expansion based solely on climatically stable areas will not be possible: 56% of ecoregions have current protected area targets that exceed the area projected to remain stable in those ecoregions (Fig. 5). If we integrate all three filters and consider areas that are not currently protected, intact, and projected to remain stable, we find that 87% of ecoregions do not have enough area that meet these criteria to achieve current protected area targets. Thus focusing conservation efforts towards stable and intact areas is necessary but insufficient to meet current conservation targets.
Expanding the PA estate will thus require investing in PAs that will experience changing patterns of biodiversity and even regionally novel conditions. At the global scale, where should these investments be made? Outside of high latitude regions, our findings identify opportunities in eastern South America, central and southern Africa, central Asia, northern Mexico, and southern Australia (Fig. 5). These areas have remnant intact landscapes and are projected to have moderate transition potential. These regions are also hotspots for biodiversity and endemism 47 which increases the imperative for developing climate informed strategies. Our results also suggest that beyond stable and intact areas we must consider investing in matrix habitats and human modified landscapes for maintaining biological diversity [48][49][50] .
In addition to a global-scale assessment, climate analogs can be used at regional scales to identify stable areas which can serve as conservation anchor points-large stable PAs that are projected to have minimal transitions-and transitional zones between current conservation targets and their suitable climates in the future. Climate analogs are likely to be sources of biota that are climatically preadapted for focal locations 51,52 . This presents an opportunity to focus conservation efforts on ensuring climatic connectivity between existing and candidate PAs and the sources of biodiversity that may colonize them 19 . Recent research has combined climate analog and connectivity modeling approaches to develop novel metrics of climate change exposure. These include distance to analogs (sensu climate change velocity 53 ) and exposure to dissimilar climates and anthropogenic stressors along potential migration pathways [53][54][55] . Combined with observed biodiversity patterns, these exposure metrics could be used to strategically expand the PA network in order to promote connectivity between refugia via the delineation of "climate corridors" 55 . Moreover, dynamic spatial conservation planning methods have been developed using complementary fine-filter approaches that account for shifting species distributions under climate change 37,56 . These approaches identify areas that are expected to remain climatically suitable (stable) for individual species over a given planning horizon and areas that provide dispersal pathways between stable habitat and locations that will be climatically suitable in the future 38 .
To provide a concrete regional example of how a spatial prioritization scheme can leverage analog-based projections, we provide a brief case-study for the Yellowstone to Yukon (Y2Y) region in western North America (Fig. 6) using analogs from the database we make publicly available. Large portions of Yellowstone, Glacier, Waterton, Banff and Jasper National Parks, as well as surrounding unprotected areas, are projected to remain stable and act as persistent refugia over the timeframe of our analysis. These PAs, oriented in a south to north sequence along the Rocky Mountain Cordillera, are interspersed with other areas-both protected and not protected-that are poised to transition to more warm-adapted ecosystems. Our projections can inform decisions on how to regionally expand PAs to promote climate connectivity via latitudinal and elevational gradients 44 while indicating whether additions would act as transient conservation stepping stones 57 or more permanent refugia. Additionally, the location of reverse analogs (some of which occur in other PAs) represents potential sources of biota that may colonize these PAs in the future (Fig. 6). This can provide information to PA managers for anticipating shifts in biotic assemblages as the climate warms. Importantly, these analog-based forecasts of biodiversity Static (not protected) Ecoregion transition (not protected) Biome transition (not protected) Static (protected) Ecoregion transition (protected) Biome transition (protected) Insufficient data 800 km Fig. 6 Regional case study of the use of spatial climate analogs for spatial prioritization along the Yellowstone to Yukon corridor in western North America. Jasper National Park (NP), Banff NP, Waterton NP, Glacier NP, The Bob Marshall Wilderness, and Yellowstone NP are outlined in red from north to south. For a random subset of 10 sites within each of these protected areas that are projected to transition, we show the location of 10 reverse analogs (black points). change are inherently adaptive and could be developed over multiple time steps and could iteratively incorporate new climatic and biodiversity data to support dynamic spatial prioritization routines (sensu 37 ).
Increasingly there are calls for integrated and scalable approaches for monitoring biodiversity trends and identifying risks to socio-ecological systems from global change. A spatial climate analog approach provides a means to anticipate climatedriven changes to patterns of biodiversity for spatially discrete planning units, which feature prominently in proposed PAexpansion schemes. By projecting changes to biogeographically informed planning units such as ecoregions, we generate information about the relative climate change exposure of such units. By delineating climatically stable areas within each ecoregion, we retain the value of the ecoregional classification system while acknowledging that ecological communities will change as species shift. Additionally, we posit that analog modeling approaches have key traits that make them amenable for conservation applications. First, the approach is generalizable: climate analogs can be used across multiple sectors for adaptation planning, for modeling changes in ecosystem state variables, and for conservation planning. Second, the approach is scalable: comparable methods can be applied at both broad global scales and finer scales within country or administrative boundaries where more finely resolved biodiversity and climate data are available. Third, the approach is interpretable: analog modeling approaches provide intuitive means for contextualizing and visualizing climate change impacts 58 . Finally, analog approaches can reduce barriers to use because analog locations can be provided to end users who can use these along with gridded data to execute impact projections based on simple spatial queries.
Climate analog modeling approaches also have multiple potential strengths for anticipating climate change impacts that are underappreciated and understudied. Analog sites are often spatially proximal to focal locations and thus retain latent spatial information that vary with climate such as disturbance regimes, edaphic properties, and regional species pools. In contrast, other commonly used empirical approaches such as regression start with geographic information (e.g., georeferenced plots), translate this into aspatial information (e.g., climate space), and then transfer model forecasts back onto the geographic landscape. This approach loses information about latent spatial patterns and processes that structure the systems of interest 59,60 . However, the extent to which an analog approach improves forecasting of climate impacts is unknown; to our knowledge there has been no systematic comparison between analog approaches and more traditional means for forecasting impacts.
Lastly, analog approaches explicitly address the extrapolation problem associated with novel climates. Recent studies 61 suggest that novel climates may be more prevalent and spatially extensive than previously appreciated. Novel climates challenge our capacity to forecast impacts because relationships identified under current conditions cannot be extrapolated into regions with novel conditions 62 . Analog approaches explicitly identify areas with non-analogous conditions and thus highlight where caution should be applied when interpreting forecasts. By contrast, other forecasting methods like regression or machine learning 63 do not inherently guard against extrapolation into novel conditions.
Our approach to modeling climate change-induced ecoregion shifts relies on multiple assumptions. First we assume that ecoregions are effective surrogates for ecosystem-level biodiversity features. It is clear that individual species within an ecoregion may have divergent responses to climatic changes, and thus ecoregions may change as patterns of biodiversity change. However, the ecoregion concept was developed to provide a global classification of spatial units that affords efficiency and consistency at multiple spatial scales and can be used in regional planning-for example, to provide an ecologically coherent boundary for regional planning processes. Further, representation of such spatial units is one of the few measures that can be reliably tracked with current data. Of course, there are multiple ecoregion-based classification systems that can be used. These have been defined based on major environment or habitat types as well as biogeographic discontinuities 64 . We chose the system used by Dinerstein 11 not only because it is commonly used, but specifically because it has been used in proposals for area-based protection targets 7,11 . Moreover, Smith et al. 65 found that range boundaries of~10,000 species were correlated with the ecoregion boundaries of this classification system suggesting that they effectively represent not only habitat types but also species distributions in aggregate. Nevertheless, like all coarse-filter targets, ecoregions are imperfect surrogates for the underlying ecological communities and processes they represent. Moreover, the extent to which ecoregions-or any coarse-filter planning unit-captures anticipated shifts in the distribution of species and habitats remains unknown. Given their wide use, we use ecoregions here, but we expect that we would find similar climate driven shifts in representation in other similarly scaled conservation features such as threatened ecosystems 9 or key biological diversity areas 15 . Our intent is to highlight the potential for novel geographic patterns in the planet's biodiversity regardless of the particular set of polygons used as planning units and to develop a flexible means to anticipate changes to a given portfolio of conservation features however they are defined.
Second, ecoregions can be associated with an a priori set of temperature and water balance variables. Climatic water balance is widely recognized as a driver of the distribution of plant physiognomic types at continental scales and species at local scales 40,66 . Similarly, temperature is a critical predictor of multiple taxa 39 . Nonetheless, individual ecoregions can span broad climatic gradients, and other non-climatic factors like biogeographic history shape present-day biota within ecosystems. Related, like many empirical approaches, analog projections assume climate and biodiversity patterns are in equilibrium. However, rates of climate change in many regions are likely to exceed the ability of plants to keep pace 67 as evidenced by disequilibrium dynamics in many plant systems 68 . Compositional and structural changes may be slow and steady driven by incremental and idiosyncratic changes in species distributions or they may exhibit rapid threshold responses catalyzed by disturbance 69 . Indeed, intensifying disturbances (e.g., wildfire) may catalyze changes in existing biodiversity patterns 70,71 and bring biotic communities more into equilibrium with current and future climate patterns.
Despite these shortcomings, a spatial climate analog approach focused on ecoregions is an effective means for predicting broadscale changes in biodiversity patterns. Given interest in ecoregion-and ecosystem-based approaches for conserving biodiversity, future research should focus on refining and assessing the approach.

Conclusions
There is increasing evidence that many of earth's terrestrial ecosystems are already on a trajectory of extensive compositional and structural changes due to climate change. The global protected area network is our frontline defense against these changes and the potential loss of biodiversity associated with human impacts more broadly. And yet, we show here that climate changes by mid-century may substantially shift the degree to which ecoregions and other spatial planning units remain protected in the global PA network. This is critical to acknowledge, as these planning units are the basis for quantifying representation in multiple PA-expansion schemes proposed in anticipation of upcoming renegotiations of global conservation targets.
An important first step for strategically expanding the PA network is to identify those portions of the planet's terrestrial ecoregions and PAs that will remain compositionally stable and those that are poised for change. Conservation efforts within stable or transitioning areas could have different objectives: the former could focus on conserving large-scale refugia whereas the latter could focus on promoting climatic connectivity between stable areas. In both cases, we must accommodate the strong likelihood of reshuffling under changing conditions. Moreover, we must look beyond PAs because the intervening matrix, especially semi-natural or working landscapes, will also be critically important for protecting and supporting the movement of biodiversity under climate change 72 . Striking a balance between a focus on existing assets and being forward-looking is critical. Above all, we must deploy strategies that can be readily implemented, as the effects of climate change continue to unfold. Indeed, the unprecedented conservation challenges we face require that we enhance the resilience of the global PA network, and all lands, to the emergence of novel patterns in the planet's biodiversity.

Materials and methods
Reverse climate analogs (e.g., backward analogs 52 ) use space-for-time substitution to identify present day sites (herein 'analogs') that are analogous to the future climatic conditions projected for a focal location. To identify analogs, we characterized future climate conditions for focal cells (i.e., all terrestrial pixels) and reference climate conditions for candidate analogs.
Reference period climate data. We characterized climate with four biophysically relevant variables: average minimum temperature (Tmin), warmest monthly average maximum temperature (Tmax) as well as two metrics that characterize the surface water balance: annual cumulative actual evapotranspiration (AET) and climatic water deficit (CWD). We considered these variables to be a minimal set which are broadly relevant for ecosystem properties such as primary productivity 73 , plant distributions 39,40 , and disturbance regimes 74 .
Gridded summaries of AET, CWD, Tmin, and Tmax were acquired from TerraClimate 75 . TerraClimate is a relatively high-spatial resolution (~4-km) global dataset of monthly climate and water balance variables covering the period from 1958-present. We used a 30-year baseline period of 1961-1990 to describe reference conditions during which the anthropogenic climate change signal was more subdued relative to a more recent 30-year period 22 .
Future climate data. Future climate conditions correspond with global mean temperatures that are 2°C above pre-industrial levels per policy targets ("+2°C"). While modeled projections of climate change vary across models and policyemission pathways, many estimates suggest that global mean temperatures will exceed 2°C above pre-industrial conditions by mid-century 76 .
The 2°C scenario was developed using a pattern scaling approach that allows one to express modeled rates of regional change with respect to the global average 77 . The pattern scaling approach allows for the explicit use of observed variability and covariance across variables from TerraClimate, while also accommodating changes in statistical properties of means and variance from global climate models (GCMs) 78 . We developed pattern scalings on a monthly basis for maximum and minimum temperature, specific humidity, 10-m wind speed, precipitation, and surface downward shortwave radiation from 23 CMIP5 climate models for two 30-year time periods, a pre-industrial period (1850-1879) using a historical forcing experiment and a 2070-2099 period using the RCP8.5 forcing experiment. For changes in means, we considered absolute differences for all variables except precipitation which used relative differences (i.e., absolute change divided by pre-industrial average). The pattern scaling for changes in interannual climate variability (ICV) used the same protocol that accounts for changes in the standard deviation. Pattern scaling was calculated separately for each model; but herein we used the 23-model median. Further details on the pattern scaling procedures can be found in Qin et al. 78 .
To account for the fact that we use observations over the period 1961-1990, we multiplied scaling coefficients by 1.65 given that both observations and climate model experiments show that global mean temperature for the 1961-1990 period was 0.35°C above pre-industrial levels.
Identifying analogs. We use Mahalanobis distance (D) and its standardization based on the Chi distribution 61 for quantifying climatic dissimilarity (distance) between focal cells and candidate analog sites within a 2000 km radius (Fig. 1). We chose a 2000 km threshold to balance computational requirements, to limit analog searches within continents, and because it exceeds dispersal rates that can be expected for most terrestrial biota over the study period [79][80][81] . Mahalanobis distance scales multivariate mean climate conditions between a focal cell and those within a search radius (i.e., candidate analogs) by the local magnitude of interannual climate variability (ICV). Because we identify reverse analogs (versus forward analogs), we use annual data from the future period to calculate the future ICV and mean conditions of each focal cell. We note that our approach based on gridded data will likely underestimate ICV compared to station-based observations 61 .
To account for data dimensionality (number of variables), we standardized D using the Chi distribution to calculate sigma dissimilarity, (σ d ; 61 ). σ d can be interpreted as a multivariate Z-score which quantifies the similarity between the future climate projected for the focal cell under a +2°C scenario and the reference climate at each candidate analog. Inferences can be made from the percentiles of the Chi distribution (i.e., 1σ, 2σ, and 3σ for~68th, 95th, and 99.7th percentiles, respectively). The analog with the lowest value of σ d is the closest match to the focal cell in climate space. However, for the dataset we distribute with this paper, we retained the 100 analog sites with the lowest σ d values for each terrestrial focal cell. For the online interactive tool (plus2C.org), we display analog sites with the single lowest σ d value identified within the search radius. Lastly, we considered analogs with σ d > 2 to represent sites that will experience regionally novel climates.
Ecoregions, protected area data, and human modification data. We used global ecoregion data from Dinerstein et al. 11 which was derived from Olson et al. 82 . This dataset consists of 847 ecoregions organized by 14 biomes. We obtained a geospatial dataset representing protected areas from the World Database on Protected Areas (WDPA). We included all protected areas identified as IUCN (International Union of Conservation for Nature) Management Categories I-VI and excluded protected areas categorized as 'marine' or 'proposed'. A large number of protected areas identified by the WDPA were not assigned an IUCN category (i.e., designated as 'Not Reported') but have high levels of conservation (e.g., Kruger National Park in South Africa). Therefore, we included additional protected areas where the level of human modification was similar or less than that observed within IUCN category I-VI protected areas. Human impacts were estimated using the human modification gradient (HMG), a gridded raster dataset (resolution = 1 km) representing cumulative human modification to terrestrial lands 46 . We calculated the mean HMG for each IUCN category I-VI protected area, identified the 80 th percentile of this distribution, and included un-assigned protected areas in our study that had a mean HMG less than or equal to the identified threshold. This procedure increased the protected lands included in this study by 30% (compared to that of strictly using IUCN I-VI according to the WDPA). We converted protected area polygons to a binary raster (protected and not protected) at the resolution of TerraClimate data (~4 km). We retained only cells >75% protected. Each pixel in our analyses was thus associated with a specific ecoregion and biome and was characterized as either protected or not protected. We then further categorized pixels as 'intact', if their HMG values were 'low' 46 (HMG < 0.1).
Projecting shifts in ecoregions. For each focal pixel, we determined the terrestrial ecoregion that was most likely to be associated with its future climate conditions. We did so by assigning, for each focal pixel, the plurality vote of ecoregions associated with the 100 top analogs with σ d < 2 (Fig. 1). Focal cells that retain their current ecoregion class were considered stable whereas those that transitioned to other ecoregions were noted. We did not assign analogs for areas with insufficient data, over inland water bodies, regions of rock and ice with no annual AET, and Tropical Areas where the ICV is too low to effectively calculate σ d . Areas of insufficient data comprised 14.1% of the total land area, with the majority in Greenland and Antarctica. Accordingly, we did not consider the current or future ecoregions within these areas in our analyses. For this reason, we used 803 current ecoregions. Finally, ecoregions within biomes that are hydrologically (e.g., Flooded Grasslands and Savannahs, Mangroves), rather than climatically defined, were included in global calculations, but excluded from ecoregion and biome level analyses leaving 759 ecoregions in the reference period.
Uncertainty in projections. One source of uncertainty in climate projections is climate sensitivity among GCMs and climate forcing scenarios. The climate data used in this study differs from typical climate projections used in studies which are tied to a specific future time-frame. By tailoring our projections to an amount of warming of the global mean temperature as opposed to a specific future time period, we limit this source of uncertainty leaving regional variability and internal climate variability as the predominant sources of uncertainty. Our approach employs a multi-model median to provide a central estimate of change. This approach minimizes biases originating from individual models-yet also does not explicitly allow for an assessment of climate model uncertainty. Inter-model variability in scaling factors is evident across variables, months, and locations. This is consistent with the different climate sensitivities of models and how they manifest variations in atmospheric circulation (that then impact precipitation, humidity, radiation). To help constrain this uncertainty, we provide additional information on variability in inter-model scaling factors in Supplementary Fig. S3.
In addition to climate model uncertainty, uncertainty arises based on our use of spatial climate analogs. We use two approaches to characterizing the uncertainty of these projections. The first relies on the dissimilarity in climate (σ d ) between projected future conditions for focal locations and their candidate analogs. For instance, sites with low σ d values have analog locations that are climatically similar to conditions projected for a focal location whereas high σ d values suggest dissimilar conditions (Fig. S1). The second approach to characterizing uncertainty relies on the plurality of "votes" used to assign ecoregions to a given focal cell. We assume greater confidence in a projected ecoregion shift if the plurality is large (i.e., if a large proportion of the 100 reverse analogs agree on the same projected change) (Fig S2).
Representation and conservation targets. PA area is neither gained nor lost within our analyses. Thus, changes in percent protection for a given ecoregion reflect the redistribution of a fixed amount of protected area-either a gain or a loss -between different ecoregions. To assess the impact of shifts in ecoregions on conservation targets, we quantified the percent of each ecoregion that was protected (percent protected; PP e ) under reference and future conditions. To examine how PP e changes under warming scenarios, we calculated ΔPP e , which is the change in PP e for a specific ecoregion between the reference and the +2°C scenario. When aggregating PP e values for a biome, we reported means using the variable PP. Following this approach, ΔPP represents the mean ΔPP e for the aggregate of ecoregions. This approach is not area-weighted to ensure changes in geographically small ecoregions are not discounted.
If climatic conditions disappear-either regionally or globally-ecoregions may be "lost". Consequently, the total number of ecoregions under reference conditions (759) exceeded the number of ecoregions in the +2°C scenario (703). ΔPP e values were therefore only defined for ecoregions that were represented in both reference conditions and the future scenario.
To assess whether ΔPP values significantly differ from 0, we used a bootstrap method based on sampling with replacement ecoregions within each biome for the current and future climate scenario. The number of resamples for a given biome is equal to the number of ecoregions in the biome. We replicated the sampling procedure 10,000 times to create an empirical distribution of the difference between the sample means and report the 10th and 90th percentiles of the distribution.
Focal sites with no-analog climate conditions complicate ecoregion assignments. We did not include them in PP calculations. No-analog climates are distinct from disappearing ecoregions in that the former was caused by nonanalogous climate conditions (σ d > 2) and hence represented an undefined reverse analog for a specific location. Disappearing ecoregions signify that a specific ecoregion was not found to be a reverse analog for any location. We considered pixels with no-analog conditions to have an unknown future that could not be assessed in the context of our analysis.

Data availability statement
Global ecoregion data used in this study was downloaded from https://ecoregions. appspot.com/. Terraclimate data is available at http://www.climatologylab.org/ terraclimate.html, and Human Modification Gradient data was accessed from https:// doi.org/10.6084/m9.figshare.7283087.v1. A netcdf file of latitude, longitude, and σ d values for 100 reverse climate analogs for global terrestrial lands is available at https://climate. northwestknowledge.net/ACSL/analog2c.tar. The expanded protected area database used in this analysis and an ecoregion summary table of changes in area and protection under a +2°C scenario can be found at https://doi.org/10.6084/m9.figshare.15106233.