Climatic, topographic, and anthropogenic factors determine connectivity between current and future climate analogs in North America

As climatic conditions shift in coming decades, persistence of many populations will depend on their ability to colonize habitat newly suitable for their climatic requirements. Opportunities for such range shifts may be limited unless areas that facilitate dispersal under climate change are identified and protected from land uses that impede movement. While many climate adaptation strategies focus on identifying refugia, this study is the first to characterize areas which merit protection for their role in promoting climate connectivity at a continental extent. We identified climate connectivity areas across North America by delineating paths between current climate types and their future analogs that avoided nonanalogous climates, and used centrality metrics to rank the contribution of each location to facilitating dispersal across the landscape. The distribution of connectivity areas was influenced by climatic and topographic factors at multiple spatial scales. Results were robust to uncertainty in the magnitude of future climate change arising from differing emissions scenarios and general circulation models, but sensitive to analysis extent and assumptions concerning dispersal behavior and maximum dispersal distance. Paths were funneled along north‐south trending passes and valley systems and away from areas of novel and disappearing climates. Climate connectivity areas, where many potential dispersal paths overlapped, were distinct from refugia and thus poorly captured by many existing conservation strategies. Existing protected areas with high connectivity values were found in southern Mexico, the southwestern US, and western and arctic Canada and Alaska. Ecoregions within the Isthmus of Tehuantepec, Great Plains, eastern temperate forests, high Arctic, and western Canadian Cordillera hold important climate connectivity areas which merit increased conservation focus due to anthropogenic pressures or current low levels of protection. Our coarse‐filter climate‐type‐based results complement and contextualize species‐specific analyses and add a missing dimension to climate adaptation planning by identifying landscape features which promote connectivity among refugia.


| INTRODUCTION
In coming decades, much of the earth's biota will experience climatic conditions outside the range to which they are adapted (Mora et al., 2013). Some populations will be able to persist in place in climatic refugia within their current range (Keppel & Wardell-Johnson, 2012).
However, persistence of many populations will hinge on their ability to disperse and colonize habitat which has become newly suitable for their climatic requirements. Opportunities for such dispersal and range shifts may be limited unless areas that play a key role in facilitating dispersal under climate change (which we term "climate connectivity areas" after McGuire, Lawler, McRae, Nuñez, and Theobald (2016)) are identified and protected from land use that may impede movement (Heller & Zavaleta, 2009).
Identifying areas that promote climate connectivity is challenging given the diverse and multifaceted nature of biotic response to climate change, as well as uncertainty as to the pace and pattern of change arising from contrasts between projections based on different atmosphere-ocean general circulation models (AOGCM) and emission pathways. In order to allocate conservation resources to maximize the persistence of species in the face of climate change, planners need robust answers to such questions as: What areas will provide key dispersal routes as climate shifts? Do such routes overlap with areas heavily altered by humans? Where do dispersal routes for different species or routes connecting different climatic niches overlap? Are these areas associated at local, regional, and continental extents with topographic and other features that could help us predict their locations? And finally, will areas that promote connectivity in the near future remain important as climate change strengthens in the more distant future?
The concept of centrality provides a novel and informative approach to addressing such questions. Many habitat connectivity analyses apply models which treat landscapes as graphs, defined in this sense as a networks of nodes linked along edges (Carroll, McRae, & Brookes, 2012;Newman, 2010) (Supporting information Figure S1, Table S1). Edges, which represent functional connections (e.g., dispersal) between nodes, are assigned weights that represent an attribute such as habitat quality. A sequence of nodes connected by edges forms a path. Although abstracted depictions of landscape pattern, graphs may reveal emergent aspects of landscape structure that are not otherwise discernible.
Most graph-based analyses focus on delineating paths between individual source and target patches (Adriaensen et al., 2003).
Betweenness centrality (henceforth centrality) metrics instead evaluate paths between all or a subset of pairwise combinations of sites on a landscape to rank the contribution of each site to facilitating dispersal across the network of sites (Ahuja, Magnanti, & Orlin, 1993;Newman, 2010). Areas of high centrality emerge as conservation priorities because the loss of sites where many paths overlap would disproportionately lengthen dispersal distances or transit times between nodes and thus disproportionately reduce dispersal across the network (Álvarez-Romero et al., 2018;Carroll et al., 2012). This paper develops a novel approach that extends previous work in several ways (Table S2). Initial efforts to identify climate connectivity areas made the simplifying assumption that such areas can be identified based solely on current climate or other environmental characteristics. For example, Nunez et al. (2013) identified climate connectivity areas by connecting patches of intact habitat with other patches that had a cooler climate under current conditions, in order to provide migration routes as climate warmed. In contrast, subsequent work by McGuire et al. (2016) connected source patches to destination patches projected to be similar or cooler under future climates. Littlefield, McRae, Michalak, Lawler, and Carroll (2017) extended this approach by evaluating dispersal potential between climatically analogous current and future patches via paths that avoided human-modified habitat.
Climatic velocity is a metric often used to assess geographic patterns of climate exposure (Loarie et al., 2009). Climatic velocity values are often based on the straight-line distance between a climate type and the nearest location of the analogous climate type under future climates (Hamann, Roberts, Barber, Carroll, & Nielsen, 2015;Ordonez & Williams, 2013). Because the assumption of straight-line dispersal is often unrealistic for terrestrial flora and fauna, Dobrowski and Parks (2016) developed a velocity metric based on the shortest or least-cost path between current and future climate analogs that avoided nonanalogous climates, which is typically longer than a straight-line path. The rationale for this approach is similar to that underlying climatic niche models: if climate is a key habitat factor for an organism, climatically hostile areas that deter movement to newly suitable habitat will threaten persistence under climate change.
In this study, we developed a centrality metric which combines the climatic-resistance-based approach of Dobrowski and Parks (2016) with methods of evaluating dispersal potential similar to those of Littlefield et al. (2017). We identified priority areas for connectivity under climate change by connecting current and future climate analogs along paths that minimize exposure to nonanalogous climates (Dobrowski & Parks, 2016). Centrality based on climatic resistance has previously been used to identify key dispersal areas for marine systems (García Molinos et al., 2017). This study represents the first application of climatic-resistance-based centrality to terrestrial systems under climate change. We synthesized results from all climates occurring in North America to identify areas of high centrality whose loss could potentially impact many organisms, and to discern what aspects of continental and regional topography help explain broad-scale connectivity patterns.
We evaluated connectivity across the landscape as a whole rather than solely between patches with low degree of human modification. We compared priorities identified by two methods (shortestand current-flow-based paths) that make contrasting assumptions about the ability of dispersers to find and follow optimal paths (Carroll et al., 2012). We compared connectivity areas important in the near future to those important in the more distant future to assess whether priorities are robust to time horizon. Finally, we assessed the vulnerability of key climate connectivity areas by CARROLL ET AL.
| 5319 evaluating the extent to which they fall within existing protected areas, or conversely within areas where high human land use intensity might impede movement. Our approach, which can be applied at a range of scales, avoids many of the simplifying assumptions of previous studies, allowing planners to more realistically assess where conservation resources should be directed to protect climate connectivity areas and facilitate persistence of species under climate change.

| Climate data
Climate data used here cover the North American continent (excluding areas south of Mexico), a study area which includes a diverse set of ecoregions that vary widely in their topographic, climatic, and biotic attributes. To represent current climate, we used modeled estimates of 1981-2010 monthly temperature and precipitation (Lambert Azimuthal Equal Area projection) as derived by the ClimateNA software (Wang, Hamann, Spittlehouse, & Carroll, 2016).
ClimateNA data at 1 km resolution were resampled to 5 km resolution to increase computational feasibility as described below. Cli-mateNA current climate data are in turn based on data developed by the PRISM project, which used regression to interpolate weather station data based on location, elevation, coastal proximity, topographic facet orientation, vertical atmospheric layer, topographic position, and orographic effectiveness of the terrain (Daly et al., 2008).
Projections of future climate (as monthly mean temperature and precipitation) were derived from the Coupled Model Intercomparison Project phase 5 (CMIP5) database which includes AOGCMs corresponding to the 5th IPCC Assessment Report (Wang et al., 2016).
For the main analysis, future projected temperature and precipitation were calculated as anomalies from the current (1981-2010) reference period to the future 2071-2100 period (hereafter "2080s"), based on an ensemble mean of 15 representative CMIP5 AOGCMs for representative concentration pathway (RCP) 8.5, a "business-asusual" scenario characterized by relatively high greenhouse-gas emissions (Wang et al., 2016). We also compared results from the main analysis with results from additional future time periods and RCPs as part of the sensitivity analysis described below. Anomaly grids were downscaled via local regression and the difference was added to the baseline climate normal data to arrive at the final climate surface (Wang et al., 2016).
From the monthly temperature and precipitation gridded datasets, we developed 11 bioclimatic variables (Table S3) that were more directly related to ecological factors (Wang et al., 2016). We used principal components analysis (PCA) to reduce the dimensionality of the 11 climate variables, which reduced collinearity and increase computational feasibility in subsequent stages of the analysis. Following Carroll et al. (2017), we used climate data based on the PCA's first 2 axes (hereafter PC1 and PC2), which explained 89.3% of variance (Table S3).

| Centrality metrics
We adapted the methods of Dobrowski and Parks (2016) to delineate paths between current and future climate analogs that jointly minimized geographic distance and exposure to increasingly dissimilar climate based on a cost penalty parameter. However, we used multivariate climate surfaces rather than mean annual temperature as did Dobrowski and Parks (2016). We defined 100 equal-width climate types for each of the first two PCA axes. Approximately one third of the 10,000 potential unique combinations of PC1 and PC2 occur under current North American climates. Previous sensitivity analyses (Carroll, Lawler, Roberts, & Hamann, 2015;Hamann et al., 2015) suggested that classifications resulting in several thousand types represent a good balance between the extremes of loss of information due to excessively broad types vs. domination of the results by disappearing or no-analog climates due to excessively narrow types. An approach based on equal-width climate types does not incorporate information on the historical range of local climatic variability at a site (Mahony, Cannon, Wang, & Aitken, 2017); incorporating such information would represent a useful future extension of our method.
For all pixels of each current climate type (i.e., unique combination of PC1 and PC2) in North America, we identified all pixels with the same climate type in the future time period. A pixel is used here to represent a unit of geographic space, as distinguished from a climate type, which represents a subdivision of climatic space (Table S1). We developed cost surfaces (one for each climate type) in a manner that accounts for spatiotemporal changes in climate, allowing isoclines to be traced over landscapes and through time. To do so, we first interpolated current and future PC1 and PC2 grids to represent incremental changes over the current to future time period (Supporting information Figure S2). These interpolated grids represent incremental changes in multivariate climate conditions between current (1981-2010) and future period (2071-2100) and assume a linear trend in changes in climate over this time period.
We created cost surfaces (Supporting information Figure S2) for each of these interpolated climate grids such that costs were proportional to the level of dissimilarity between PCA values of the pixel of interest under current climate and all other pixels. We adapted the methods of Dobrowski and Parks (2016) to our multivariate approach by adding the dissimilarity values from the two PCA axes and multiplying by a cost penalty value of two dimensionless units such that cost increased by two for each increase in dissimilarity of one scaled PCA units (see Supporting information Text S1 for code). Lakes and ocean were arbitrarily assigned a cost value of 5,000 to force trajectories to avoid water when possible. We then calculated, among all intermediate cost surfaces, the minimum cost for each pixel; this minimum cost is used to generate the final cost surface for each climate type.
By basing resistance on the minimum cost for the entire period, we sought a computationally efficient approximation that would allow trajectories between source and destination locations to follow the climate as it changes and avoid traversing dissimilar climates. For example, assuming a source pixel in the southern Great Plains region and a destination pixel in the northern Great Plains region with little intervening topographic variation, an organism following the minimum cost surface could shift at the same rate as climate changed without experiencing dissimilar climates.
Although the minimum-resistance approach provides a good approximation in flatter terrain, it is less accurate in topographically complex terrain where there may be a temporal mismatch between the timing of potential range expansion through a given pixel and the resistance value assigned to that pixel. For example, an area of similar climate (and hence low resistance) might be present in a distant valley in the near future, but minimum resistance would not accurately reflect the climate experienced by an organism when that area was transited in the distant future, and might identify overly simplified (i.e., straighter) paths.
An alternative approach explicitly considers the temporal ordering of climate-based resistance in a step-wise fashion, by identifying paths for each intermediate timestep and linking these paths into continuous chains across the entire time period of interest using network flow algorithms (Alagador, Cerdeira, Araújo, & Saura, 2014;Graham, Vanderwal, Phillips, Moritz, & Williams, 2010;Phillips, Williams, Midgley, & Archer, 2008). However, such algorithms are currently too computationally intensive for analysis at the extent and resolution considered here (Alagador, Cerdeira, & Araújo, 2016;Carroll et al., 2012). In addition, our approach implicitly assumes that organisms are adapted to the current climate at a location, and any departure from the current climate type should be avoided. We do not consider factors such as competition or the relationship of local climate to the species climatic niche as a whole, which could result in an organism experiencing increased fitness under altered climate.
For each climate type, a raster of centrality values was calculated using the set of source and target pixels for each climate type identified previously. Values derived from this type of input are termed subset centrality because source and target nodes are a subset of the network as a whole, although in total across all subsets (climate types) they encompass the entire network (Carroll et al., 2012). We summed centrality results across all climate types in order to give equal weight to all pixels. An alternate approach which weighted common and rare climate types equally irrespective of their geographic extent gave similar results (Spearman's rank correlation (ρ) of centrality values = 0.93).
We evaluated centrality using two contrasting types of graphbased algorithms (shortest-path and current flow; Newman, 2010).
The primary approach used current-flow algorithms based on electrical circuit theory that treat landscapes as conductive surfaces, that is, networks of nodes connected by resistors. These models are widely applied in ecological research and conservation planning, and model predictions have been shown to correlate with dispersal routes and gene flow (McRae & Beier, 2007). When current is injected into a source node and allowed to flow across a network until it reaches a target node, the amount of current flowing through each intermediate node reflects the likelihood that a random walker leaving the source node, and moving along edges with probabilities inversely proportional to their resistance weights, will pass through the intermediate node on its way to the target node.
Current-flow-based methods assume that dispersers have limited knowledge about distant habitats, and thus a high degree of randomness in their dispersal routes (McRae & Beier, 2007). The models integrate the contributions of all possible pathways across a landscape or network, rather than a single shortest path (Supporting information Figure S1). As in electrical circuits, the addition of new pathways increases total connectivity but may reduce the flow on any particular path by distributing flow across more routes (McRae & Beier, 2007). Areas of high current-flow centrality correspond to areas of constrained connectivity where few alternate routes exist . Areas of low centrality may either experience no dispersal flow or diffuse dispersal along many alternate routes.
We calculated current-flow centrality using the passage function of the gdistance package in R (van Etten & Hijmans, 2010), with the theta parameter set to 0, the default corresponding to current flow. Although van Etten and Hijmans (2010) refer to the resulting values as "passage density," they are equivalent to the more widely used term "centrality." Calculation of centrality, especially current-flow centrality, is computationally intensive, as the complexity of the problem increases quickly with an increase in the number of nodes or pixels considered. Given n source pixels and m target pixels, complexity is on the order n 2 (n + m) for current-flow methods, vs. nm + n 2 /log(n) for the shortest-path methods described below (Ahuja et al., 1993). A single current-flow centrality raster was calculated for each climate type using all current locations of that type as sources and all future locations of that type as targets ( Figure 1).
Calculation at the extent of North America at 5 km resolution required approximately one week on a 16-core server with 256 GB RAM.
We characterized the results of the current-flow centrality analysis in terms of their distribution in climatic space. We first plotted in climate space (i.e., in terms of PCA1 and PCA2 axes) the sum of the centrality values (for all climate type rasters) found in the spatial locations characterized by each climate type. Each pixel in climate space in these plots represents the connectivity value of the geographic areas which exhibit this climate type. We then plotted the maximum values from the entire raster (i.e., for all geographic areas) of centrality results for each climate type onto a plot of the two PCA axes, both as an unscaled value and as a value divided by the total centrality values for that climate type raster. Each pixel in climate space in these latter plots represents the maximum value across all North America of the centrality map connecting current and future analogs of a single climate type. Climate types with high maximum centrality may have constrained connectivity with their future analogs (i.e., dispersal is limited to a few alternative routes).

| Sensitivity of current-flow results to alternative model structures and extents
We analyzed three alternative current-flow-based models or scenarios to discern the effects of different aspects of the model in producing observed centrality patterns. Firstly, what we term the "full model," which connects present-day climate types to future analogs across a climate-type-specific resistance layer as described above, was used to identify areas of high centrality resulting from constraints imposed by climatic gradients in the intervening landscape.
Secondly, the "uniform-resistance model," which connects present-CARROLL ET AL.
| 5321 day climate types to future analogs across a uniform resistance layer, was used to identify areas of high centrality resulting from the geographic distribution of climate types without consideration of the intervening landscape. The contrast between the uniform resistance and full model is analogous to the contrast between Euclidean (straight-line) distance-based climatic velocity metrics and the minimum-exposure-distance-based velocity metric of Dobrowski and Parks (2016). Finally, the "null model" which connects all points across a uniform resistance layer, was used to identify areas of high centrality resulting from the shape of the North American continent.
We assessed sensitivity of results to the spatial extent of analysis by comparing results from the continental-extent current-flow analysis with those from an analysis at the extent of an example ecoregion, at both 1 and 5 km resolutions. We selected for this comparison an ecoregion (the Central Basin and Range ecoregion of the interior western US; C.E.C., 1997) which encompasses a large temperature and precipitation gradient. We buffered the ecoregion by 200 km to reduce edge effects.

RCPs
To evaluate sensitivity to parameter uncertainty, we also ran the full model using eight representative CMIP5 AOGCMs (CCSM4, CNRM- the 2080s periods. We also ran the full model based on AOGCM ensemble predictions for current to near-future (2050s) projections, and near-future to distant-future (2080s) projections. We compared resultant centrality values using Spearman rank cross-correlation coefficients which provide the mean correlation (i.e., across the entire geographic space) between two rasters expressed as matrices.
Parameter sensitivity analyses were performed at 10 km resolution to increase computational feasibility, a decision supported by our finding that results from the full model at 5 km were highly correlated with those obtained at 10 km (ρ = 0.96).

| Comparison of current-flow centrality results with other climatic and nonclimatic variables
We hypothesized that current-flow centrality results would not be strongly correlated with climatic and topographic variables previously used to identify refugia, but that nonlinear relationships might exist that would be informative in interpreting centrality patterns at the continental extent. We compared centrality output from the full model to seven additional variables using: cross-correlation coefficients and generalized additive models (GAM; Hastie & Tibshirani, F I G U R E 1 Map of shortest-path and current-flow analysis results for an example climate type located in western North America. Forward shortest paths connecting each current climate location to its closest future analog differ from backward shortest paths connecting each future climate location to its closest current analog. Current-flow centrality results connect all current and future locations of the climate type by means of diffuse flow rather than single paths 1990). The seven variables included four nonclimatic variables (latitude, elevation, elevational diversity, and topographic position index).
Elevation data for North America were resampled to 1 km resolution from an original resolution of 1 arc-second (~30 m) (ASTER; ASTER GDEM Validation Team, 2009) to 3 arc-second (~90 m) (SRTM ;Farr et al., 2007). We then derived elevational diversity values as in Carroll et al. (2017), using a form of Rao's quadratic entropy (Rao, 1982), by measuring the mean elevation contrast between all pairs of pixels within a spatial neighborhood defined by the moving window of the specified extent (here 27 by 27 km after Carroll et al. (2017)). Topographic position index (TPI) measures the elevational position of a 5 km pixel in relation to its surrounding neighborhood (here 81 by 81 pixels after Stralberg et al. (2018)) as the pixel elevation divided by the mean elevation of the surrounding neighborhood (Wilson, OConnell, Brown, Guinan, & Grehan, 2007).
The seven variables also included three climatic variables (climatic dissimilarity, forward climatic velocity, and backward climatic velocity). We defined climatic dissimilarity as the Euclidean distance in PCA climatic space between current and future climate at a pixel.
Forward and backward climatic velocities were defined using the analog-based method and represent the straight-line distance between analogs . To evaluate the relationship between connectivity areas and refugia, we also categorized EPA Level III ecoregions (n = 182) into four quadrants based on high and low centrality and backward climatic velocity values.

| Comparison of current-flow and shortest-path centrality
We compared results from the full current-flow-based model with shortest-path centrality values derived using methods adapted from Dobrowski and Parks (2016). Shortest-path models represent cost (e.g., energetic cost or mortality risk) of movement through different habitat (here climate) types as distance (i.e., points in less permeable habitat are conceived as farther apart), in order to identify the route between two predetermined endpoints that has the shortest total distance (least total cost) in terms of the summed resistance of all pixels on the path (Newman, 2010). Shortest-path methods effectively assume that dispersers have complete knowledge about values of distant habitats. Although the dispersal model underlying shortest-path metrics may be less realistic than that underlying current-flow metrics, the greater computational feasibility of shortest-path algorithms allowed us to calculate paths on a pairwise basis between each source location of a climate type in time period A and the target location in the time period B which was "closest" in terms of summed resistance ( Figure 1). This allowed us to assess the sensitivity of centrality values to assumptions regarding maximum dispersal distance thresholds, by comparing results from all path lengths with results which excluded paths of greater than a threshold value (termed boundeddistance centrality; Newman, 2010). We also compared shortest-path centrality results with current-flow centrality and the seven environmental variables described above.
We calculated shortest paths in both forward and backward directions (Figure 1), which are analogous to the forward and backward versions of climatic velocity . Backward climatic velocity represents the distance and rate at which organisms adapted to a location's future climate will need to move to reach that location, and reflects a location's ability to serve as a refugium . In contrast, forward climatic velocity, which represents the rate at which an organism currently at a location must move to find future suitable climate, is more relevant to measuring threats to organisms themselves. Forward velocity will often be high in alpine areas because reaching the nearest analogous future climate may require dispersal to distant higher elevation mountaintops . Conversely, backward velocity is generally low in alpine areas, because adapted organisms can reach the site from nearby downslope locations. Values are often high in valley bottom habitat because organisms must travel longer distances to colonize these locally new habitat conditions. Because pairwise paths between all current and future locations of a climate type are analyzed jointly during the current-flow calculations, the distinction between the forward and backward directions is not relevant to current-flow centrality results ( Figure 1). Therefore, we aggregated forward-and backward shortest path output when comparing shortest-path centrality to current-flow centrality.

| Distribution of high centrality areas in climatic space
Mean annual temperature and precipitation were the variables with the highest loadings on our PCA axes 1 and 2, respectively (Supporting information Table S3, Figure S3). However, several other bioclimatic variables had similarly high contributions, resulting in PCA axes that incorporated more information than would have occurred if only mean annual temperature and precipitation had been used. Areas CARROLL ET AL. | 5323 which had high centrality (as summed over all climate types) were primarily located in the wettest and coldest portions of climate space under both current and future climates (Supporting information Figure S4a,b). The climate types which were most constrained in terms of connectivity between current and future climates (as measured by their maximum centrality values) were more broadly distributed in climate space, and were more common in hot semi-arid climate types than in wetter regions (Supporting information Figure S4c,d).  Table S4). All climatic metrics showed an increasing trend with latitude. Climate dissimilarity and forward velocity showed a stronger trend than did centrality and backward climatic velocity (Supporting information Figure S6a), due to polar amplification of the pace of climate change (Huang et al., 2017) in the case of climate dissimilarity, and to the northern continental edge leading to disappearing climate types in the case of forward velocity. In interpreting the GAM results, note that planners would typically seek to conserve areas with high values for centrality (i.e., potential connectivity areas) but low values for climatic dissimilarity and velocity (i.e., potential refugia).

|
We found a quadratic relationship between current-flow central-  Figure S6d).
The geographic distribution of current-flow centrality values was consistent with GAM results, in that high values were found along major north-south-oriented passes (Figure 2b). Long north-south oriented mountain chains also showed high importance, although centrality was often greatest in midslope areas rather than along summits (Figure 2d). Paths tend to avoid areas of high forward velocity and disappearing climates, especially in lower latitudes, with centrality increasing in the interstices between such areas as paths were constrained to a few potential routes (Figure 2c).

| Robustness of priorities to AOGCM, RCP, time period, and extent
Results of the sensitivity analysis suggested that current-flow centrality values were highly robust to parameter or data uncertainty However, current-flow centrality results for our example ecoregion (Central Basin and Range) were sensitive to analysis extent.
Only moderate rank correlation (0.66) was shown between results at 5 km resolution from the continental-and ecoregional-extent analyses, although general patterns were similar (Figure 3). Ecoregionalextent analyses were highly correlated at resolutions of 5 km vs.

| Contrasts between alternative current-flow models
Results from the null model (Supporting information Figure S7a), which connected all points across a uniform resistance layer, showed areas of high centrality resulting from the shape of the North American continent. Limited path options and hence high centrality values were generated along coastal isthmuses (e.g., the Isthmus of Tehuantepec) and peninsulas (e.g., Aleutian), and between large lakes, as well between Hudson Bay and the Arctic Ocean.
Results from the uniform-resistance model (Supporting information Figure S7b), which connected current climate types to analogous future types across a uniform resistance layer, showed areas of high centrality resulting from the geographic distribution of current climate types relative to their future analogs, but without consideration of climatic heterogeneity of the intervening landscape. In addition to the areas evident in the null model, high centrality was shown more widely across western Canada and southern Alaska, especially along the Pacific coastal ranges. Additional areas were evident in the interior of the continent, especially the western US. Results from the full model (Figure 2), which connected current climate types to analogous future types across a type-specific climatic resistance layer, showed areas of high centrality resulting from constraints imposed by climatic gradients in the intervening landscape. Fine-resolution variation in centrality was added to the patterns evident in the null and uniform-resistance models, especially in topographically complex regions of western North America.
Rank correlation between results of the uniform-resistance model and full model was 0.84. Partial rank correlation between full and uniform-resistance after accounting for the common influence of continental shape (null model) was 0.82, indicating that most commonality between results was attributable to the distribution of climate types, and little to continental shape. GAM relationships between shortest-path centrality and topographic variables were weaker than those shown by current-flow centrality (Supporting information Figure S6, Table S4), but indicated that backward shortest-paths tended to fall at lower elevations and along valleys rather than ridges (Supporting information Figure S6b,c).

| Shortest-path centrality
The minimum path length within a 5 km pixel was relatively low for both forward and backward centrality within montane regions outside of high latitude regions, and relatively high in nonmontane and high-latitude regions (Supporting information Figure S9). Target pixels (destinations of dispersers) for forward shortest-paths were concentrated in coastal and interior montane regions and the high Arctic, whereas target pixels (origins of dispersers) for backward shortest-paths were more widely distributed (Supporting information Figure S10).
Comparison of the mean values by EPA Level III ecoregion between current-flow and shortest-path centrality analyses (Figures   2 and 4) suggested that there was moderate rank correlation (ρ = 0.52, n = 182) between the rank of ecoregions in terms of the two metrics.

| Centrality, protected areas, and the human footprint
Although mean current-flow centrality values were slightly higher within protected areas than outside of such areas (6.18 vs. 5.17, respectively), protected areas, which covered 11.4% of the continent, held only 13.3% of total current-flow centrality value. Shortest- and southeastern California (Supporting information Figure S11).

| DISCUSSION
The persistence of many species under climate change will hinge on their ability to disperse and colonize habitat which has become newly suitable for their climatic requirements. Many conservation plans focus on promoting climate adaptation by identifying climatic refugia (Keppel & Wardell-Johnson, 2012), but few have considered which areas merit protection for their role in facilitating dispersal under climate change (Littlefield et al. (2017), McGuire et al., 2016).
This study, the first to characterize the location of climate connectivity areas at a continental extent, can help planners identify and protect such areas from land use that may impede movement (Heller & Zavaleta, 2009).
Our results demonstrate that the priority areas for conservation of climatic connectivity are distinct from refugia, and thus poorly captured by many existing conservation strategies. Many climate connectivity areas lie outside of protected areas and face pressure from anthropogenic land use and habitat conversion. Our results suggest that connectivity priorities are robust to uncertainty as to the magnitude of future climate change arising from differing emissions pathways and general circulation models, but vary with analysis extent and differing assumptions concerning dispersal behavior and distance.

| Connectivity areas differ from refugia
Climate connectivity areas (areas which showed high current-flow centrality because large numbers of dispersal paths overlapped) did not coincide with areas of low climatic velocity (i.e., potential refugia; Keppel & Wardell-Johnson, 2012;Carroll et al. 2017), supporting the relevance of grouping sites into categories of conservation targets based on high and low centrality and refugia values (Supporting information Figure S5). The underlying cause of this contrast is that, as terrestrial organisms disperse in response to shifting climates, they must often traverse intervening nonrefugial areas to reach their targets, and are constrained by broad-scale topographic features that do not affect velocity metrics based on straight-line distances. While velocity metrics are useful in capturing factors determining macrorefugia location, centrality captures the broad-scale factors that determine how such macrorefugia are connected by dispersal.
Like velocity, centrality integrates factors operating at a range of spatial scales . For example, topographically diverse areas may facilitate climate connectivity over short distances, but broad-scale topographic features such as mountain passes may be key areas for facilitating longer distance dispersal.
Both the generalized additive model (GAM) results and mapped output suggested that high centrality areas were associated with north-south trending passes and valley systems, as well as northsouth trending ridgelines at higher latitudes. This latitudinal factor may be due to the fact that, in a warming climate, tropical montane  Figure S4). Current-flow centrality values were also high in areas of the far north where polar amplification of climate change increased the distance between a current climate and its future analogs (Huang et al. 2017), and the shape of the continent funneled dispersal paths into limited areas.
The current protected area system does not capture climatic connectivity areas much above what would be expected by random placement. However, existing protected areas in the southwestern US and north Pacific coast of British Columbia and Alaska include areas of high centrality. Expansion of the human footprint associated with exurban growth in these regions (e.g., in southeastern California and British Columbia) could result in land use changes that reduce the permeability of the landscape to native species (Batllori, Parisien, Parks, Moritz, & Miller, 2017). Similarly, the high climate connectivity value of the western slope of the Appalachians is potentially impacted by intensification of human land use. The high importance of the Arctic for climate connectivity may also be threatened by future expansion of the human footprint in that region.
Some approaches to climate adaptation planning avoid use of future climate projections due to the inherent uncertainty in forecasting future emissions and their effect on global climate (Beier, Hunter, & Anderson, 2015). However, we found that climate connectivity priorities were robust to uncertainty caused by variation in input data between alternative emission pathways (RCP 4.5 and 8.5) and AOGCMs. Priorities based on centrality derived from different AOGCMs showed less between-AOGCM variation than did the input climate data or climatic velocity metrics based on that climate data . Conservation planners may place greater focus on sensitivity analyses of the effect of assumptions concerning dispersal behavior and analysis extent, rather than on parameter uncertainty arising from differing climate projections.
The high correlation between centrality results based on near-future and distant-future climate projections should offer planners more confidence that climate corridors protected based on shorter time horizons will remain valuable over longer time frames, assuming a linear trend in changes in climate. However, identification of continuous paths tracking climate over multiple intermediate time steps, as is possible using network flow approaches, would be a useful refinement of our approach that avoids such assumptions (Graham et al., 2010;Phillips et al., 2008).

| Relevance of graph-based models to planning for dispersal under climate change
Although highly simplified representations of dispersal, graph-based models are often more computationally feasible than more realistic dispersal simulations, allowing analyses at finer resolution or broader spatial extent (Adriaensen et al., 2003). While our approach is applicable at spatial extents ranging from continents to local regions, continental-extent analyses are an important initial step in discerning patterns that would not be evident in an analysis limited to a smaller extent, for example, an ecoregion (Rudnick et al., 2012). For example, because future climate analogs for some climate types will only occur outside an ecoregion, an analysis limited to the extent of a single ecoregion would not be able to evaluate connectivity for these climate types. Results from an analysis limited to the extent of a single ecoregion were only moderately correlated with values for the same region derived from the continental-extent analysis, suggesting that planners should compare metrics derived at broader and more limited extents to assess robustness of conclusions to spatial scale.
Because the two types of graph-based models (shortest-path and current-flow-based centrality) applied here make contrasting assumptions about the ability of dispersers to find optimal paths, comparison of their output provides complementary information about priority areas and helps assess robustness to assumptions about dispersal behavior (Carroll et al., 2012;Carroll, Fredrickson, & Lacy, 2014;García Molinos et al., 2017). Shortest-path results can help identify a minimal "skeleton" of landscape linkages which may not be evident from current-flow results. Conversely, current-flow approaches identify both specific pinchpoints along this minimal network and broader regions with high importance for connectivity (Carroll et al., (2014); McRae et al., (2016)). Current flow results are also more biologically realistic in that they are more influenced by flow between nearby points than between distant pairs of points ( Figure 1). However, the greater computational feasibility of shortest-path algorithms allows more detailed analyses comparing forward and backward climate connectivity areas, and results based on differing maximum dispersal distances. The relatively coarse resolution of our continental-extent analysis, however, limits our ability to identify fine-scale refugia important to species with limited dispersal capabilities (Stralberg et al., 2018).  (Dobrowski & Parks, 2016) and climate centrality become relevant in assessing conservation priorities.
However, several factors limit the comparability of the shortestpath lengths from our models with dispersal distances recorded from field observations. To the extent that published summaries of dispersal rates (Settele et al., 2015) are based on measurements of straight-line displacement, they may be more comparable to standard climatic velocity values than to circuitous climatic-resistance-based paths, which averaged 45% longer than straight-line displacements in our results (Supporting information Figure S8). Both climatic-resistance-based and straight-line velocity values are sensitive to how finely continuous climatic values are divided into types (Dobrowski & Parks, 2016;Hamann et al., 2015). The climate type width used here allow assessment of climatic tolerances narrower than those shown by most species, and can help inform conservation of locallyadapted populations, an often-ignored facet of the extinction crisis (Ceballos & Ehrlich, 2002). Graph-based assessments of climate connectivity derive their practical relevance not from direct comparability with field data on dispersal but because they provide a relative measure of climate exposure that is more realistic for many organisms than is straight-line climatic velocity, and allow planners to identify and protect landscape features of importance for broad-scale connectivity. Climate connectivity areas represent important conservation targets even in regions of high climatic velocity (red areas in Supporting information Figure S9) where unassisted migration will be feasible for a smaller subset of the biota.

| Coarse-and fine-filter approaches to assessing connectivity
Our approach resembles many previous studies in defining source and target areas based on climate data (i.e., climate analogs), effectively using climate types as a "coarse-filter" surrogate for biodiversity targets (Noss & Cooperrider, 1994). Alternatively, several studies have instead identified species-specific source and target areas using future species distributions based on climatic niche models. Schloss, Nunez, and Lawler (2012) used niche model projections to measure straight-line distance and velocity between current and future species climatic niches across the Western Hemisphere. Lawler, Ruesch, Olden, and McRae (2013) extended this niche-model-based approach by mapping current-flow-based paths between current and projected futures species niches. More recent work (Miller & McGill, 2018) has simulated dispersal of temperate tree species under climate change using detailed data on species-specific life histories.
Coarse-and fine-filter approaches are often seen as complementary (Tingley, Darling, & Wilcove, 2014), because coarse-filter analyses provide context by allowing a comprehensive depiction of patterns across climate space rather than being limited to the climatic niches of the subset of taxa for which distribution data are available (Brito-Morales et al., 2018). For example, coarse-filter analyses, especially at the resolution permitted by graph-based methods, allows identification and characterization of topographic features associated with climate connectivity areas in a manner that would be difficult using species-specific models.

| 5329
The robustness of our results to uncertainty from alternative AOGCMs and emissions pathways should give planners increased confidence in applying conservation priorities based on future climate projections. The sensitivity of results to model uncertainty should be expected given that the different approaches used here make contrasting assumptions concerning dispersal behavior. Planners should see these contrasting results as complementary information sources that must be viewed in context of the strengths and limitations of graph-based models. Although the metrics used here are necessarily simplified representations of real-world dispersal processes, they can add a missing dimension to climate adaptation planning by identifying landscape features which promote connectivity among refugia, a key management strategy under climate change (Heller & Zavaleta, 2009).

ACKNOWLEDGMENTS
The authors dedicate this paper to the late Brad McRae, in appreciation of his contributions to advancing both this work and the broader field of connectivity conservation science. The manuscript benefited from suggestions by P. Elsen and two anonymous reviewers. Co-authors are listed in the order of relative contribution (SDC system). The Wilburforce Foundation provided support for the authors as part of the AdaptWest Climate Adaptation Planning Project (https://adaptwest.databasin.org; see site for data associated with this paper).