Climate change likely to reshape vegetation in North America's largest protected areas

Climate change poses a serious threat to biodiversity and unprecedented challenges to the preservation and protection of natural landscapes. We evaluated how climate change might affect vegetation in 22 of the largest and most iconic protected area (PA) complexes across North America. We use a climate analog model to estimate how dominant vegetation types might shift under mid‐ (2041–2070) and late‐century (2071–2100) climate according to the RCP 8.5 scenario. Maps depicting vegetation for each PA and time period are provided. Our analysis suggests that half (11 of 22) of the PAs may have substantially different vegetation by late‐21st century compared with reference period conditions. The overall trend is toward vegetation associated with warmer or drier climates (or both), with near complete losses of alpine communities at the highest elevations and high latitudes. At low elevation and latitudes, vegetation communities associated with novel climate conditions may assemble in PAs. These potential shifts, contractions, and expansions in vegetation portray the possible trends across landscapes that are of great concern for conservation, as such changes imply cascading ecological responses for associated flora and fauna. Overall, our findings highlight the challenges managers may face to maintain and preserve biodiversity in key PAs across North America.


| INTRODUCTION
Protected areas (PAs), including national parks (NPs), wilderness areas, and nature reserves, are legally designated and managed to conserve natural ecosystems or to protect landscapes comprising unique biological, geological, physical, or ecological features (Dudley, 2008). About 15.4% of the world's terrestrial areas are sheltered within PAs under governmental jurisdiction (Pringle, 2017). Yet, this network of PAs, and their associated flora and fauna, face unprecedented threats from climate change and other human pressures (e.g., land-use conversions; IPCC, 2014). Even the most remote PAs and those with the highest level of protections (e.g., Strict Nature Reserves and Wilderness Areas; Dudley, 2008) are not safeguarded from climate change (Batllori, Parisien, Parks, Moritz, & Miller, 2017;Gonzalez, Wang, Notaro, Vimont, & Williams, 2018).
The creation of PAs commenced over a 100 years ago in North America with the intent to maintain a sample of pristine wildlands; managers endeavored to preserve or restore landscapes to conditions prevailing prior to the arrival of Europeans (Leopold, King, Cottam, Gabrielson, & Kimball, 1963). Today, managers must consider a far more dynamic future (Colwell, Avery, Berger, & Davis, 2014) as the Earth warms at a rate unmatched in the Holocene (Hansen et al., 2006). With continued emissions of greenhouse gases, climatic isotherms are projected to shift-often poleward and upslope-with a 3 C change in mean annual temperature in temperate zones. This corresponds to an isothermal shift of 300-400 km latitude or 500 m in elevation (Hughes, 2000). Such climatic shifts are projected to trigger major changes in the distribution of terrestrial ecosystems (Davis & Shaw, 2001;Eigenbrod, Gonzalez, Dash, & Steyl, 2015;Gonzalez, Neilson, Lenihan, & Drapek, 2010) and alter disturbance regimes (e.g., wildland fire, insect outbreaks), droughtinduced vegetation mortality, and their interactions (Allen, Breshears, & McDowell, 2015;Bentz et al., 2010). A key challenge for PA managers is deciding how to respond to the accelerating change. Possible strategies include stewardship measures to promote adaptation to a new climatic environment (e.g., increasing connectivity among PAs to enhance species' movement) or measures to maintain existing ecosystems (e.g., decreasing forest density to increase drought resistance) (c.f., Millar, Stephenson, & Stephens, 2007). In either case, managers would benefit from a greater understanding of ecosystems' potential responses to climate change.
Predicting future ecological conditions, especially the distribution of dominant vegetation, is typically achieved using either correlative or process-based model approaches.
Correlative approaches, such as species distribution modeling (SDM), relate the presence of species with environmental covariates to predict distributions (Elith & Leathwick, 2009). Process-based approaches, such as dynamic global vegetation models (DGVMs), more explicitly incorporate physiological and population-level processes (e.g., growth, regeneration, mortality) and their interactions (Hartig et al., 2012;Hickler et al., 2012). A third approach, the climate analog model, has recently been introduced as an alternative to SDMs and DGVMs (Parks, Holsinger, Miller, & Parisien, 2018;Pugh et al., 2016). The climate analog model was developed based on the literature pertaining to climate velocity, in which the locations with the best climatic match between one time period (e.g., future) and a different time period (e.g., contemporary) are considered climate analogs (Batllori et al., 2017;Hamann, Roberts, Barger, Carroll, & Nielsen, 2015;Wuebbles & Hayhoe, 2004). The climate analog model assumes that future ecosystem characteristics (e.g., dominant vegetation or crop yield) at a given site can be characterized by the current ecosystem characteristics of other sites that are identified as climate analogs. Because climate variables, rather than species data, are used to describe potential ecosystem shifts, the climate analog model avoids many of the species-level assumptions inherent to SDMs or the "data-hungry" parameterization required by DGVMs. It also offers an advantage over correlation-based approaches, which fit the central tendency of observations and in so doing, forfeit the description of much of an ecosystem's variability. The analog-based approach lacks this constraint and more directly integrates the intrinsic spatial variability of climate-vegetation interactions and to a certain extent, ecological feedback from environmental drivers . The climate analog model has been used to infer potential changes in crop yields and vegetation patterns under changing climate (Parks, Dobrowski, Shaw, & Miller, 2019;Pugh et al., 2016).
This study uses the climate analog model to examine how vegetation in some of the largest PAs across North America (México, United States, Canada; n = 22) may change as the climate warms. Specifically, we used a bivariate combination of climate variables known to influence vegetation distribution (Stephenson, 1990) and gridded data sets representing vegetation type to evaluate the potential for vegetation change from contemporary conditions into the mid-and late-21st century. Batllori et al. (2017) conducted a similar continent-wide assessment but aggregated results at the continental scale (i.e., land cover changes within PAs were not reported). Here, we advance this work by narrowing the focus to an evaluation of potential vegetation shifts within individual PAs. We present fine-scale spatial projections (i.e., maps) for a subset of PAs across North America and quantify the potential for existing vegetation types to shift under a changing climate for each PA.

| Methods
We inferred climate-induced shifts in natural vegetation communities within North American PAs at a 1-km spatial resolution. We used gridded vegetation maps for each country to define reference conditions for each PA, and estimated the potential for change for two future time periods representing mid-and late-21st century. Vegetation changes were evaluated using the climate analog model, where, for any given pixel in a PA, we identified the vegetation in the nearest locales with current climate conditions that match its projected future climate (c.f., Batllori et al., 2017;Parks et al., 2018). This approach assumes that, at large spatial scales, climate is the dominant factor controlling current and future vegetation distributions, and ecosystem interactions with natural disturbances are implicitly represented (Prentice et al., 1992;Stephenson, 1990).
Protected areas. We chose 22 study sites (encompassing nearly 500,000 km 2 ) composed of individual or contiguous PAs across North America based on the World Database of  Figure 1) were selected based on their size (>10,000 km 2 ) and international recognition for conservation and social interests (i.e., IUCN category I or II). However, IUCN status and size criteria were relaxed in eastern Canada, eastern United States, and México to provide broader geographical representation. Gridded vegetation data. We used gridded vegetation maps specific to each country to represent the best data in terms of classification accuracy, spatial resolution, and in equilibrium with climate (see Supporting Information I). The gridded vegetation data set for Canada was derived from land-use/land-cover grids based on Moderate Resolution Imaging Spectroradiometer (MODIS) (Pouliot, Latifovic, Zabcic, Guindon, & Olthof, 2014) and Advanced Very High Resolution Radiometer (AVHRR) (Latifovic & Pouliot, 2009) imagery; for the United States, it was derived from LANDFIRE Biophysical Settings (Rollins, 2009); and for México it was derived from a land-use/land-cover grid (INEGI, 2017). Because the number of vegetation classes varied widely by country (12-523 per country), we aggregated these classes into broader groupings (8-24 per country) to make projections more robust and to facilitate comparisons (Supporting Information I). We chose not to use available continental-wide vegetation maps, such as the North American Land Change Monitoring System (NALCM) data set (USGS Land Cover Institute, 2010), based on our observations of classification inaccuracies at fine spatial scales (in particular, those associated with natural disturbances). For example,~50% of Yellowstone NP was classified as grass or shrubland in the NALCM data set rather than the early seral lodgepole pine that re-colonized those areas after the 1988 fires (Romme et al., 2011). Consequently, we chose maps individual to each country that best demonstrated vegetation patterns driven by long-term historical climate rather than short-term events (i.e., fire, insect outbreaks).
Climate. We obtained gridded climate data (1-km resolution) from Wang, Hamann, Spittlehouse, and Carroll (2016) characterizing climatic normals (annual data summarized into three-decade averages) for three time periods: 1961-1990 (reference), 2041-2070 (2055; mid-century), and 2071-2100 (2085; late-century). We chose Hargreave's climatic moisture deficit (CMD, mm/year) and Hargreave's reference evaporation minus CMD (hereafter evapotranspiration [ET, mm/year]). CMD and ET characterize site-level evaporative demand and water availability (Stephenson, 1998) and are recognized as strong predictors of vegetation across North America (Dyer, 2002;Stephenson, 1990). Our primary analysis uses the MPI-ESM-LR global climate model under the RCP 8.5 radiative forcing scenario (Moss et al., 2010), which is representative of "median" climate change for North America and has high validation statistics (Knutti, Masson, & Gettelman, 2013; projected changes in CMD and ET shown in Supporting Information II, Figure S1, S2). Recognizing variability in Global Climate Model (GCM) projections, we also evaluated two additional GCMs which characterize a future with less change in temperature and precipitation (INM-CM4) and higher increases in temperature and precipitation (GFDL-CM3) (see Supporting Information II, Figure S3), and we provide the full complement of results in Supporting Information III.
Projected vegetation change. We assumed that existing vegetation in a given location may also persist in other locations with analogous climatic conditions across space and time (Ohlemüller, Huntley, Normand, & Svenning, 2012). Following the rationale and general methods in Parks et al. (2018), we first identified analogs by evaluating where unique climate combinations under future time periods (midand late-21st century) were located under the reference time period . We assumed that the future vegetation at a site is represented by the reference period vegetation at the locations identified as climate analogs. Our approach to inferring future vegetation in the 22 PAs started with a characterization of climate (i.e., CMD and ET) for each pixel in North America for all three time periods (reference, mid-and late-century). Then, for each pixel in a PA, we identified the nearest (in geographic distance) seven pixels whose reference period climate was within ±1 mm/year CMD and ET (after the square-root transformation) of the future time period climate of the pixel of interest. These seven pixels represent climate analogs, and we used the majority vegetation type among them to represent the pixel of interest's future vegetation. Where ties occurred, we incrementally increased the number of considered analogs until the tie was broken. If the focal pixel is an analog with itself, it was counted as one of the seven nearest analogs. Climate bin width influences the precision of analog detections (Carroll, Lawler, Roberts, & Hamann, 2015); the bin widths for this study were selected based on the highest accuracy for classifying forest versus non-forest groups across the U.S. West Supporting Information II).
We evaluated upland vegetation only, excluding hydrologically controlled land cover types (i.e., wetlands) from the analysis. Although freshwater wetlands can convert to upland vegetation (and conceivably vice versa), incorporating this complex dynamic was beyond this study's scope. We also excluded water bodies and a small proportion of F I G U R E 1 Protected areas for which we projected vegetation shifts in North America pixels classified as urban or agriculture within the Mexican PAs. Thus, pixels classified as water, wetland, urban, or agriculture were not evaluated in the pool of analogs; and no change was projected for these pixels.
Differences in vegetation classification schemes created some challenges when analogs for a pixel occurred in a different country. For example, the climate analogs for some Canadian PA pixels were located in the United States (and vice versa, but to a lesser degree). We took several approaches to improve the correspondence of climate analogs to the classification framework for each country. In Canada and Alaska, it was fairly straightforward to crosswalk between the classification schemes of each country. Glacier and North Cascades complexes, however, were less straightforward. Future climate analogs for these U.S. PAs could be located in Canada, and because the Canadian vegetation classes were broader, there was no clear crosswalk of these coarsely defined classes to the finer scaled U.S. vegetation classes (e.g., evergreen needleleaf not easily assigned to either cold montane conifer or western mesic conifer). In these limited cases, we constrained the search for climate analogs to the conterminous United States. In total, six PAs had cross-border incongruences but only two were affected to a large extent: in year 2085, Algonquin Provincial Park had~90% of its vegetation as non-resident (i.e., originating in the United States) and Wabakimi Provincial Park had~20%, whereas the remaining PAs were nominally affected (<3%; see Supporting Information I for methods to resolve cross-border incongruences).
Our primary objective was to produce maps depicting possible changes within PAs; however, we also quantified potential change in vegetation composition across PAs based on a variant of standard Euclidean distance which is sensitive to relative proportions among vegetation types (Ludwig & Reynolds, 1988), as follows: where RED is the relative Euclidean distance, A i is the area of vegetation type i during the reference period, and B i is the area of vegetation type i in the future time period, for the n vegetation types occurring in each PA. RED ranges from 0 to √2; we normalized RED and multiplied by 100 (to represent units in percent), such that 0 indicates the distribution of vegetation does not change and 100 signifies complete dissimilarity.
We also quantify the percent vegetation composition for each time period, and plot vegetation transitions from reference conditions to year 2085 as chord diagrams. In addition, we evaluate: (a) the degree of agreement among analogs meeting climate bin width criteria (i.e., ±1 mm/year CMD and ET) to the final vegetation type assignment (by majority) for each focal pixel, (b) central tendency statistics for the Euclidean distance of those analog pixels to focal pixel; and (c) correspondence between agreement among analog pixels and Euclidean distance from analog to focal pixel.

| Results
For our primary analysis using the median climate scenario (MPI-ESM-LR model and RCP 8.5 emissions), we found the potential for substantial changes in vegetation across most PAs by mid-to late-21st century ( Figure 2). By mid-F I G U R E 2 Change in vegetation composition from the reference condition, based on the relative Euclidean distance index, across all 22 protected areas and both time frames. Colors differentiate Alaska, Canada, contiguous United States, and México century, over a third of the 22 PAs could experience vegetation composition changes of more than 25% and up to 55%. By late-21st century, 11 of 22 PAs could undergo changes greater than 25% and up to 62%, with the southernmost PAs in Canada and the northernmost PAs in the conterminous United States showing the largest potentials for vegetation change as estimated by RED (Figure 2).
Our analysis using the median climate scenario projected a variety of changes to the spatial distribution and composition of vegetation within each of the PAs. By late-21st century, the overall composition of Denali NP is projected to be relatively stable, but this park may see an increase in boreal mixed forests to the west and north of the highelevation alpine-covered Alaska Range; this increase is primarily at the expense of boreal conifer forest and arctic tundra (Figure 3). Wabakimi Provincial Park in Canada may also see an increase of mixed forest, but primarily at the expense of evergreen needleleaf forest (Figure 4). By late-21st century, Wabakimi Provincial Park could also experience the emergence of deciduous broadleaf forest throughout its north-central region, a vegetation cover not currently mapped in the park. The Sierra Nevada complex in the western United States likewise could experience shifts in forest type with transitions away from cold montane conifer to other western forest types; these changes are broadly distributed throughout the complex (Figure 5). At the same time, the overall amount of cold montane conifer forest in the Sierra Nevada complex stays fairly stable due to projected transitions from alpine cover. Tehuacán-Cuicatalán in México may experience gains in tropical deciduous lowland dry forest at the expense of near total losses in desert scrub and oak and pine-oak forest vegetation by late-21st century ( Figure 6). Tehuacán-Cuicatalán may also experience a shift to vegetation associated with a novel climate along its eastern boundary. The other PAs we analyzed provide similar examples of potential shifts in the distribution and composition of vegetation. In several of the farthest north PAs (Gates of the Arctic, Wrangell-Kluane, Tuktut Nogait, Nahanni, Banff-Jasper), anticipated trends include decreasing alpine cover and increasing forest and/or tundra (Supporting Information III, Figure S1, S3, S4, S6, S7). Trends toward deciduous, broadleaf, or oak forest types are widely anticipated (Supporting Information III, Figure S8-S11, S16, S18, S19), as are shifts among conifer forest types (Supporting Information III, Figure S12-S15). We also found numerous examples of PAs that are projected to gain new or novel vegetation types by late-21st century (Supporting Information III, Figure S2, S4, S17-S22).
We found variability in these vegetation projections when evaluating potential trends under climate conditions with greater temperature and precipitation changes, in particular by the GFDL-CM3 GCM under the late-21st century time frame. Some PAs may experience a much greater extent of vegetation turnover with this more extreme climate projection, such as the Glacier complex where cold montane conifer forest may transition largely to western mesic conifer forest (Supporting Information III, Figure S13) or the Sierra Nevada complex where cold montane forest may be replaced by dry conifer forest, deciduous forest, and shrubland (Supporting Information III, Figure S16). In other PAs, vegetation types not currently present within PAs may become widespread or vegetation from no-analog climates may form. In the northern latitudes, Denali NP may see western mesic conifer forest and shrubland spread into the park (Supporting Information III, Figure S2). Tuktut-Nogait NP and Thelon Wildlife Sanctuary may have evergreen needleleaf extend into their landscapes (Supporting Information III, Figures S4, S5). The Great Smoky Mountains NP may transition from mainly deciduous forest to a mixture of temperate conifer forest, sparsely vegetated/barren, and vegetation from a novel climate (Supporting Information III, Figure S18). PAs may also undergo fundamental ecosystem changes such as transitioning from forest to non-forest landscapes. In Canada, the Algonquin Provincial Park may switch from mixed forest to shrubland (Supporting Information III, Figures S9, S10). In México, Cumbres de Monterrey may lose its pine/oak forest to become dominated by scrub, shrub and grassland (Supporting Information III, Figure S20). The level of agreement between possible analogs (i.e., closest pixels meeting climate bin width criteria per focal pixel) and the final vegetation type assignment for a pixel ranged from median values of~60 (Wood Buffalo complex) to 100% (North Cascades complex) across the 22 PAs (Supporting Information III, Figure S23). The median Euclidean distance of analog pixels to focal pixel ranged from 7 (Death Valley NP) to 293 km (Algonquin Provincial Park) with the greatest distances generally occurring in Canadian PAs (Supporting Information III, Figure S24). We found no clear relationship between the level of agreement among analog pixels and Euclidean distance from analog to focal pixels (Supporting Information III, Figure S25).

| Discussion
Globally, climate isotherms are projected to redistribute during the 21st century, with serious implications for range shifts in species' distributions (Kelly & Goulden, 2008). Within 22 of the largest and most iconic PAs in North America, vegetation communities will likely respond to this climatic redistribution since organisms will be pushed beyond their physiological tolerances (Woodward, 1987). Based on our analysis, the majority of the PAs we analyzed may experience substantial departures (25-62%) from current vegetation communities. Directional trends toward vegetation associated with warmer or drier climate conditions, or both, could not only shift vegetation patterns but render some PAs climatically unsuitable for key species and habitats (Araújo, Cabezas, Thuiller, Hannah, & Williams, 2004).

| VEGETATION CHANGES
Projected changes vary geographically as expected, with an overall trend toward vegetation associated with warmer or drier climates (or both). Most PAs in North America, including the 22 we examined, are situated in mountainous highelevation terrain (Michalak, Lawler, Roberts, & Carroll, 2018) that currently supports extensive alpine communities. Where it currently exists, alpine cover may be replaced by coniferous forest, given suitable soil conditions, which is consistent with expected upward tree-line shifts in mountainous environments (Harsch, Hulme, McGlone, & Duncan, 2009). Another possible trend is that forest types currently associated with cold climates may transition towards more temperate conifer forests. In turn, dry coniferous or mixed forests across Canada and the western United States may trend toward deciduous forest. In the southernmost PAs, we expect mesic deciduous and pine-dominated systems to F I G U R E 5 Potential changes in vegetation distribution in the Sierra Nevada complex: (a) maps of vegetation under reference conditions, and years 2055 and 2085, (b) percent vegetation composition for each time period, excluding agriculture, urban, water, wetlands, and any vegetation types that occupied <1% of the protected area, (c) projected shifts among vegetation types from reference conditions to year 2085, where the sector size in the outer circle indicates the distribution of change over vegetation types, and the thickness of connecting flows depict the amount of transition from one vegetation type to another (vegetation names abbreviated for simplicity; any vegetation types representing 5% or less of the landscape were excluded) transition towards oak-dominated or dry deciduous forests. In México, low elevation desert scrub, mesquite forest and chaparral may shift upward in elevation, and into mountainous PAs. Within individual PAs, exceptions and variations to these patterns exist, but overall patterns are consistent with observations and predictions of others (Kelly & Goulden, 2008;McIntyre et al., 2015).
The extent and rate at which such potential shifts in vegetation actually unfold will highly depend upon the nature of climate change. We examined potential changes under the emission scenario RCP8.5, which assumes no climate mitigation target; yet even within this one trajectory, a range of possible futures was evident, as represented by various GCMs. Under a more extreme climate variant (i.e., GFDL-CM3 GCM), more widespread and rapid changes of in situ vegetation classes were apparent, as was amplified incursion of new vegetation types or unknown assemblages from a novel climate (see Supplemental Information III). Clearly, uncertainties exist in any projection based on climate modeling. However, the potential for broad-scale changes across landscapes evidenced here suggests that climate conditions could become dissimilar enough to push a PA's ecosystem beyond an ecological threshold and towards an entirely new state (Gunderson, 2000). The potential for abrupt threshold responses can alert managers to the array of scenarios to consider when developing planning approaches for the future (Stephenson & Millar, 2012).

| DISAPPEARING EMBLEMATIC VEGETATION AND NEW SPECIES ASSEMBLAGES
Significant among our projections is the potential for many North American PAs to lose the community features that define them (Barrows & Murphy-Mariscal, 2012). One example is the biosphere reserve Tehuacán-Cuicatlan, established to preserve the endemic cacti of this region (Sáenz-Romero et al., 2010); our results suggest that the habitat for these succulent species may decline by over 80% by 2100 in this PA. A second example is Glacier NP, home to 150 glaciers when established in 1910; in this study (where glaciers were classified as snow, ice, barren, and so on), all are projected to disappear by century's end (Hall & Fagre, 2003). Third, old-growth eastern white pine (Pinus strobus) stands are protected in Algonquin Provincial Park, but potential conversion from coniferous to deciduous forest, as projected here, may herald the decline and disappearance of such remnant stands (Joyce & Rehfeldt, 2013). All these potential losses and extirpations highlight the vulnerability of species to become reduced to isolated pockets or pushed out of areas established to protect them.
Several PAs may develop climates not found on the continent today, and entirely novel assemblages of vegetation may form . As the planet warms, we anticipate the warmest PAs in North America, located at lower elevation or latitudes, to become even hotter and drier. For example, Death Valley NP, a PA widely recognized for its high endemism (Baldwin et al. 2017), may develop climatic conditions not observed in North America today. Sierra Gorda and Tehuacán-Cuicatalán may also develop novel vegetation due to hotter, drier conditions with no analogs within North America. These results are consistent with a general expectation for novel climates at the warmer, southern margins of the continent and at lower topographic positions (Mahony, Cannon, Wang, & Aitken, 2017;Williams, Jackson, & Kutzbach, 2007).

| CASCADING ECOLOGICAL EFFECTS
The projected vegetation shifts imply a cascade of ecological changes. For example, we anticipate tundra vegetation to expand in the Gates of the Arctic complex, a trend already observed in this region (Myers-Smith et al., 2011). If this transition is toward tundra that is tall shrub-dominated, highquality forage habitat for caribou (Rangifer tarandus granti) could be greatly reduced (Fauchald, Park, Tømmervik, Myneni, & Hausner, 2017) and, conversely, habitat for moose (Alces alces gigas) could be enhanced (Tape, Sturm, & Racine, 2006). In high-elevation mountainous PAs in western North America, loss of alpine habitat would affect plant species adapted to a short growing season and extreme cold and wind (tussock grasses, dwarf trees, small-leafed shrubs, and heaths) and the animals dependent on this vegetation, such as pikas, marmots, mountain goats, bighorn sheep, and ptarmigan (Saunders, Easley, Logan, & Spencer, 2006). Already, the American pika (Ochotona princeps) has experienced climate-mediated extirpations and upslope range contraction, especially in the Great Basin  but also in Yosemite (Moritz et al., 2008) and Glacier NPs (Moyer-Horner, Beever, Johnson, Biel, & Belt, 2016). In Yellowstone NP, a declining snowpack has increased browsing impacts of elk (Cervus elaphus) on quaking aspen (Populus tremuloides), presaging greater potential for aspen decline (Brodie, Post, Watson, & Berger, 2011). In the Rocky Mountains of Canada and United States, loss of cold montane conifer forest may further accelerate the decline in whitebark pine (Pinus albicaulis) populations in multiple PAs (Banff-Jasper, Glacier, Yellowstone, and Idaho Wilderness), affecting bird and mammal species, such as Clark's nutcracker (Nucifraga columbiana) and grizzly (Ursus arctos) and black bears (Ursus americanus) dependent on this at-risk, keystone species (Keane, Holsinger, Mahalovich, & Tomback, 2017). It would be impossible to address the full array of potential changes to individual species and their interactions, but these examples illustrate that we can anticipate cascading effects on flora and fauna dependent upon existing vegetation communities.

| COMPARISON TO OTHER MODELS
Our projections generally parallel those from other models for future vegetation in North America, although direct comparisons are challenging due to disparate methods (i.e., DGVM, correlation-based), GCMs, emissions scenarios, vegetation frameworks, and spatial scales (coarser, or often of regional extent). In Alaska, a trend toward greater woody cover has been projected (Pearson et al., 2013), similar to our expectations for the Gates of the Arctic complex and Denali NP. Across Canadian PAs, temperate forests and woodland are expected to increase at the expense of tundra (Rowland, Fresco, Reid, & Cooke, 2016;Scott, Malcolm, & Lemieux, 2002), a trend we also found for the country's northern parks including the Wrangle-Kluane complex. In the Canadian Rocky Mountains, Shafer, Bartlein, Gray, and Pelltier (2015) projected shifts away from alpine and cold forest and toward cool forest and open cool forest/woodland, similar to our projections for the Banff-Jasper complex. In the western Pacific Northwest, Sheehan, Bachelet, and Ferschweiler (2015) projected a shift from predominantly conifer to warmer mixed forests, fairly consistent with our findings for dry conifer forests infiltrating the North Cascades complex. For Glacier and Yellowstone NP, Hansen and Phillips (2015) projected a decline in alpine tundra and subalpine/montane conifer forest, with increases in mesic conifer forest, deciduous forest, and scrubland-as does our study. In the eastern United States, Rehfeldt, Crookston, Sáenz-Romero, and Campbell (2012) projected relative stability for evergreen-deciduous forests, generally corresponding with our results for oak-forests to endure across the Great Smoky Mountains NP-at least in the near term. In México, projections include an expansion of tropical dry deciduous forest (Gómez-Mendoza & Arriaga, 2007;Rehfeldt et al., 2012), a trend consistent with our projections especially for Sierra Gorda and Tehuacán-Cuicatalán. Most of these models, including our own, underscore the potential for substantial shifts in vegetation patterns in North American PAs over the next century.
While results from our analysis share commonalities with those using other modeling methods, the climate analog approach offers advantages over others. For one, it carries no assumptions about the climatic tolerances of species (Veloz et al., 2012) and may better capture a broader potential for species to migrate rather than relying on evidence from realized niches in current distributions . Also, the climate analog model may better represent the potential for range shifts because it does not depend on curve fitting which may mute important spatial variability. Finally, the climate analog model is relatively simple, efficient, and flexible to implement, and its intuitive approach can be readily understood by public and policy makers.

| STUDY LIMITATIONS
The climate analog model produces results that contain an inherent degree of uncertainty which is essentially immeasurable and challenging to characterize. For one, modeling future climate as an average over multiple years based on one representative concentration pathway simplifies the seasonal, inter-annual variability and climatic extremes that are known to drive population dynamics and dispersal processes (Walther et al., 2002). Two, we represented climatic analogs based on two water balance variables, but other energyrelated variables such as radiation also drive vegetation dynamics, especially in boreal and some temperate regions (Barichivich et al., 2014;Papagiannopoulou et al., 2017). Three, we portrayed the climate tolerances of vegetation based on a bin width to define climate analogs; the precision of such climate stratifications is known to affect analog computations Hamann et al., 2015). Fourth, we assigned future vegetation based on a majority of potential analogs; agreement among analogs was generally high (~60-100% per PA) but also variable. These necessary choices potentially propagate unknown levels of error and could explain certain incongruities in our results; for example, we project that dominant vegetation could switch back and forth among time periods (notably, oak-forest in the Great Smoky Mountains NP). Our aim here was to help envision the nature and magnitude of possible changes that PAs may face, but we emphasize that our projections reflect uncertainties and should not be interpreted as if they were precise predictions (Stephenson, Millar, & Cole, 2010), especially at the pixel scale.
An inconsistency in thematic resolution of vegetation data (i.e., number of classes) among countries also challenged our interpretations across national boundaries, particularly between Canada and the United States. For example, in the Glacier complex, we project conversion of cold montane conifer forest to dry and mesic conifer forest. Although similar conversions might be expected in the Banff-Jasper complex (Langdon & Lawler, 2015), corresponding distinctions were impossible because all conifer forests in Canada were described with one broad class (i.e., evergreen needleleaf). Furthermore, coarse classification of non-forest vegetation in both the United States and Canada prevented characterization of important shifts among non-forest types. For example, our singular categorization of tundra may mask the ecosystem changes anticipated for those communities (Walker et al., 2006). Indeed, we project only a relatively small degree of change in many high latitude PAs. In addition, our findings are sensitive to the number of vegetation types per country; the greater number of transitions modeled in the United States and México are likely a function of more classes. Despite the constraints imposed by the inconsistencies in thematic resolution, our results consistently point to directional vegetation shifts across PAs over the upcoming decades.
Also, not directly addressed in our framework is that climate change rates will likely exceed that to which vegetation can migrate, creating a disequilibrium between prevailing vegetation composition and changing climate conditions (Thom, Rammer, & Seidl, 2017). Some species and landscape characteristics will infer a slow response or high inertia to climate change. For example, species that are limited by dispersal (e.g., wind dispersed vs. animal dispersed) or that mature slowly may take longer to establish in new environments (Svenning & Sandel, 2013). Climate analogs located at long distances infer a lengthier migration time (i.e., Algonquin Provincial Park, Wabakimi Provincial Park, Thelon Wildlife Sanctuary). Treeline expansion into alpine areas may be highly protracted depending upon microsite conditions (Batllori, Camarero, Ninot, & Gutiérrez, 2009), species characteristics (Harsch & Bader, 2011), or the slow development of ecosystem structure such as soil properties (Skre, Baxter, Crawford, Callaghan, & Fedorkov, 2002;Svenning & Sandel, 2013). Or, individuals of long-lived species such as red fir (Abies magnifica), western white pine (Pinus monticola), and whitebark pine (P. albicaulis) may continue to survive in PAs providing high inertia to the system, although they may not produce viable seedlings in situ (Grubb, 1977).
In contrast, disturbance processes such as wildland fire can serve as catalysts for rapid changes (e.g., opening sites for seed establishment or competitive release), accelerating the response of vegetation to climate change (Thom et al., 2017). The paleo-record demonstrates that fire events have induced significant and rapid vegetation change during periods of high climate velocity in the past (Crausbay, Higuera, Sprugel, & Brubaker, 2017). In the future, disturbance dynamics themselves are likely to intensify with climate change, possibly triggering even more rapid or profound ecological changes. For example, increased fire activity may facilitate conversion of dry forests to non-forest in the western United States , and conifer to deciduous forest and grasslands in boreal ecosystems . Such fire-induced changes may have a greater impact on vegetation than the direct effects of climate change (Bond & Keeley, 2005). Drought stress on trees may increase forest vulnerability to attacks by insects and pathogens (Allen et al., 2010;Bentz et al., 2010). In northern latitudes, thawing of the permafrost layer may trigger ground surface subsidence, killing trees (via waterlogging), and collapsing permafrost plateaus (Baltzer, Veness, Chasmer, Sniderhan, & Quinton, 2014). Such climate-mediated increases in disturbance have the potential to exceed the ecological resilience of forests, inducing broad-scale die-off (Allen et al., 2010), and shifts to non-forest ecosystems as tipping points are crossed (Batllori et al., 2018;Davis et al., 2019;Reyer et al., 2015;Stevens-Rumann et al., 2018;Tepley et al., 2018).
That vegetation could lag behind climatic changes is particularly relevant at the species level. Species demography, seed dispersibility, phenotypic plasticity, interspecific competition, gene flow, and mutation rates may all affect how well individual species track climate (Brubaker, 1986;Corlett & Westcott, 2013;Davis, Shaw, & Etterson, 2005;Sittaro, Paquette, Messier, & Nock, 2017), potentially leading to asynchronous migrations that produce novel species assemblages (Alexander, Diez, & Levine, 2015). As such, the combinations of species assimilating in the future are uncertain and may diverge from the vegetation communities that occupy today's climate space. Therefore, our forecasts may differ from realized shifts, and they are best interpreted as the potential for existing communities to redistribute under future climates over long timescales.

| Conservation implications and future work
The PAs assessed in this study will likely experience fundamental changes in their biogeography, which, in some cases, could compromise their capacity to embody the very ecological characteristics for which they were established. Some parks are already at this precipice, such as Wood Buffalo NP, currently at risk of losing its UNESCO World Heritage status due to climate change and human activities from neighboring areas (Parks Canada, 2019). The fine-scale portraits of vegetation presented here can perhaps serve to better envision the possible future and critical challenges for not only native vegetation but the animal species dependent upon such habitat and food resources. In particular our projections might better enable managers to build conservation strategies and climate adaptation measures for their respective PA. Millar et al. (2007) advocate for a portfolio of management approaches including restraint (i.e., no intervention), resilience/resistance (efforts that sustain ecosystems in the nearterm), and realignment (long-term adaptation). Specifically, intervention may not be needed for some species in large PAs where populations may shift their distributions within the circumscribed landscapes, such as moving upward in elevation (Lawson, Bennie, Thomas, Hodgson, & Wilson, 2014) or to cooler slopes that provide micro-refugia (Maclean, Hopkins, Bennie, Lawson, & Wilson, 2015). However, strategic measures that resist change or build resilience may be needed, such as sustaining highly vulnerable species for as long as possible to buy time for new populations to establish elsewhere (Hansen & Phillips, 2015). For example, planting rust-resistant whitebark pine seedlings may help sustain the species where cold-montane forests are projected to decline. Conversely, if a strong range expansion is indicated, managers could consider strategies to facilitate natural colonization of associated species (i.e., realign). Examples include planting species better adapted to warmer conditions following disturbance (Stephenson & Millar, 2012), or enhancing habitat quality to increase potential dispersers, thereby improving connectivity between PAs (Hodgson, Thomas, Wintle, & Moilanen, 2009;Lawson et al., 2014). Whether species' ranges remain stable, contract, or expand, it is likely that PAs will remain important for conservation under a changing climate (Johnston et al., 2013). Already, the majority of species are disproportionately colonizing PAsevidence that PAs will continue to deliver high biodiversity benefits into the future (Thomas & Gillingham, 2015). PA managers will also need a comprehensive vision for North America's network of protected lands to plan for shielding and connecting important habitats that may provide organisms with safe havens under changing environmental conditions. All PAs, including the large and iconic PAs we studied, interact within a surrounding matrix of unprotected lands, where natural vegetation is increasingly becoming fragmented or reduced (Piekielek & Hansen, 2012). Conservation strategies to facilitate the ability for species to persist, adapt, and migrate across this matrix, might include expanding the sizes of PAs, creating new areas that conserve critical refugia (Michalak et al., 2018), and connecting habitats along species' migratory routes (Carroll, Parks, Dobrowski, & Roberts, 2018). Future work based on the climate analog model could inform such efforts by integrating more species specific information (e.g., dispersal distance), identifying areas important for species movement among PAs (sensu Littlefield, McRae, Michalak, Lawler, & Carroll, 2017) or critical refugia (sensu , and by including habitat fragmentation impacts (Batllori et al., 2017). Conservation actions might in turn be better equipped to cope with the coming changes in the magnitude and direction of species' shifts with climate change.