Forecasting Distributional Responses of Limber Pine to Climate Change at Management-Relevant Scales in Rocky Mountain National Park

Resource managers at parks and other protected areas are increasingly expected to factor climate change explicitly into their decision making frameworks. However, most protected areas are small relative to the geographic ranges of species being managed, so forecasts need to consider local adaptation and community dynamics that are correlated with climate and affect distributions inside protected area boundaries. Additionally, niche theory suggests that species' physiological capacities to respond to climate change may be underestimated when forecasts fail to consider the full breadth of climates occupied by the species rangewide. Here, using correlative species distribution models that contrast estimates of climatic sensitivity inferred from the two spatial extents, we quantify the response of limber pine (Pinus flexilis) to climate change in Rocky Mountain National Park (Colorado, USA). Models are trained locally within the park where limber pine is the community dominant tree species, a distinct structural-compositional vegetation class of interest to managers, and also rangewide, as suggested by niche theory. Model forecasts through 2100 under two representative concentration pathways (RCP 4.5 and 8.5 W/m2) show that the distribution of limber pine in the park is expected to move upslope in elevation, but changes in total and core patch area remain highly uncertain. Most of this uncertainty is biological, as magnitudes of projected change are considerably more variable between the two spatial extents used in model training than they are between RCPs, and novel future climates only affect local model predictions associated with RCP 8.5 after 2091. Combined, these results illustrate the importance of accounting for unknowns in species' climatic sensitivities when forecasting distributional scenarios that are used to inform management decisions. We discuss how our results for limber pine may be interpreted in the context of climate change vulnerability and used to help guide adaptive management.


Introduction
Protected areas are a primary means to conserving biodiversity, but species afforded local protection are often distributed considerably beyond administrative boundaries [1][2][3]. This challenge to place-based conservation and management is further exacerbated by climate and other forms of environmental change that operate over broad spatial scales [4][5][6]. In the US, the Department of the Interior recognized the need for new, landscape-level approaches to resource management and responded by establishing a system of Landscape Conservation Cooperatives (LCCs) [7]. LCCs are intended to facilitate the comanagement of shared species and other resources across jurisdictional boundaries, but management decisions are still ultimately made and enacted by individual management units belonging to agencies and institutions with sometimes very different missions [8], [9]. Given these challenges -and new opportunities afforded by LCCs -managers of protected areas require a detailed understanding of how their key, defining species will respond at multiple scales to ongoing and future environmental change. Here, we forecast possible distributional responses of limber pine (Pinus flexilis) to climate change in Rocky Mountain National Park (Colorado, USA) using species distribution models parameterized at management-relevant scales.
Limber pine is a species of white pine (subgenus Strobus) in North America that influences three major ecosystem processes important to managers: (i) post-fire succession [10], (ii) food provisioning for wildlife [11][12][13], and (iii) snow accumulation and retention [14], [15]. Limber pine is considered a species of management concern [16], [17], primarily due to recent and widespread tree mortality events caused by native mountain pine beetle [18], [19] and invasive white pine blister rust [20], [21]. Climate change is anticipated to interact with and exacerbate these threats in complex ways [22], [23], and it may further shorten the fire return interval in ways that either favor (increased disturbance) or disfavor (unable to reach reproductive age) limber pine, but at a more proximate and tractable level for forecasting it is also expected to exert direct effects on particular dominance and age classes [24][25][26].
Focusing on these direct effects, limber pine has the broadest elevational range of any tree species in Colorado, ranging from 1600 m on the plains to tree-line at 3300 m [27]. The species also occurs throughout most western states [28]. Such broad geographical and elevational limits to distribution raise important and unanswered questions about the spatial extent at which climatic sensitivity (i.e., a species' ability to tolerate a change in climate [29]) should be inferred from limber pine occurrence data and used to model distributional responses to climate change in individual parks and other protected areas. On the one hand, if the traits governing distribution are truly conserved at the species level, inferring climatic sensitivity from limber pine occurrences in and near Rocky Mountain National Park will underestimate the ability of local populations to respond to future change. On the other, if local adaptation or community dynamics largely determine limber pine distribution in and near the park, inferring sensitivity rangewide will overestimate the ability of local populations to respond to future change. Because we do not know the precise mechanisms that govern distributional limits in limber pine, a practical compromise approach for developing forecasts useful to managers is to parameterize models using a range of plausible climatic sensitivities that are inferred from multiple spatial extents [30].
Our analyses bracket local and rangewide spatial extents to evaluate how the size, position, and shape of limber pine distribution in Rocky Mountain National Park may change throughout the 21 st century in response to climate change. We address two questions important to park managers: (i) how long will areas within the current limber pine distribution in the park remain climatically suitable, and (ii) when and where will areas outside this current distribution become more climatically suitable than present? We evaluate two representative concentration pathways (RCP 4.5 and 8.5), designed to capture uncertainty in future greenhouse gas (GHG) emission rates, as expressed by their equivalent radiative forcings in the year 2100 in watts per square meter [31], and ask whether future climate uncertainty or spatial extent has a greater effect on our confidence in distributional responses of limber pine to climate change. We conclude with a discussion of how results may be used to inform assessments of climate change vulnerability [29] and guide possible management strategies.

Species Occurrence Data
Limber pine occurrences were obtained from two data sources that differ in spatial extent: (i) rangewide from the Whitebark and Limber Pine Information System (WLIS) [32], and (ii) locally within Rocky Mountain National Park from the National Park Service (NPS) Vegetation Inventory Program (VIP) [33]. WLIS is a database that compiles from the scientific literature all limber pine observations throughout the species' geographic range. While WLIS is not a single comprehensive survey across the entire geographic distribution of limber pine, the occurrence data compiled by WLIS are broadly distributed and include marginal areas of the species' range with respect to latitude, longitude, and elevation. Approximately 93% of the limber pine range occurs in the US [28] -the maximum geographic extent of our climate data (described below) -and a total of 260 georeferenced WLIS limber pine point localities from the US were considered in the rangewide analysis. Although 14 localities from Canada were excluded, these observations are within 24 to 248 km (mean 68 km) of the northernmost localities in the US (US localities span approximately 1500 km); Canadian localities are also encompassed both longitudinally and elevationally by the US localities. Hence, climates of Canadian localities are thought to be reasonably approximated by the US localities. After removing duplicate point observations at the resolution of the climate data (30 arc-seconds or ,800 m; see below) -to reduce the effects of geographic biases in sampling -a total of 237 spatially unique point localities were used to train the rangewide limber pine models.
VIP data were collected within 1731 km 2 , including Rocky Mountain National Park plus a 1.6 km buffer to the north, west and south, and a 6.4 km buffer to the east, which included an extensive urban interface. In total, 632 vegetation plots were sampled in 2002, where plot size in habitats containing limber pine was 400 m 2 . Limber pine was dominant in two woodland National Vegetation Classification (NVC) types: Limber Pine/ Kinnikinnick Woodland (CEGL000802) and Limber Pine/Common Juniper Woodland (CEGL000807). In the plots, limber pine also commonly co-occurred with Subalpine Fir-Engelmann Spruce (CEGL000986) and Subalpine Fir-Krummholz Shrubland (CEGL000985), but when mapped as a dominant class using aerial imagery it was limited to CEGL000802 and CEGL000807, which encompassed 309 contiguous polygons for a total of 26.81 km 2 (mean: 0.087 km 2 ; 95% CI: 0.004-0.309 km 2 ; estimated accuracy 82.4-96.6%). From these polygons we computed the percentage area in each ,8006800 m climate grid cell (388 grid cells encompassing limber pine; 2646 grid cells encompassing the extent of the park vegetation map) and used these percentages to assign weights to point localities used to train the local limber pine models. Assigning weights in this fashion assumes that grid cells with a higher percentage of limber pine are climatically more suitable or favorable than grid cells with lower percentages, due to fine-scale topographic heterogeneity and the distribution of limber pine suitable land facets that influence microclimates.

Climate Data
Contemporary climate data were obtained from the parameter regression on independent slopes model, PRISM [34]. We selected PRISM 1981-2010 normals (30 arc-seconds, ,800 m spatial resolution) for analysis; this 30-year period provided an estimate of contemporary climate (i.e., smoothing over interannual variability in weather) that was also temporally concomitant with the vegetation inventory. We obtained data for mean minimum monthly temperature, mean maximum monthly temperature, and total monthly precipitation.
Future climate projections were obtained from the NASA Earth Exchange (NEX) Downscaled Climate Projections (DCP) for the conterminous US (NEX-DCP30) [35]. The NEX-DCP30 dataset includes more than 100 climate simulations conducted as part of the Coupled Model Intercomparison Project Phase 5 (CMIP5), downscaled to 30 arc-seconds using the Bias-Correction Spatial Disaggregation (BCSD) approach [36][37][38]. Monthly PRISM data from 1950 through 2005 were used to produce the NEX-DCP30 dataset, as a training reference in BCSD. From this dataset, we selected for analysis a yearly sequence of 30-year normals (2006-2035, 2007-2036 … 2071-2100) designed to encompass both near-and long-term management considerations. We also selected two representative concentration pathways (RCP 4.5 and 8.5) designed to reasonably bracket the lower and upper extremes of modeled changes in temperature [31]. Individual model simulations from a total of 31 general circulation models (GCMs) were downscaled under both RCP 4. From these, we used the ensemble averages for mean minimum monthly temperature, mean maximum monthly temperature, and total monthly precipitation to provide overall consensus estimates of future climate.
The gridded contemporary and future monthly climate variables (30-year normals) were used to calculate a series of more biologically meaningful variables, termed bioclimatic variables. A total of 19 bioclimatic variables were considered in developing the limber pine distribution models [39] (Table 1). Future climate data used in projecting the limber pine models were first calculated as delta surfaces from the NEX-DCP30 1981-2010 baseline; we then added these deltas to the PRISM data used in model training so as to ensure that predicted changes in limber pine distribution were not due to any pixel-level artifacts between the training and projection climate data.

Species Distribution Models
Distribution models of limber pine were developed using MaxEnt [40], version 3.3.3k, although analyses could utilize any of the common presence-only or presence-absence modeling methods [41], as well as more complex dynamic range simulation methods that rely on niche theory [42], [43]. Models were developed at two spatial extents: rangewide (conterminous US) and local (Rocky Mountain National Park). Importantly, this comparison holds resolution constant (,800 m) while varying the size of the geographic domain. Seo et al. [44] show that -holding the spatial extent constant -this grid cell size generally results in higher measures of model performance, compared to models developed at coarser spatial resolutions.
The geographic domain associated with the rangewide model was delineated by all Commission for Environmental Cooperation Level III ecoregions [45] containing the limber pine range [28]. The geographic domain associated with the local model was delineated by the geographic extent of the VIP map. Pseudoabsence data (i.e., background points) were thus constrained to each of these two spatial extents, and we generally accepted default MaxEnt settings as tested and recommended by Phillips & Dudík [46], but we did not render occurrence data spatially unique because: (i) our rangewide point occurrence training data had already been rendered spatially unique at the resolution of the climate data, and (ii) we used the percentage of limber pine area in each ,8006800 m climate pixel to weight point occurrence training data, thus allowing climate pixels with higher cover percentages to be more influential in training the model, under the premise that they also have more suitable microclimates driven by fine-scale topography.
Following model training and comparison, we used MaxEnt to develop logistic probability projections for the top performing models. All projections were restricted in spatial extent to the VIP map. Hence, although models were trained at two spatial extents, and the rangewide model was developed using occurrence data collected throughout the US, both local and rangewide map predictions were evaluated for areas inside and proximate to the park (i.e., the spatial scale of influence for park managers). Temporally, projections were developed for the training period

Statistical Analysis
Our analyses of the size, position, and shape of limber pine distribution over time necessarily required us to render the logistic projections binary. We achieved this using the MaxEnt-provided logistic threshold equating the entropy of thresholded and original distributions, calculated separately for the rangewide (0.159) and local (0.174) models. To further test whether the dominant limber pine class in the park was randomly vs. non-randomly predicted by the rangewide model, we queried the rangewide predictions against all and dominant-only VIP grid cells and compared the mean and confidence interval (CI) of logistic probabilities; evidence of a non-random distribution of dominant logistic probabilities was assessed based on differences between the means and the extent to which CI's were non-overlapping. The two statistical distributions were considerably different, suggesting that dominant limber pine was non-randomly predicted by the rangewide model, so we calculated a new rangewide threshold (0.434) that yielded the same sensitivity (true positive rate) as the dominant model trained at the local extent. For each model6threshold combination (local 0.174, rangewide 0.159, rangewide 0.434), we calculated the area (size), mean elevation (position), and core patch index (shape) of each projection; these metrics in different ways affect the ability of limber pine to influence key ecosystem processes important to managers. Elevations were obtained from the National Elevation Dataset [50]; digital elevation model (DEM) data were obtained originally at 30 m resolution and bilinearly resampled to match the resolution of model predictions. The core patch index is a measure of the degree to which the observed core patch area approaches an idealized patch area, where core area is maximized while controlling for total area. It is calculated as A c /(!A22) 2 , where A c is the total area of core patch, estimated here with an edge width of 1 pixel or grid cell, and A is the total area of the predicted limber pine distribution (i.e., including edge pixels).
We then quantified over time the extent to which limber pine distribution moved or shifted outside the current observed range of the species, as documented and recorded by the vegetation inventory. This was achieved by separately intersecting and taking the union of the binary projections with the grid cells containing dominant limber pine polygons. For each projection, we calculated: (i) area shift as the total area intersect divided by the total area union, (ii) elevation shift as the total elevational range intersect divided by the total elevational range union, and (iii) core patch area shift as the core patch area intersect divided by the core patch area union; all ratios were multiplied by 100 and expressed as percentages.
Finally, we used an additional feature offered by MaxEnt to understand when and where model projections were being extrapolated beyond the training space and into novel climates. The multivariate environmental similarity surface (MESS) measures, for each pixel, the degree to which training and projection environments are the same, where negative values identify dissimilar pixel pairs [51], [52]. Negative MESS values are important in management considerations because they identify areas of extreme uncertainty. For each training spatial extent (rangewide and local) and each RCP (4.5 and 8.5), we intersected each binary projection with its associated MESS and calculated the MESS mean across binary grid cells to estimate if, or when, in the future model projections turned negative.

Results
At the local extent, the full model considering all 19 bioclimatic variables yielded the lowest estimates for AICc (Table 2). Models based on the quarterly variables had the second lowest estimates of AICc, and those based on monthlies had the third lowest, but neither of these yielded acceptably low values for DAICc (Table 2). Hence, all rangewide and local analyses were based on full models.
Current predicted probabilities of limber pine occurrence in and near Rocky Mountain National Park differed dramatically between the models trained at the rangewide and local extents, especially throughout areas where limber pine is not the dominant vegetation class (Figure 1). In these areas, the rangewide model predicted occurrence at relatively high probability (0.1-0.5), whereas the local model predicted occurrence at extremely low probability (,0.05). Within areas where limber pine is the dominant vegetation class, both the rangewide and local models predicted occurrence at high probability (0.5-1.0, Figure 1). This correspondence suggests that, even at a rangewide extent, dominant limber pine is non-randomly distributed in climatic niche space. At the local extent of projection, this pattern is indeed confirmed by the non-random association of logistic probabilities for dominant-only grid cells (mean: 0.540, 95% CI: 0.258-0.716) vs. all grid cells (mean: 0.378, 95% CI: 0.122-0.678). Despite overlapping CIs, the differences are quite pronounced when compared to the logistic probabilities for the local model of dominant limber pine (mean: 0.384, 95% CI: 0.021-0.763), which has a relatively low mean probability due to some marginal grid cells receiving low values ( Figure 1B).
Projecting the distribution of limber pine in the park through 2100, uncertainty surrounding the appropriate spatial extent for training the models is greater than the uncertainty in future climate conditions associated with the use of different RCPs ( Figure 2). However, even in the face of these major sources of uncertainty, one pattern is consistent. Under both spatial extents and both RCPs, the distribution of limber pine in the park is projected to increase in elevation ( Figure 2C). Projections for changes in area ( Figure 2A) and patch fragmentation ( Figure 2E), however, are highly variable and result in trends that occasionally even differ in direction.
Relating these projected changes to the current distribution of dominant limber pine in the park, the spatial extent of model training has a substantial effect on the question of whether, and when, to manage for stasis versus change (Figure 2). According to the model trained at the local extent, the percentage of total area ( Figure 2B) and elevational range ( Figure 2D) decreases to approximately 60% of the current observed distribution by 2035, while core patch area ( Figure 2F) decreases to approximately 20%. By 2100, the RCPs have diverged, and these same measures decrease to approximately 50% (elevational range, RCP 4.5) to 0% (total area, RCP 8.5; core patch area, both RCPs). Meanwhile, according to the rangewide model, the measures decrease to a low of 90% of the current observed distribution by 2035, and by 2100 range anywhere from 100% (total and core patch area, RCP 4.5) to 0% (core patch area, RCP 8.5). If we restrict the rangewide model to dominant limber pine predictions, then the measures decrease approximately 85-0% by 2100. This restriction is reasonable since dominant limber pine in the rangewide model is predicted non-randomly throughout the park ( Figure 1A). Thus, under a scenario that considers this thresholding of the rangewide model, plus the local model, the distribution of limber pine in the park is expected to move upslope and outside of its current elevational range during the 21 st century. Although both the magnitude and rate of change are uncertain ( Figure 2D), a net effect is that dominant limber pine will also shift into new areas ( Figure 2B, F).
Finally, novel future climates did not constitute a major source of uncertainty in the projections. Focusing on areas within the park where limber pine was predicted to occur, according to the different thresholds, mean MESS values from 2035 through 2100 were greater than zero under both spatial extents and RCPs, except for models trained at the local extent and projected to RCP 8.5. Under this model and scenario, mean MESS values turned negative in the year 2091, indicating that projections from this year forward should be interpreted with caution ( Figure 2).

Discussion
A distinct advantage of the species distribution models developed here is that they consider uncertainty in the spatial extent at which climatic sensitivities should be inferred and used to project local responses to climate change in individual parks and other protected areas. When applied to limber pine in Rocky Mountain National Park, results show: (i) the distribution is consistently predicted across models to move upslope in elevation and into new areas, and (ii) magnitudes and rates of change are more uncertain due to questions surrounding the appropriate spatial extent for inferring climatic sensitivity than they are to either future climate scenario (RCP 4.5 or 8.5) or the existence of novel future climates (MESS values). Combined, these results suggest that reducing uncertainty in the correlative model of limber pine distribution is presently more dependent on improving our integration of the species' geographical ecology than on refining our predictions of future environments.
Species often vary by age in their sensitivity to climate [53], [54]. In our analysis of limber pine in Rocky Mountain National Park, we had reason to believe that age classes were mixed within the dominant class, and most importantly, included seedlings and young trees that were recruited during the period of model training . Although older (200-300 year old) limber pine can succumb to drought stress [24], seedlings are an Table 2. Akaike Information Criterion, with finite sample size correction (AICc), for distribution models of limber pine trained at the local extent using six different combinations of bioclimatic variables.   especially climate-sensitive life-stage [25]. In Rocky Mountain National Park, limber pine seedlings were found in 57% of plots surveyed and throughout all five major conifer vegetation types [55]. Hence, recruitment in the park appears to have been ongoing during the 1981 to 2010 period in which climatic sensitivity was inferred, suggesting that -despite the long-lived nature of limber pine at higher elevations [56] -it is reasonable to train a distribution model in contemporary time. An alternative approach would have been to reconstruct past climatic conditions of trees or stands based on their age [24], [57], but this was not possible due to a lack of available dendrochronological data at the necessary spatial extents. Importantly, extension of these analyses to other species should similarly consider both potential dominance and age effects on the forecasts, ensuring that any bias in either factor is accounted for in model interpretation. This is especially important given that disturbance history may not only impact community dominance, but also age class distributions, thus potentially causing age class disequilibrium with respect to climate, which in turn biases climate sensitivities inferred from occurrence data. Our focus on areas where limber pine is the community dominant tree species was motivated by the fact that this structural-compositional vegetation class is relevant to park management because of the ecosystem processes it mediates [10][11][12][13][14][15]. The local model was trained on the dominant class while the rangewide model was indiscriminate with respect to dominance. This difference in dominance between the two models was intentional in that it allowed us to bracket the range of climatic sensitivities inferred from the observed species-climate associations. The low to moderate probabilities of occurrence predicted by the rangewide model throughout most areas of the park are plausible because limber pine is outcompeted by more shadetolerant and rapidly maturing conifers [58][59][60], but exists rather pervasively as a minor component of other structural-compositional vegetation classes [33], [55]. Furthermore, the rangewide model predicted dominant limber pine almost exclusively at high probabilities, suggesting that predictions associated with this dominance class scale geographically and thus offer insights into future distributional responses that are independent of the model trained at the local extent. For managers seeking to integrate forecasts of species' distributional responses to environmental change into management scenarios, the benefits of this are twofold: (i) when inferring environmental sensitivities from species' occurrences, it reasonably brackets plausible estimates of this important source of biological uncertainty, and (ii) when needing to look beyond protected area boundaries (e.g., to broader LCCpartnerships), it provides an opportunity to hierarchically evaluate distributional responses of particular management classes or categories at multiple spatial extents.
Despite opportunities to reduce uncertainty in this fashion in the present analysis, it is noteworthy that elevation emerged as the only distributional metric where predictions were consistent enough across spatial extents and climate scenarios to draw conclusions useful to park managers. This finding underscores the importance of understanding the scales at which species or taxa are in states of distributional equilibrium with respect to the environment, and how such uncertainty can have a dramatic effect on the utility of forecasts made from correlative models. In the case of limber pine, more certain predictions of total and core areas were desirable because these distributional metrics are believed to relate strongly to the potential for establishment in new areas [59], including especially those predicted at higher elevations beyond the current upper elevational limit of the species in the park. In an assessment of climate change vulnerability, such insights are valuable for quantifying the adaptive capacity of limber pine to colonize new environments and areas. However, important insights into adaptive capacity may also be gleaned from local field studies confirming the climatic changes, soils, and land facets needed for upslope establishment in areas that are currently above tree-line [25], [26], [61]. These studies collectively suggest that limber pine seedling establishment is not precluded by alpine soils and microenvironments, but also that moisture stress will impose major limits on upslope recruitment. Under the most extreme climate scenario considered in the present analysis (RCP 8.5), temperatures are projected to increase in the park, but so are annual and most seasonal measures of precipitation (except for Bio 14 (stable) and Bio 18 (decreasing); see also [62]). How exactly these translate to soil moisture content and seedling moisture stress -especially given interactions with nurse objects that further influence seedling microenvironments [63] -merits further study, but mechanistically, the increasing precipitation seen in the ensemble mean of the climate forecasts suggests the possibility that soil moisture conditions will remain -or perhaps become more -favorable. Soils are not presumed to be limiting because soil classes associated with limber pine inside the park also exist at higher elevations in areas that are presently alpine [64]. Of further note are transplantations made at Niwot Ridge, Coloradoimmediately south of Rocky Mountain National Park -that show how limber pine seedlings can maintain key physiological functions at elevations above 3500 m [25]. This elevation is higher than all forecasted mean elevations for the year 2100, except for the most extreme forecasts captured by the local model projected to RCP 8.5, which, due to novel climates, are uncertain. Hence, projected future upslope movements of limber pine above tree-line are plausible given what is known about the highelevation biology of the species.
Future increase in the mean elevation of dominant limber pine in the park is further mediated by a concomitant increase in the lower elevational limit, as evident by a decrease in distributional area. Importantly, this does not equate to a loss of limber pine per se, but rather a loss of dominance status from lower elevations. Furthermore, the loss of dominant limber pine at lower elevations in the park is not explained by sampling constraints preventing immigration because even the rangewide model under RCP 8.5 predicted loss. Mechanistically, such loss is facilitated by succession and limber pine gradually becoming out-competed by more shade-tolerant and rapidly maturing conifers (Douglas-fir, Pseudotsuga menziesii; Engelmann spruce, Picea engelmannii; lodgepole pine, Pinus contorta; subalpine fir, Abies lasiocarpa). According to the vegetation inventory [33], some of this succession at low and high elevations may already be underway. At higher elevations, limber pine occurs in multiple strata of three of the alpine vegetation classes (herbaceous upland alpine .9600 ft, herbaceous upland alpine fellfield, shrub upland alpine). At lower elevations, but within the dominant limber pine vegetation class, all plots and strata containing limber pine also contain one or more of the shade-tolerant competitors (Douglas-fir, Engelmann spruce, lodgepole pine, subalpine fir), and in a subset of these the percentage cover of the competitor species combined represents 50-100% of the percentage cover reported for limber pine. Field photos from limber pine plots that span the compositional gradient further show what successional dynamics might look like on the ground ( Figure 3). Combined, these observations from the vegetation inventory suggest the possibility that limber pine is already in the process of moving upslope, and that projected distributional increases in elevation from 2035 to 2100 are an acceleration of ongoing successional dynamics.
If the projected upslope movements occur, there may be some future need to manage limber pine for change, yet, from a management standpoint, it is useful to place both the magnitude and direction of elevational change in a longer term paleoclimatic perspective. In the early Holocene, limber pine was more common throughout lower elevations (,2500 m) than it is today [65], and many of the now isolated populations are believed to be remnants left behind by recent range retraction from lower elevations [66]. In Colorado, compared to present day, upper tree-line was generally higher in elevation from the early to middle-late Holocene, and generally at comparable or lower elevation over most of the last 3500 years [67]. More recently, since the end of the Little Ice Age (c. 1890), upper tree-lines have been readvancing in elevation [61], [67], [68]. Combined, these past changes in elevation illustrate the degree to which tree-line has been in a state of flux, and suggest perhaps that some of the upslope movement predicted in the park throughout the 21 st century will constitute a re-colonization of alpine habitats by subalpine species. However, it is important to note that most of the paleoclimatic changes in tree-line elevation [67] are smaller in magnitude than those forecasted here, especially for the more extreme model projections. Hence, undoubtedly some of the late 21 st century projections reflect elevational changes that -should they occur -exceed a historical range of variability.
Given the elevational changes forecasted for limber pine in Rocky Mountain National Park, what might be some reasonable and feasible next steps to consider in an adaptive management context? One relates to monitoring and a test of the model prediction that limber pine is moving upslope. Even looking to the year 2035, the magnitudes of elevational change are sufficiently large that one would expect some limber pine recruitment to already be occurring in the lower alpine. The NPS Rocky Mountain Inventory and Monitoring Network has implemented an alpine vegetation monitoring protocol [69]. By c. 2090 under RCP 8.5, climates are projected to favor limber pine establishment in two of the four alpine monitoring sites [70], and forecasts presented here could be used to identify and prioritize additional sites to monitor for possible upslope movements. Another opportunity relates to co-opting management strategies that are more practically developed and implemented in response to some other issue of management concern. In the face of climate change, limber pine management may require efforts to ensure enough individuals establish in diverse land facets or microclimates in order to increase the odds of some persisting as refugia and future seed sources. A very similar management strategy has been advanced in response to white pine blister rust, where the enrichment planting of seedlings from rust resistant parents is considered a potential approach to reducing vulnerability [71]. Hence, implementation of this management strategy to reduce limber pine vulnerability to the rust pathogen would have the added benefit of reducing vulnerability to climate change. Two possible modifications would need to be considered: (i) ensuring a planting coverage of heterogeneous environments that may serve as future climate refugia, and (ii) adoption of seed sources that -in addition to exhibiting rust resistance -have known physiological responses to climate change [25]. This co-opting of management strategies is also noteworthy in that it allows climate change and rust impacts to be jointly managed without having to explicitly forecast their complex yet weakly understood interactions. Ultimately, the value judgments that direct management actions at Rocky Mountain National Park remain to be seen. The modeling presented here offers significant insights into future scenario planning efforts and will assist in dialog necessary to reconcile disparate social values related to active management and intervention in a national park that is 95% designated wilderness.
Finally, stepping back from a focus on limber pine and Rocky Mountain National Park, the conceptual and analytical framework presented here affords new opportunities to evaluate the potential impacts (i.e., exposure6sensitivity [29]) of future climate change on species' distributions at management-relevant scales. Although the use of correlative distribution models to quantify such impacts is not new, an important advancement for management is the ability to explicitly quantify uncertainty introduced by not knowing the spatial extent at which species' climatic sensitivities should be inferred and used to project responses to climate change. In this context, the framework for better linking correlative species distribution models to management needs benefits in full from all of the various modeling methods available [41][42][43]. Applications in other species and geographies have the potential to encourage the use of correlative distribution models in developing 21 st century management scenarios for protected areas, LCCs, and other landscape-scale partnerships.