Toward a climate‐informed North American protected areas network: Incorporating climate‐change refugia and corridors in conservation planning

Global and national commitments to slow biodiversity loss by expanding protected area networks also provide opportunities to evaluate conservation priorities in the face of climate change. Using recently developed indicators of climatic macrorefugia, environmental diversity, and corridors, we conducted a systematic, climate‐informed prioritization of conservation values across North America. We explicitly considered complementarity of multiple conservation objectives, capturing key niche‐based temperature and moisture thresholds for 324 tree species and 268 songbird species. Conservation rankings were influenced most strongly by climate corridors and species‐specific refugia layers. Although areas of high conservation value under climate change were partially aligned with existing protected areas, ∼80% of areas within the top quintile of biome‐level conservation values lack formal protection. Results from this study and application of our approach elsewhere can help improve the long‐term value of conservation investments at multiple spatial scales.

species and their underlying ecosystems become moving targets (Hannah, Midgley, & Andelman, 2007;Lovejoy & Hannah, 2019). Meanwhile, naturally functioning ecosystems are critical for mitigating greenhouse gas emissions via carbon sequestration (Dass, Houlton, Wang, & Warlind, 2018;Seddon, Turner, Berry, Chausson, & Girardin, 2019). Thus, substantial increases in the protection of natural ecosystems above Aichi Target levels can both mitigate emissions via carbon sequestration and respond to the combined effects of climate and land-use change (Dinerstein, Vynne, & Sala, 2019). Although recent surveys have shown public support for ambitious increases in protected area targets (Wright, Moghimehfar, & Woodley, 2019), systematic conservation planning approaches are needed to identify the most efficient investment of conservation and restoration resources.
Indicators of ecological exposure to climate change can be used to identify areas of high climate-change vulnerability (e.g., Batllori, Parisien, Parks, Moritz, & Miller, 2017), but they can also be used as inverse indicators of ecosystem persistence, aka refugia potential (Carroll, Roberts, & Michalak, 2017). To assist climate-informed conservation planning, several such metrics have been proposed and mapped across large areas of North America, enabling the explicit consideration of climate-altered futures within systematic conservation planning frameworks (Belote, Carroll, & Martinuzzi, 2018;Carroll et al., 2017). Climatic macrorefugia represent areas where climate-threatened species are more likely to persist in place or can more readily disperse to within decadal timeframes (Michalak, Lawler, Roberts, & Carroll, 2018;Stralberg et al., 2018). New macrorefugia indices are based on the concept of climate velocity, that is, the speed at which an organism must migrate to keep pace with climate change (Hamann, Roberts, Barber, Carroll, & Nielsen, 2015;Loarie et al., 2009) and can be derived from climate analogs Michalak et al., 2018) or species' niches (Stralberg et al., 2018). These "climate-tracking" metrics (Michalak, Stralberg, Cartwright, & Lawler, 2020) inherently encompass two simpler climate exposure metrics for a given spatial location: the magnitude of projected climate change and the slope of surrounding climatic gradients.
However, macrorefugia metrics are limited by the coarse resolution of global climate models (Ashcroft, 2010). Thus, it is also important to consider microrefugia generated by climatically buffered terrain features such as valley bottoms and north-facing slopes (Dobrowski, 2011). Although microrefugia potential is difficult to quantify, metrics of environmental diversity can serve as surrogates (Ackerly, Loarie, & Cornwell, 2010;Lawler, Ackerly, & Albano, 2015).
Moreover, a focus on refugia could exclude species with distributions that are geographically limited, fragmented, or characterized by high levels of climate exposure. Thus, it is also important to identify areas that may serve as climate "corridors" to facilitate long-distance movement (McGuire et al., 2016). Although connectivity in a climate-change context can be assessed in a number of different ways (Keeley, Ackerly, & Cameron, 2018), network theory principles are widely used to estimate climate corridors by delineating paths between current and future distributions of a species or specific climate types (Carroll, Parks, Dobrowski, & Roberts, 2018;Littlefield, McRae, Michalak, Lawler, & Carroll, 2017). Predictions from such models have been shown to correlate with dispersal routes and gene flow (McRae & Beier, 2007).
One way to assess multi-species refugia or corridor potential is to develop composite indices based on weighted sums of results from individual species (Stralberg et al., 2018) or climate types . However, at a continental scale this leads to a strong emphasis on mountainous regions, given the low velocity and high environmental diversity associated with steep gradients in terrain. Furthermore, species-specific refugia depend on environmental niche characteristics not generally captured by combined metrics (Michalak et al., 2020). Thus, systematic conservation planning tools, such as the widely used Zonation software (Moilanen, 2007), are necessary to consider species complementarity and regional priorities.
We applied the Zonation algorithm to identify priority conservation areas across the continental United States and Canada based on indicators of macrorefugia, microrefugia, and climate corridors. Our primary goals were to (a) identify areas that constitute the most efficient and complementary candidates for broad-scale land conservation in the face of climate change, and (b) assess how well the current protected areas network represents these conservation priorities in different regions. Secondarily, we aimed to analyze how priorities change geographically and quantitatively when considering species-specific refugia, climate corridors, and microrefugia potential; when stratifying by biome; and when accounting for human development. To accomplish this, we developed an approach that considers complementarity of multiple conservation objectives, explicitly incorporates climate exposure for a diversity of species, and provides a relatively high-resolution analysis over a broad spatial extent. Results from this study and application of our approach elsewhere can help improve the long-term value of conservation investments at multiple spatial scales.

METHODS
Our study area consisted of the continental United States and most of Canada (∼16.1 million km 2 ). We omitted portions of North America in the far north (Canadian high Arctic) and south (Mexico, Central America, and the Caribbean islands) based on data availability. We conducted our analysis at a 5-km grid cell resolution and used a Lambert Azimuthal Equal Area projection.
Our objective was to rank individual grid cells according to their combined conservation value with respect to five classes of complementary, climate-informed conservation objectives (Figure 1). In doing this, we also considered constraints imposed by human development, and evaluated priorities at the scale of (a) our entire northern North America study area and (b) individual biomes within this study area. Rather than generating scenarios consisting of all possible combinations of objectives, constraints, and scales, we started with a full scenario considering all seven factors and then evaluated the contribution of each individual factor by generating seven additional scenarios, each omitting one of the factors. To quantify the contribution of each individual factor to overall rankings, we calculated and averaged grid cell-level differences between the full scenario and each of the seven additional "leave one out" scenarios.
We used the Zonation software, which ranks the value of each cell based on user-defined criteria and iteratively removes cells with the lowest proportional value across all criteria (Moilanen et al., 2014). The result is a ranking of all cells, with the cells that are last removed having the highest rank. For all scenarios, we used the "core-area Zonation" (CAZ) cell-removal rule, which maximizes core habitat for each individual conservation objective, rather than treating them as interchangeable. To improve optimality of solutions, we specified that Zonation removes a small number of 5-km grid cells from any location within the study area during each iteration (i.e., warp factor = 10). To exclude highly human-modified areas, we used a 2009 built environment layer (Venter, Sanderson, & Magrach, 2016) as an analysis mask. Zonation settings files are provided in Appendix S1.
Our first conservation objective was a species-neutral index of macrorefugia potential based on the inverse of logtransformed backward climate velocity  for 3,700 unique multivariate climate types as described in Carroll et al. (2017) (Figure 1a). Climate types, each on the order of 100-1,000 km 2 , generally represent regions smaller than an average species' range, and thus could represent locally adapted genotypes or range-restricted endemic species. The index was based on the ensemble mean of 15 representative global climate model (GCM) projections (averaged over up to five multiple runs) for end-of-century (2071-2100) climate conditions given continued high greenhouse gas emissions, as defined by representative concentration pathway (RCP) 8.5 (Wang, Hamann, Spittlehouse, & Carroll, 2016).
The second and third conservation objectives consisted of macrorefugia potential based on current species' niches for 268 songbird species (Figure 1b) and 324 tree species (Figure 1c). We considered these partial surrogates for other relatively wide-ranging terrestrial species, recognizing that the representation of many species is likely inadequate (Rodrigues & Brooks, 2007). The species-specific refu-gia index, also derived from backward climate velocity, is a negative-exponential function of the straight-line distance from areas of end-of-century climatic suitability (RCP 8.5) to the nearest area of current climatic suitability, averaged across four representative GCMs (Stralberg et al., 2018). We first used Zonation with the CAZ cell-removal rule to rank conservation value across the study area for all tree and songbird species combined, and used that as an input to the combined scenario. This allowed us to consider requirements of all species individually, which translates into an emphasis on future rare species, effectively down-weighting the importance of species with large areas of high refugia value.
Recognizing that many species and subspecies with restricted ranges will have to migrate large distances to keep pace with climate change, our fourth conservation objective was climate corridors (Figure 1d): areas through which multiple species will need to migrate to reach future suitable climate space . We used climate corridors for the same 3,700 climate types, using a current flow centrality metric based on circuit theory, as described in Carroll et al. (2018). Although macrorefugia indices were based on the straight-line distance between current and future climate analogs, climate corridors follow more circuitous paths that minimize exposure to climatically hostile areas (Dobrowski & Parks, 2016). Again, we used results from separate corearea Zonation analyses as inputs to the combined scenario.
All of the macrorefugia and corridor metrics that we included have been shown to be robust to the choice of GCM Carroll et al., 2018;Stralberg et al., 2018). Although we focused on an end-of-century, high-emission climate-change scenario (RCP 8.5), the conservative assumptions of the refugia and corridor metrics make them broadly applicable across a range of less extreme change scenarios.
Our fifth and final conservation objective was environmental diversity (Figure 1e), which serves as a proxy for microrefugia potential in that topographically heterogeneous landscapes support a range of microclimate conditions, some of which are decoupled from regional climates (Lawler et al., 2015). Specifically, we used an elevational diversity metric based on the mean elevation difference between all pairs of cells within a given spatial neighborhood (Ackerly et al., 2010), in this case within a 3-km moving-window based on a 100-m digital elevation model generalized to 1 km .
Pairwise correlations (Pearson's R) among these five input layers were generally low, ranging from -.09 between tree macrorefugia and climate corridors to .55 between climatetype macrorefugia and elevational diversity (Table S1).
In terms of scale, we imposed geographic stratification of conservation priorities by biome using Zonation's 'administrative units' feature. Biomes were defined as level 1 ecoregions (Figure 1) as mapped by the Commission for Environmental Cooperation (CEC, 1997) and were weighted  ; (b) tree macrorefugia (Stralberg et al., 2018); (c) songbird macrorefugia (Stralberg et al., 2018); (d) climate corridors ; and (e) elevational diversity ). (f) A 2009 human development index (Venter et al., 2016) was used as a cost layer equally to ensure proportional distribution of priority areas (see Appendix S1 for Zonation settings files).
Finally, as a constraint we added a 2009 human development index (Venter et al., 2016) with scores ranging from 0 (no development) to 50 (most developed) as a cost layer in Zonation (Figure 1f). We squared the development index to increase the contrast between and high-and low-development areas, divided it by 100, and made it nonzero by adding 1, resulting in an index ranging from 1 to 26.
For the full, biome-stratified scenario, as well as the nonstratified version of the full scenario, we summarized Zonation rankings within the protected areas network for each biome. We defined protected areas according to the IUCN's definition and included all categories (I-VI). To obtain the most comprehensive and up-to-date map of protected areas, we merged layers from the Canadian Protected and Conserved Areas Database 1.x (2019), PAD-US Conservation Biology Institute 2.1, and the Commission for Environmental Cooperation (CEC, 2017). We converted protected area polygons to 5-km rasters based on a maximum combined area assignment in ArcGIS 10.6.1. Spatial analyses were conducted in R version 3.5.2 (R Core Team, 2018) using the "raster" package (Hijmans, 2019).

RESULTS
Zonation-based conservation rankings for the biome-stratified full scenario were well distributed across the study area, by definition (Figures 2a and S1a). Comparing this full scenario with those in which one factor was omitted, mean pixel-level differences ranged from 0.019 for climate-type macrorefugia to 0.110 for biome stratification, with various patterns of positive (blue) and negative (red) influences (Figure 2b-h). Climate-type macrorefugia had the largest influence in the North American Desert and Northwestern Forested Mountains biomes (Figures 2b and S1b). Songbird-driven increases were primarily in the Northern Forests and Great Plains (Figures 2c and S1c), whereas tree-based increases occurred mostly in the Northern Forests and Eastern Temperate Forests (Figures 2d and S1d). Climate corridors resulted in the largest gains through the center of the continent (Figures 2e and S1e), and elevational diversity was responsible for gains mostly in western mountains (Figures 2f and S1f). Finally, stratification by biome yielded the largest gains in the Eastern Temperate Forests (Figures 2g and S1g), whereas consideration of human development cost led to largest shifts within the Great Plains and Eastern Temperate Forests (Figures 2h and S1h).
Highest priority areas identified in the full scenario included several existing large national parks, but also many large areas with no current protection, especially within the Taiga, Northern Forest, Hudson Plain, and Northern Great Plains biomes (Figure 3).
With respect to conservation values within protected areas, the full scenario without biome stratification resulted in highly variable results across biomes (Figures 4a and 4b). In general, western biomes had the highest conservation value within protected areas, with median values much higher than would be expected by chance (0.5). However, all except the Hudson Plain biome had higher conservation values within versus outside protected areas. With biome stratification, the median conservation value varied less across biomes, and was greater than or equal to 0.5 for all biomes except for the Hudson Plain (Figure 4c). Averaged across biomes, conservation value within protected areas was greater than 0.5 for both full scenarios, but highest for the version without biome stratification (Figures 4b and 4c).
Focusing on the top-ranked quintile (20%) of cells (∼3.2 million km 2 ), 26.0% (∼0.84 million km 2 ) was contained within the current protected areas network under the scenario without biome stratification (Table S2), whereas 18.9% (∼0.61 million km 2 ) of top-quintile conservation value was protected under the biome-stratified scenario (Table 1). Comparing individual biomes within the biome-stratified scenario, the very small (within our study area) Tropical Wet Forests biome had the largest percent of high-value conservation areas protected (80.4%), followed by the Marine West Coast Forest biome (45.3%) ( Table 1). The biomes with the smallest percent of high-value conservation areas protected were the Great Plains (4.4%) and the Hudson Plain (5.4%). Areawise, the biome with the greatest area of top-quintile protection was the Northwestern Forested Mountains (∼0.16 million km 2 ), followed by the Taiga (∼0.08 million km 2 ). The biomes with the largest unprotected top-quintile areas were the Northern Forests (∼0.52 million km 2 ) and Great Plains (∼0.49 million km 2 ) biomes, followed closely by the Taiga (∼0.45 million km 2 ) and Eastern Temperate Forests (∼0.40 million km 2 ) biomes.

DISCUSSION
In a combined evaluation of multi-species climate-change refugia and corridor potential, we systematically identified areas with high potential for persistence of multiple biodiversity elements given climate change. This builds upon previous refugia mapping efforts that either did not contain information about individual species' requirements  or did not account for complementarity among species and metrics (Michalak et al., 2020).
The inclusion of climate corridors had the strongest influence on conservation rankings, shifting priorities toward central regions with high potential for species to migrate in response to shifting habitat suitability. Rankings were also strongly influenced by the addition of speciesspecific macrorefugia layers, which capture key niche-based   (Stralberg et al., 2018). This improved the biological relevance of resulting conservation priorities, which could be fur-ther refined by including additional taxonomic groups as they become available. The consideration of elevational diversity as a proxy for microrefugia potential resulted in redirection of conservation priorities to better capture terrain differences that may provide local refugia. Taken together, these key biodiversity components can help address the need for protected areas strategies that are representative of biodiversity, not just geographies and landforms, thereby providing greater opportunities to meet and exceed Aichi targets (Dinerstein et al., 2019). By optimizing conservation values for multiple species and local climate types simultaneously, our approach distributed rankings across major biomes and identified specific areas of greatest combined conservation value within each biome. Thus, our results are more regionally relevant than simple combinations of these metrics Stralberg et al., 2018), which may overlook areas of high unique importance to multiple species, especially within high-velocity areas such the Northern Forests and Great Plains.
Areas of highest regional conservation value under climate change were partially aligned with existing protected areas, which had higher conservation value on average than the biome as a whole. This was particularly true in western biomes, where protected areas are naturally biased toward areas of high elevation and steep terrain, largely due to a combination of remoteness, inaccessibility, and high scenic and recreational value, as well as low land-use value (  Pfaff, 2009). Our analysis largely excluded tropical regions, where many effectiveness gaps have been identified with respect to biodiversity representation (Rodrigues, Andelman, & Bakarr, 2004). However, we also found that many large areas of high regional and continental conservation value lack formal protection. Within the highest-ranked 20% of the study area, over 80% of land remains unprotected. Biomes with the lowest rates of protection for high-value refugia and corridors (∼5%) were the Hudson Plain and Great Plains. The former is especially important from a climate corridor perspective, as many northern species' northward-shifting distributions are bifurcated by the Hudson and James Bay division of north-central North America (Murray, Peers, & Majchrzak, 2017;Stralberg, Matsuoka, & Hamann, 2015). It also contains a large portion of global carbon stores (Tarnocai et al., 2009) and high ecological inertia  in peat-dominated wetlands.
The consideration of human development further constrained conservation priorities by emphasizing areas with greater broad-scale ecological integrity. This is especially valuable from the standpoint of identifying large intact conservation areas of continental significance, especially for wide-ranging species. However, recognizing the value of small and fragmented areas for biodiversity (Wintle, Kujala, & Whitehead, 2019), we suggest that rankings not considering human development are useful for identifying subregional conservation and restoration priorities, especially in the context of comprehensive, multi-stakeholder planning processes.
The spatially explicit, climate-informed conservation rankings generated by our systematic assessment may serve as direct inputs to broad-scale conservation planning efforts. However, given the sensitivity of our rankings to each of the individual conservation objectives considered, we caution against a prescriptive interpretation of these results, especially at subregional scales. Rather, we offer a general approach that can be repeated or enhanced-with refined species inputs or in combination with other conservation objectives and constraints.
Finally, although our analysis was focused on identifying areas that promote species persistence, as candidates for conservation and restoration activities, we do not suggest that areas of high vulnerability and turnover potential should be ignored. In such areas, change monitoring and active intervention may be appropriate conservation strategies even within protected areas (Belote, Dietz, & Jenkins, 2017;Gillson, Dawson, Jack, & McGeoch, 2013). Indeed, the development of guidelines concerning the relative importance and geographic distribution of resistance versus transformation strategies will be a key step in increasing the relevance of broad-scale spatial analyses such as ours in conservation planning.

CONCLUSION
Given the complex task of conserving biodiversity in a changing climate, conservation priorities and decisions benefit from systematic, quantitative planning using a suite of climateinformed metrics. Our wide-ranging, comprehensive analysis identified multiple regions of high conservation value with no formal protection. In terms of both achieving percentage area targets and prioritizing areas of high biodiversity importance as part of the Aichi Targets and beyond, these findings can provide valuable guidance and insight for climate-conscious assessments of conservation and management priorities.