Topographic heterogeneity explains patterns of vegetation response to climate change (1972–2008) across a mountain landscape, Niwot Ridge, Colorado

ABSTRACT

The distributions of biomes worldwide are predicted to shift as vegetation tracks climate change. Ecologists often use coarse-scale climate models to predict these shifts along broad elevational and latitudinal gradients, but these assessments could fail to capture important dynamics by ignoring fine-scale heterogeneity. We ask how the elevational ranges of vegetation types have changed in a mountainous landscape, and investigate the influence of fine-scale topographic, snowpack, and soil properties on vegetation change. We manually classified vegetation from high-resolution repeat aerial photographs from 1972 and 2008 at Niwot Ridge, Colorado, USA, and generally found that trees and shrubs colonized tundra, while tundra colonized barren soils. Only shrubs expanded their elevational range. Several fine-scale topographic, soil and snow characteristics, including elevation, slope, solar radiation, soil bulk density, and interannual snowpack variability, modulated where plant establishment occurred. Each vegetation type had a unique suite of variables best predicting its establishment in new areas. We suggest that fine-scale heterogeneity may strongly control how plants in mountainous regions respond to climate change, and different vegetation types may be sensitive to different aspects of this heterogeneity. An improved understanding of the factors controlling vegetation change gives us a broader understanding of ecosystem response to climate change, nitrogen deposition, and release from grazing.

KEYWORDS:

• Climate change
• alpine vegetation
• topography
• snowpack
• solar radiation

Introduction

Rising temperatures have dramatically altered growing conditions for plant species in many regions of the globe (IPCC 2014; Walther 2003). Global temperatures have increased by 0.6°C since 1951, and are expected to rise between 0.3°C and 4.8°C by the year 2100, depending on emissions scenarios (IPCC 2014). The effects of climate change are amplified in mountain ecosystems, where plants are generally more susceptible to habitat loss (Elsen and Tingley 2015; Engler et al. 2011). In addition to climate change, many sites are also experiencing higher levels of nitrogen (N) deposition, and in some cases, cessation from grazing (Dullinger et al. 2003; Galloway et al. 2008). Because mountain ecosystems have a high degree of topographic heterogeneity, which in turn leads to heterogeneity in microclimates, snowpack, and soil properties, plant responses to environmental change may not be uniform across the landscape (Ackerly et al. 2010; Ford et al. 2013; Geiger, Aron, and Todhunter 2012). For example, over the scale of tens of meters, there can be a shift from a dry, windswept knoll-top with little to no snowpack and shallow, rocky soil, to a leeward knoll-slope with a meters-deep snowpack, and deep, moist, organic soils (Bowman and Seastedt 2001; Bruun et al. 2006). Thus, to make accurate predictions about changes in species distributions, it is crucial to examine how these finer-scale factors mediate responses to environmental change (Austin and van Niel 2011; Barrows and Murphy-Mariscal 2012; Luoto and Heikkinen 2008; Randin et al. 2009; Trivedi et al. 2008).

One of the primary impacts of climate change on mountain ecosystems has been a shift in the distribution of herbaceous plants, shrubs, and trees. Several studies have reported dramatic upward shifts in herbaceous alpine plant species (Grabherr et al. 2010; Lenoir et al. 2008; Parmesan and Yohe 2003; Parolo and Rossi 2008; Peñuelas and Boada 2003; Pockley 2001; Walther, Beißner, and Burga 2005). Shrub encroachment into both arctic and alpine tundra has been widely reported as a response to warming (Elmendorf et al. 2012; Hallinger, Manthey, and Wilmking 2010; Myers-Smith et al. 2011; Sturm, Racine, and Tape 2001; Tape, Sturm, and Racine 2006), and a global manipulative warming experiment provides strong evidence that warming is the driver of this response (Walker et al. 2006). In addition, tree range expansion toward higher elevations has been associated with increases in air temperature (Beckage et al. 2008; Lenoir et al. 2008; Peñuelas and Boada 2003).

Despite these results, shifts in the distribution of vegetation have not been uniform in the context of warming, indicating that other factors are at play. Importantly, fine-scale variables such as topography, snowpack, and soil properties could play a role in modulating vegetation responses to warming, especially in mountainous regions (Bourgeron et al. 2015). It is well established that all of these factors influence plant growth, phenology, and community composition in alpine tundra (Bowman and Seastedt 2001; Dearborn and Danby 2017; Körner 2013; Malanson, Bengtson, and Fagre 2012; Rose and Malanson 2012; Suding et al. 2015; Theobald, Breckheimer, and HilleRisLambers 2017; Vonlanthen, Kammer, and Eugster 2006), subalpine meadows (Loneragan and del Moral 1984), treeline communities (Bader and Ruitjen 2008; Grafius, Malanson, and Weiss 2013; Holtmeier 2009; Weiss, Malanson, and Walsh 2015), and montane forests (McKenzie et al. 2003). Thus, they could also play a role in influencing vegetation change over time, but this has been less well studied.

For example, some studies have found that topographic factors such as the slope and aspect of a site play a greater role than atmospheric temperature in governing the distributional responses of low-growing vegetation (Bennie et al. 2006; Scherrer and Körner 2011). Others have identified plants shifting to lower elevations as a result of climate-based changes in plant water balance (Crimmins et al. 2011). For shrubs, both habitat conditions and biotic interactions make expansion into tundra meadows spatially heterogeneous (Dullinger, Dirnböck, and Grabherr 2013). Lastly, factors such as high wind speeds, which often occur in mountainous regions, can negate the warming effect and prevent shifting of the treeline (Holtmeier and Broll 2010). In other cases, upwards shifts in treelines have occurred, but only in certain aspects or slopes (Danby and Hik 2007; Treml and Chuman 2015). In a global review of treeline studies, which had an average length of 59 years, Harsch et al. (2009) it was found that treelines had only shifted to higher latitudes or elevations in half of the studies examined, and krummholz treelines, typical of elevational gradients, were less likely to have shifted with warming. The non-uniformity of elevational vegetation responses to warming in mountain regions warrants further investigation into the factors that can influence these responses.

Central to the discussion of this variation in responses to climate change is the difference between microclimate and macroclimate. Weather stations do not capture the fine-scale climatic variations at the landscape scale (Ashcroft, Chisholm, and French 2009; Ashcroft and Gollan 2012; De Frenne et al. 2013; Scherrer and Körner 2010), which can be influenced by elevation, radiation, moisture, and exposure (Ashcroft, Chisholm, and French 2008; Lookingbill and Urban 2003). If plants are already confined to specific microclimates, then they will vacate microclimates at the trailing edge at the same rate as they would colonize new microclimates at the leading edge. However, the distribution of alpine microclimates and rates of warming are highly variable (MacLean et al. 2016), and therefore local vegetation responses to climate change differ significantly from predictions by macro-scale models.

In this study, we investigated the extent to which the distributions of high-elevation plant vegetation types (alpine tundra, subalpine forest, shrubs) have responded spatially to climate change over the past four decades. We utilized fine-scale topographic data, satellite-derived snow water equivalent estimates, and field-derived soil measurements to assess the influence of topography, snowpack, and edaphic properties on vegetation change. Because temperature and late melting snow are often thought to be limiting factors at the upper elevation ranges of all these vegetation types, we expected (i) all vegetation types to show movement uphill in response to observed warming over the last several decades. Additionally, we expected (ii) a suite of fine-scale characteristics to also influence vegetation establishment, with greater establishment in microsites with greater exposure to warming (greater solar radiation; MacLean et al. 2016), shallower snowpacks that lengthen the growing season (Franklin et al. 1971), flatter slopes (that experience less disturbance from rock slides and receive more soil moisture that drains from hillslopes; Suding et al. 2015; Walker et al. 1996) and on deeper and less dense soil (reflecting greater soil development, water holding capacity and nutrient concentrations; Aina and Periaswamy 1985; Chaudhari et al. 2013). Thus, we expect that vegetation can move upwards in elevation, but only given certain conditions of other fine-scale variables.

Materials and methods

Study site

Our study was conducted at the Niwot Ridge Long Term Ecological Research (LTER) site and the adjacent valleys to the north (Brainard Lakes) and south (Green Lakes Valley), in the Front Range of the Rocky Mountains, Colorado, USA (Figure 1, 40° 3′ 20′ N, 105° 35′ 22′ W). The site borders the continental divide on its western end, and is about 25 km west of Boulder, CO (Figure 1). Figure 1 There is evidence of human activity at this site, including game-wall systems for hunting, that ranges in age from 7650 to 500 years before present (Bowman and Seastedt 2001). The site was also heavily grazed by sheep in the 1940’s (Bowman and Seastedt 2001). Nitrogen (N) deposition has been increasing over the last few decades to current estimated rates of 6.2 kg N ha⁻¹ yr⁻¹ (Formica et al. 2014; Simkin et al. 2016). Average precipitation from 1952 to 2012 in the alpine at our site was 1090 ± 230 mm yr⁻¹, with a 60 mm yr⁻¹ increase over that time period, driven mostly by increases in winter precipitation (Kittel et al. 2015). Summer temperature data from 1972 to 2008 at our site show a warming trend, especially since 1980 (Figure 2), and overall temperatures also increased between 1989 and 2008 (McGuire et al. 2012). This has led to increased positive degree days and earlier snow meltout times (Caine 2010; McGuire et al. 2012; Preston et al. 2016). Climate at our site is affected by larger scale oscillation patterns in the Pacific, and changes in climate in the 1980’s and 1990’s may be partly attributable to shifts in the Pacific Decadal Oscillation (Kittel et al. 2002). Treeline, or the upper limit of tree life including krummholz (stunted, windblown tree mats) (Tinner and Theurillat 2003), is mostly a gradual krummholz form, where the trees at their uppermost limit are stunted and windblown. Treeline species here are primarily Picea engelmanni (Parry ex Engelmann) (Engelmann Spruce), Abies lasiocarpa (Hook.) Nutt. (Subalpine Fir), and to a lesser extent, Pinus flexilis (James) (Limber Pine). Some P. engelmanni seedlings have been found 50 m above the main closed canopy timberline (Daly and Shankman 1985; Peet 1978). The most abundant shrub is the willow Salix glauca (Linnaeus). Alpine tundra on Niwot Ridge is representative of the region and includes wet, moist, and dry meadow communities, as well as fellfield and snowbed communities (Suding et al. 2015).

Figure 1. Map of the study area and 2000 randomly generated points showing where the three most common vegetation changes occurred as well as points that did not change vegetation type. The points are overlaid on the 2008 orthophoto with hill shading and 100 m contour lines. Bare to tundra n = 39, tundra to shrub n = 27, tundra to open forest n = 26

Figure 2. Mean summer (June, July, August) temperatures from the D1 Meteorological Station (3739 m.a.s.l.) on Niwot Ridge, CO, USA, showing a strong warming trend in the time between the orthophotos (1972, 2008). The linear regression line (red, p < 0.001, R² = 0.31) as well as a loess function (blue), with 95% confidence intervals (shaded) are shown

Aerial images

High resolution (0.6 m) orthographic photographs (orthophotos) taken in 1972 (color-infrared) (DEM) and 2008 (true-color) included an approximately 38 km² region that encompasses subalpine forest, alpine tundra, subnival talus areas, and high peaks, with an elevation range of approximately 3100–4100 m. The photos were topographically corrected and are available on the Niwot Ridge LTER website (niwot.colorado.edu). We characterized ground cover at 2000 randomly generated points across the extent of the orthophotos. These points became our dependent variable in logistic regression models (0 for no vegetation change, 1 for vegetation change). The cover class categories were bare ground, permanent snowpack/glacier, rock, water, alpine tundra, shrub, open forest, and closed canopy forest. Both supervised and unsupervised maximum likelihood classification did a poor job in classifying these particular images into these desired classes, so we classified points manually. We defined open forest as points touching tree vegetation that existed in patches or islands with less than 75% canopy cover within a 25-meter radius of the point. If the point landed on grassy vegetation within an open forest (i.e. a subalpine meadow), the point was classified as tundra vegetation, as these meadows still contain many alpine species. Of the 2000 points, 424 were water, snow, or rock in both years, leaving 1576 points with vegetation (Figure 1). Our analyses focused on 1532 of these points, which overlapped with our snow water equivalent (SWE) dataset (Jepsen et al. 2012, see below).

Predictor variables

Elevation, slope, and solar radiation were derived from a 2 m LiDAR-based digital elevation model (DEM) and conferred to each point using QGIS (QGIS Development Team 2015). Solar radiation was calculated using the area solar radiation tool in Spatial Analyst in ArcGIS 10.2.2 (ESRI 2015). This tool calculates solar radiation based on slope, aspect, and shading from a DEM for a particular date and time. For consistency with our SWE data (see below), we calculated solar radiation at midday on June 1^st, which is also representative of typical exposure to sun during the main part of the growing season in July. Solar radiation has a large impact on temperature and is commonly used to adjust temperature models of mountain landscapes (Daly et al. 2007; Dubayah 1994; Fridley 2009). Because it is based in part on topography, we consider it as a topographic variable.

The SWE data is from a 12-year distributed SWE reconstruction model, which integrates hydrometeorological observations, a distributed snowpack energy balance model, and Landsat-derived snow cover from 1996–2007 (Jepsen et al. 2012). The model back-calculates SWE wherever the snow is deposited, but still suffers errors due to distribution by wind (Jepsen et al. 2012). In general, the dataset adequately identifies parts of the landscape receiving more snow than others. We acknowledge that errors in the model could lead to false negative or false positive results, which should be interpreted cautiously. Snow water equivalent is a measure of the amount of water contained in the snowpack, expressed as a linear depth for each 30 m grid cell. Since late spring snow determines the length of the growing season, which affects plant growth and reproduction (Kirdyanov et al. 2003; Kudo, Nordenhall, and Molau 1999; Totland and Alatalo 2002), we focused on the 12-year mean, inter-annual variation, and 12-year trend (slope of linear regression model), of SWE on June 1^st. Inter-annual variation was calculated as the coefficient of variation (CV, standard deviation/mean) for the 12 years. These variables were calculated for each pixel, and then conferred to our classified points.

Finally, because colonization is most likely to occur near sources of propagules, we calculated a percent cover (of trees, shrubs, and tundra) variable. Cover (proportion of points) and elevation had strong but nonlinear relationships. Thus, we made polynomial models of the relationship between elevation and the percent cover of each vegetation type and then used the coefficients from the model to calculate an approximate percent cover of each vegetation type at each of the 1532 points (Figure A1). This variable, however, can only be viewed as a rough approximation of the amount of viable seed produced by each vegetation type. Some herbaceous plant, shrub, and tree species show decreases in seed production at their range edges (Jump and Woodward 2003; Vaupel and Matthies 2012; HilleRisLambers et al. 2013; Myers-Smith et al. 2011; Kroiss and HilleRisLambers 2015; Buechling et al. 2016). Two years of monitoring seed traps at our site captured no seed production at treeline (Robert Andrus, personal communication). Thus, it is likely that we overpredicted the importance of cover on vegetation change. However, it is still important to include this variable as a basic representation of propagule sources.

In June–August 2015, we ground-truthed vegetation cover at a total of 187 points. Our overall classification accuracy of land cover identity was 92%, with respect to our 2008 classifications. Our sampling scheme involved surveying points in groups of three: a point at which the vegetation cover type had changed based on the orthophotos, a nearby (~ 20–200 m away) point where vegetation had not changed and was the same cover type as the first point in 1972, and a nearby point where vegetation had not changed and was the same cover type as the first point in 2008. In two instances, one point of the three was unable to be reached. At the locations we visited, we measured soil depth by hammering in rebar until it hit bedrock, and collected three soil cores to 10 cm depth to measure gravimetric water content, pH from a 1:2 soil to deionized water slurry, and bulk density. At points with shrubs and trees present, shrub and tree species were identified within a 2 m radius of the GPS point.

Analyses

We examined the minimum and maximum elevation values for each vegetation type in each year to determine if there were range expansions or contractions. Because this could be driven by outliers and not reflect any of the dynamics within the range, we also examined the 5^th and 95^th percentile elevation values for each vegetation type in each year (Zhu et al. 2012). Then, we conducted both univariate and multiple logistic regressions to test our hypotheses about the influence of certain topographic, snowpack, and soil properties on vegetation change. In these models, a value of zero signified any vegetation class in 1972 other than the focal vegetation class that did not change to the focal vegetation class by 2008, while a value of 1 signified any vegetation class in 1972 other than the focal vegetation class that did change to the focal vegetation class by 2008. To be realistic, we limited the forest analysis to points < 3600 m, given that the treeline occurs at ~ 3550 m (Peet 1978; personal observation). The shrub analysis was similarly restricted to points lower than < 3760 m because the highest shrub identified in 2008 was at 3710 m. We limited the tundra analysis to points higher than the 3550 m treeline elevation. The tree and shrub ranges incorporate a 50 m elevation buffer where it is reasonably possible that a tree or shrub could disperse to and establish. This number is reasonable based on previous work on willow seed morphology, spruce and fir dispersal, and the winds at our site (Alexander and Edminster 1983; Noble and Ronco 1978; Uchytil 1992). Models including topography and SWE variables utilized the entire dataset of 1532 points, whereas models including soil factors only included data from the 187 ground-truthed points where soils were sampled.

To test our hypothesis about fine-scale topographic and SWE variables, we ran logistic regressions using an exhaustive, all subsets, AIC selection method (runs all combinations of predictor variables and selects models with the lowest AIC) to identify the best combination of predictors of a vegetation transition (R package ‘bestglm’; McLeod and Xu 2014), which was then used to map predicted probabilities of change across the landscape. We excluded soil variables from this analysis, as we only sampled these variables at the ground-truthed points (12% of all the points). To make our analyses more robust and account for the spatial processes at play during vegetation change, as well as potentially complex distributional relationships and latent interactions, we also ran geographically weighted logistic regressions (GWLR, R package ‘GWmodel’; Gollini et al. 2015), which essentially adds a spatial term to the logistic regression model, and random forest classification models (RFC, R package ‘randomForest’; Liaw and Wiener 2002), which is a nonparametric approach involving decision trees that can take into account interactions and nonlinearity. We ran these models to see if they gave the same results as the logistic regressions.

To specifically test for effects of soil variables on vegetation transitions, we conducted univariate logistic regressions for each predictor variable (soil depth, soil pH, soil bulk density). We interpreted terms with a p-value < 0.05 to be significant predictors, and determined the relationship between the predictor and response variables by the sign of the coefficient from the model.

To help describe the data, rates of change per decade for each vegetation class were calculated using Equation 1(1) $RateofChange=\frac{100\times \mathrm{l}\mathrm{n}\left(N_{2008}÷N_{1972}\right)}{3.6}$(1) , following Dial et al. (2007).

(1) $RateofChange=\frac{100\times \mathrm{l}\mathrm{n}\left(N_{2008}÷N_{1972}\right)}{3.6}$(1)

where 3.6 represents the number of decades and N₂₀₀₈ and N₁₉₇₂ are the number of points of each vegetation class in that year. All analyses were performed with the statistical software R (version 3.4.0, R Core Team 2017).

Results

We found strong evidence of vegetation change across the landscape, especially encroachment by woody vegetation into areas previously characterized by other cover types (Table 1, Figure 1, Figure 3). The largest increase in cover was by shrubs, which increased by nearly 8% per decade. This expansion was driven by colonization of tundra and barren soil by Salix glauca at elevations between 3400–3500 m (Figure 3). For tree cover, open forest increased by 3.9% per decade while closed canopy forest cover increased by a mere 0.6% per decade. The increase in open forest was driven by infilling of tundra in the 3200–3400 m (subalpine tundra meadows) elevation range, mostly by Picea engelmanni and Abies lasiocarpa (Figure 3). Bare ground cover decreased by 4.9% per decade, constituting the largest decrease. Alpine tundra vegetation decreased overall by 0.5% per decade, but increased at higher elevations (Figure 1, Figure 3). Tundra vegetation colonized barren soils at higher elevations, with a 6.2% increase per decade at 3700–3800 m elevation (Figure 3). Plant species richness at ground-truthed plots that switched from barren soil to tundra vegetation (n = 13) ranged from 5 to 11. The forb Geum rossii (R. Br.) Ser. most frequently colonized barren soils. Individuals representing vegetation that had established since 1972 were typically smaller in size and isolated, suggesting that new seedlings established during the time period rather than in situ growth.

Table 1. Land cover change matrix based on classification of orthophotos from 1972 and 2008. Numbers are the number of sample points out of 1576 within each category. Also shown are the rates of change per decade. Bolded numbers highlight the three most common types of vegetation change, which are the focus here

	2008
Land Cover Classification	Bare	Tundra	Shrub	Open Forest	Closed Canopy	Total (1972)
1972Bare	291	39	13	4	0	347
Tundra	0	578	27	26	1	632
Shrub	0	0	91	8	0	99
Open Forest	0	3	0	182	6	191
Closed Canopy	0	0	0	0	307	307
Total (2008)	291	620	131	220	314	1576
Rate of change (% per decade)	−4.9	−0.5	+ 7.8	+ 3.9	+ 0.6

Figure 3. Difference in the percentage of points in each land cover class between 2008 and 1972. At 3700–3800 m, tundra increased from 41.9% of the points to 52.4% of the points, which is a difference in percent of 10.5

Elevation expansion or contraction

No vegetation type experienced a range contraction at the lower end of their elevational range. Tundra vegetation and open forests did not experience range expansion at the upper end of their elevation range, while shrubs moved uphill 14 meters. The distribution of points within the ranges tended to shift upwards in elevation, with all vegetation types showing increases in the 95^th percentile (Table 2).

Table 2. Minimum, 5th percentile, 95th percentile, and maximum elevation (meters above sea level) values for each vegetation type in each year

	Minimum		5^th %		95^th %		Maximum
Vegetation	1972	2008	1972	2008	1972	2008	1972	2008
Closed Canopy	3121	3121	3147	3147	3365	3367	3412	3412
Open Forest	3124	3124	3171	3175	3420	3421	3482	3482
Shrub	3198	3198	3246	3253	3536	3593	3697	3711
Tundra	3141	3141	3255	3256	3801	3813	3972	3972

Topography, SWE, and soil

There were no consistent effects of topographic, SWE, or soil variables on vegetation expansion across vegetation types. Solar radiation was the only predictor variable common to all three of the vegetation types, but the direction of the relationship ranged from positive (tree and shrub) to negative (tundra). Each vegetation transition had its own suite of important predictor variables (Table 3).

Table 3. Predictor variables included in the best-fit logistic regression model for expansion of each vegetationtype. Variables bolded and italicized were also present in both the best geographically weighted logistic regression model and random forest classification model. Italicized variables were only included in one of the other modeling techniques. The ∆AIC was calculated as AIC_{null model} – AIC_{best model} where the null model is an intercept only model. Lower AIC’s indicate better model fits. A decrease in AIC of more than 2 typically means significant model improvement. The R² value is Nagelkerke’s pseudo R² value for logistic regression. AUC is the area under the receiver operating characteristic curve and is a measure of model accuracy that ranges from 0 to 1 with 1 being a perfectly accurate model (no false positives or negatives). SWE = snow water equivalent. SWE trend is the slope of a linear regression model of SWE over 12 years (1996–2007); SWE CV is the coefficient of variation (standard deviation/mean) of SWE

Response Variable	Predictor Variables (direction of relationship)	AIC	∆AIC	R^2	AUC
Open Forest Expansion	Tree Cover(-), Elevation(-), Solar Radiation(+)	318.88	−40.61	0.17	0.81
Shrub Expansion	Shrub Cover(+), Elevation(+), Solar Radiation(+), SWE trend (+)	352.63	−30.98	0.13	0.77
Tundra Expansion	SWE CV(+), Solar Radiation(-), Slope(-)	138.42	−9.84	0.14	0.73

For open forest expansion, tree cover (range = 0 – 93%, b = −0.26, p = 2.70e-05), elevation (range = 3121 – 3600m, b = −0.08, p = 5.71e-05), and solar radiation (range = 11 – 28 WH m⁻², b = 0.27, p = 0.00842) was the best combination of predictor variables for both logistic regression and GWLR (n = 1004). Open forest expansion was more likely in lower elevation areas with less tree cover and higher solar radiation (Figure 4c). Random forest classification models suggest that mean SWE is also an important variable driving open forest expansion (Table A1). Soil characteristics did not appear to influence open forest expansion (n = 42, p > 0.10). However, when just analyzing Engelmann spruce (Logistic regression, n = 42), there was a significant effect of soil bulk density (range = 0.11–0.87g cm⁻³, b = 297.32, p = 0.03), with spruce colonizing areas with higher bulk density.

Figure 4. Predicted probabilities of transitions to tundra (> 3550m), shrub (< 3760m), and open forest (< 3600m) from 1972 to 2008. A 0.1 means there is a 10% chance of that point changing to the vegetation type of the map. Note the different ranges of probabilities. Values are calculated using the intercept and slope coefficients from the best-fit logistic regression model for each (Table 2). Blue areas are lakes

Shrub expansion was best predicted in logistic regression by a combination of shrub cover (range = 0–12%, b = 0.48, p = 3.62e-05), elevation (range = 3121-3760m, b = 0.02, p = 7.13e-05), solar radiation (range = 8–28 WH m⁻², b = 0.14, p = 0.0614), and SWE trend (range = −0.13–0.02 m yr⁻¹, b = 12.51, p = 0.1578) predictor variables (Logistic Regression, n = 1305, Table 3). Shrub expansion was more likely in higher elevation areas with greater shrub cover, greater solar radiation, and increasing snowpack (Figure 4b). In GWLR and RFC models, mean SWE replaced solar radiation as an important predictor variable and slope was also an important variable in GWLR (Table A1). There were no significant effects of any of the soil variables (Logistic regression, n = 48, p > 0.05).

For tundra expansion (n = 215), the best-fit model included SWE CV (range = 0.25–2.59, b = 1.36, p = 0.00142), solar radiation (range = 8 – 28 WH m⁻², b = −0.22, p = 0.01003), and slope (range = 1 – 65°, b = −0.07, p = 0.01648) predictor variables for both logistic regression and GWLR (Table 3). Tundra expansion was more likely in areas with greater interannual variation in snowpack, lower solar radiation, and flatter slopes (Figure 4a). Tundra cover and elevation replaced solar radiation and slope as important predictors in the RFC model (Table A1). There were no significant effects of any of the soil variables (Logistic regression, n = 36, p > 0.05).

Discussion

Over the past four decades, we observed changes in the distributions of several plant communities in conjunction with a directional increase in summer warming. Unlike early climate change predictions for the alpine tundra and treeline (IPCC 1990), we did not observe directional shifts in vegetation moving upslope, with the exception of shrubs. Consistent with other studies, our data suggest an upward expansion of shrubs into tundra vegetation. While tundra vegetation and open forests did not increase their maximum elevations, their distributions within the elevation range shifted uphill, as evidenced by increases in the 95^th percentile of their elevation range. We observed vegetation-specific responses to suites of topographic, snowpack, and edaphic factors. These results suggest that in mountain ecosystems it is crucial to examine fine-scale factors to make accurate predictions about changes in species distributions.

Elevation expansion or contraction

Given that other studies have found uphill responses of vegetation to warming, and that plants are expected to track changing climates, we hypothesized that all vegetation types would show directional uphill shifts. Our results partially support this hypothesis, with all vegetation types increasing at the 95^th percentile. Interestingly, shrubs were the only one of the three vegetation types to experience a true range expansion (uphill shift in its maximum elevation). While tundra plants and trees are establishing in areas in which they were formerly not present (Table 1), this is occurring at the upper end of, but not beyond, their 1972 elevation range, and can thus be considered infilling.

The lack of tree colonization into higher elevation tundra could be due to low winter temperatures, which have not increased over time (McGuire et al. 2012), high wind speeds, ice abrasion, or low soil moisture (Hadley and Smith 1986; Harsch et al. 2009; Holtmeier and Broll 2010; Moyes et al. 2013; Norton and Schoenberger 1984). In addition to these abiotic factors, biotic factors such as herbivory, granivory, and lack of seed production and dispersal can also inhibit treeline shifts (Brown and Vellend 2014; Herrero et al. 2012; Kroiss and HilleRisLambers 2015). Seed trap data from treeline at our site showed no seed rain into the tundra from treeline (Robert Andrus, personal communication). Our results are consistent with the 41 of 50 other alpine krummholz treelines that have not shifted uphill over time (Harsch et al. 2009).

The lack of tundra colonization into higher elevation barren soil could be due to the late melting of snow, poor soil development, nutrient limitation, or lack of microbial mutualists, and testing the effects of these factors on plant growth beyond their range is an important avenue for future research (Chapin et al. 1994; Darcy et al. 2018). The barren soils we collected had high bulk density and low water content, suggesting poor development and water holding capacity. The effects of nutrient limitation on alpine plant range expansion depend on site characteristics such as climate and age. Recent studies in the high alpine of Perú and Alaska suggest that plant colonization is limited by phosphorus, likely due to low weathering rates (Darcy et al. 2018). Sites with greater weathering rates and phosphorus availability could be limited by nitrogen (Raffl et al. 2006), while other sites with developed soils may not be limited by nutrients. Since phosphorus levels and microbial enzyme activity in unvegetated and sparsely vegetated soils at our site are similar to those reported in Perú and Alaska, plants at our site could also be limited by phosphorus in some areas (Bueno de Mesquita et al. 2017; King, Meyer, and Schmidt 2008). We have fertilization experiments in place to test this hypothesis.

Topography, SWE, and soil

We hypothesized that all vegetation types would benefit from the same advantageous topographic, soil, and snowpack conditions. On the contrary, trees, shrubs, and tundra plants each had their own suite of predictor variables that predicted their establishment in new areas.

Solar radiation appears to be an important determinant of the distributions of all three vegetation types. There was a positive relationship between radiation and both tree and shrub expansion, which is consistent with greater tree establishment and shrub growth on more sun exposed slopes (Liu et al. 2015; Stueve et al. 2011). This suggests that higher energy inputs from the sun and the resulting warmer microclimate are important factors for woody plant establishment and growth in the tundra. Solar radiation was also included in the best fit model for tundra expansion, but the direction of the relationship was negative. One potential explanation of this trend is that plant colonization in higher elevation areas beyond intact tundra meadows is more limited by soil moisture, and the higher soil water evaporation on sun exposed slopes (Isard 1986) is thus detrimental to plant establishment. This is consistent with findings on a glacial chronosequence in Austria, where early vegetation developed faster in shaded areas (Raffl et al. 2006).

For trees, negative relationships with both tree cover and elevation suggest that establishment was likely to occur in the lower part of their elevational range in this study, but where open meadows (lower tree cover) were present. This trend of infilling over a period of time with warming summers and increases in precipitation is consistent with trends in Yellowstone National Park (Jakubos and Romme 1993). The increases in winter precipitation at our site could be particularly important, as winter snowpacks are important for insulating tree seedlings over winter (Holtmeier 2009). We expected soil moisture to be beneficial for tree establishment, as low soil moisture has been shown to increase seedling mortality in subalpine fir (Cui and Smith 1991), Engelmann spruce (Hessl and Baker 1997), and limber pine (Moyes, Germino, and Kueppers 2015), in both lower and upper subalpine forests (Moyes et al. 2013). The lack of a relationship is likely due to limitations of our single time point measurement, which does not reflect the water stress that may become more apparent later in the growing season after our sampling. No relationship with SWE variables is in contrast to the suggestion by Franklin et al. (1971) that the snow free period affects tree establishment in subalpine meadows, but is in line with the recent findings of Bader et al. (2018), who found that early snowmelt by two weeks had minimum effects on tree seedling establishment. No effect of soil depth is consistent with other studies (Butler, Malanson, and Resler 2004) and suggests that these trees are able to colonize areas with shallow soil. We found that Engelmann spruce was more likely to colonize areas with higher soil bulk density. This is surprising, given that lower soil bulk density typically is associated with less soil moisture and organic matter and more sand, but may be explained by less competition with tundra plants. Other studies have reported conifer growth on thin, rocky soils as opposed to deeper, finer soils where a thick mat of herbs can inhibit tree establishment (Malanson and Butler 2013; Peet 1988). Another explanation is that the coarser soils facilitated seed trapping and moisture retention around seeds, which facilitated the establishment of subalpine fir in a glacial chronosequence (Jumpponen et al. 1999). Interestingly, we found no relationship between bulk density and subalpine fir establishment, but our results combined with those of Jumpponen et al. (1999) suggest that coarse soils can be beneficial for both of these treeline species.

For shrubs, positive relationships with both shrub cover and elevation highlight their establishment at higher elevations and where there were shrubs to provide propagules. Importantly, we also report a positive relationship between shrub expansion and the trend in SWE. While other studies have reported shrub expansion, to our knowledge this is the first study to connect this expansion to a detailed snowpack dataset at the landscape scale. Areas with increasing snowpack over time appear to be beneficial for shrub establishment, likely because snowpack can insulate soil (Brooks, Williams, and Schmidt 1996; Myers-Smith and Hik 2013), and protect shoots from winter damage (Myers-Smith et al. 2011; Tape, Sturm, and Racine 2006), which has been shown to increase shrub survival (Formica et al. 2014). This result indicates that shrub expansion should be most rapid in areas undergoing simultaneous increases in temperature and snowfall, both of which are occurring at our site (Formica et al. 2014). Furthermore, the rates of N deposition at our site also likely promote shrub growth (Formica et al. 2014), and shrubs also likely started to increase following the cessation of grazing in 1949 (Bowman and Seastedt 2001). The lack of a significant relationship with slope suggests that shrubs are colonizing areas of tundra with both steep and flat slopes, as has been seen in other studies (Tremblay, Lévesque, and Boudreau 2012). Slope was not correlated with SWE in our dataset.

Tundra was more likely to expand on flatter slopes and in areas with greater interannual variation in snowpack. Many of the higher elevation sites in our study landscape are steeper than the main stretch of intact tundra on Niwot Ridge, which creates an unstable landscape with frequent rock slides and snow avalanches that can inhibit establishment (Walker et al. 1996). Our results show that there are areas of shallower slopes at these elevations that are more conducive to plant establishment. The relationship with greater variability in snowpack suggests that plant establishment over time benefits from a combination of years with deep snowpack and years with shallow snowpack. Germinating seeds, for example, may benefit from higher growing season moisture content (Sayers and Ward 1966) provided by later melting snow, while seedlings with established root networks may benefit from earlier melting snow and longer growing seasons by increasing their cover and reproductive output (Galen and Stanton 1993; Kudo 1991; Totland 1997). While we did not measure or include this in our analyses, the presence of larger rocks (> 20 cm diameter) can also create microhabitats suitable for plant establishment in these environments (Jumpponen et al. 1999).

Conclusion

We find that changes in vegetation distribution with time are affected by fine-scale variation, which suggests that fine-scale factors are important in mediating vegetation change in mountainous areas. In particular, the increase in cover and range expansion of shrubs was substantial and should be expected in areas experiencing both increases in temperature and precipitation. Although no one suite of variables was beneficial for all three vegetation types, the factors important for each specific vegetation type are likely to be important in other mountainous regions, but this warrants further study. The lack of a unified response across vegetation types demonstrates how environmental change can lead to a broader reorganization of vegetation. Contrary to our original predictions of similar advantageous microclimates for alpine plants, shrubs, and trees, our data show that each of these vegetation types establishes in different microhabitats as climate changes. One interesting contrast among the vegetation types is that woody species (both trees and shrubs) appear to perform better in warmer, higher energy microclimates (high solar radiation), while flat and cool microclimates are important for herbaceous species moving into unvegetated areas. This result may reflect differences among vegetation types in limitations by energy versus moisture. A key conclusion from our results is some form of heterogeneity matters for all vegetation types, but the importance of topography, versus soil, versus snow depends on type of vegetation. While we only presented findings from one site, our work supports the idea that some plant species may only need to migrate tens of meters to track climate instead of hundreds of meters uphill or hundreds of kilometers poleward (Ford et al. 2013). It is important to note, though, that our results pertain only to a limited number of plant species that we surveyed (i.e., 3 tree species, 1 shrub species, and alpine plant communities). There are certainly many other plant species in our landscape that have not expanded their range and could be experiencing range contractions (Pauli et al. 2007). Local models incorporating topographic heterogeneity lead to drastically different extinction predictions than continent-scale models (Randin et al. 2009). Important future research about the role of fine-scale heterogeneity in distribution modeling should address whether habitat loss has been overpredicted (Austin and van Niel 2011; Barrows and Murphy-Mariscal 2012; Luoto and Heikkinen 2008; Randin et al. 2009) or underpredicted (Trivedi et al. 2008). In any case, incorporating finer scale data will be crucial in shaping local biodiversity conservation planning and management in the face of climate change.

Acknowledgments

Logistical support was provided by the Niwot Ridge LTER program and the CU Mountain Research Station. We thank Dr. William D. Bowman and Dr. Jane G. Smith for helpful suggestions about our fieldwork, and Dr. Lauren M. Hallett for suggestions on the article.

Additional information

Funding

We thank the Niwot Ridge LTER program (NSF DEB 1027341) for funding CPB, the NSF REU program (NSF DBI 1460906) for funding LST, CDB and KCR, and NSF EAR 1141764 for funding NPM. Publication of this article was funded by the University of Colorado Boulder Libraries Open Access Fund.

Appendix

Table A1. Predictor variables included (noted with an X) in the best models of three different techniques, logistic regression (LR), geographically weighted logistic regression (GWLR), and random forests classification (RFC). Variables for GWLR were forward selected following the procedures of the ‘gwr.model.selection’ function in the GWmodel R package and these were compared to the variables selected in the LR model. Variables for RFC were selected using the mean decrease accuracy metric from the ‘importance’ function in the randomForest R package. The top 3 (open forest and tundra expansion) and 4 (shrub expansion) important variables based on this metric were selected to compare to the LR and GWLR results

		Model
Response Variable	Predictor Variable	LR	GWLR	RFC
Open Forest Expansion	Tree Cover	X	X	X
	Elevation	X	X	X
	Solar Radiation	X	X
	SWE Mean			X
Shrub Expansion	Shrub Cover	X	X	X
	Elevation	X	X	X
	Solar Radiation	X
	SWE Trend	X		X
	Slope		X
	SWE Mean		X	X
Tundra Expansion	SWE CV	X	X	X
	Solar Radiation	X	X
	Slope	X	X
	Tundra Cover			X
	Elevation			X

Figure A1. The relationship between elevation and percent cover of the three main vegetation types. Curves are from 3^rd order (tree, R² = 0.99), and 4th order (shrub, R² = 0.84 and tundra, R² = 0.95) polynomial functions