Climate change and micro-topography are facilitating the mountain invasion by a non-native perennial plant species

Mountainous areas and their endemic plant diversity are threatened by global climate change and invasive species. Mountain plant invasions have historically been minimal, however, climate change and increased anthropogenic activity (e.g. roads and vehicles) are amplifying invasion pressure. We assessed plant performance (stem density and fruit production) of the invasive non-native forb Linaria dalmatica along three mountain roads, over an eight-year period (2008–2015) in the Greater Yellowstone Ecosystem (GYE), USA. We evaluated how L. dalmatica performed in response to elevation, changed over time, responded to climate and how the climate of our sites has changed, and compared elevation, climate, micro-topography (slope aspect and angle), and fruit production among sites with differing temporal trends. Linaria dalmatica stem density and fruit production increased with elevation and demonstrated two temporal groups, those populations where stem densities shrank and those that remained stable or grew over time. Stem density demonstrated a hump-shaped response to summer mean temperature, while fruit production decreased with summer mean maximum temperature and showed a hump-shaped response to winter precipitation. Analysis of both short and long-term climate data from our sites, demonstrated that summer temperatures have been increasing and winters getting wetter. The shrinking population group had a lower mean elevation, hotter summer temperatures, drier winters, had plots that differed in slope aspect and angle from the stable/growing group, and produced less fruit. Regional climate projections predict that the observed climate trends will continue, which will likely benefit L. dalmatica populations at higher elevations. We conclude that L. dalmatica may persist at lower elevations where it poses little invasive threat, and its invasion into the mountains will continue along roadways, expanding into higher elevations of the GYE.


Introduction
Mountainous areas are often regarded as resistant to invasive plant species due to the harsh climate and a low human footprint (Pauchard et al. 2009;Alexander et al. 2016). However, temperatures are rising more rapidly at high elevations (Pepin et al. 2015) and a changing climate, combined with anthropogenic activity, will increase the invasion risk in these areas (Pauchard et al. 2009;McDougall et al. 2018). Non-native plant species abundance is associated with human disturbance and activity (Jauni et al. 2015) and, in mountainous regions, this is commonly roads, their construction and associated vehicular traffic (Pauchard et al. 2009). Seed dispersal by vehicles (Taylor et al. 2012;Rew et al. 2018;Yang et al. 2021) suggests that propagule pressure is likely not limiting along roadsides and the importance of road corridors for dispersal of non-native species has increased attention given to the spread of these species into mountainous areas (Stohlgren et al. 2010;Pollnac et al. 2012;Kalwij et al. 2015;Alexander et al. 2016).
Mountain invasions by non-native plant species generally follow a process called directional ecological filtering (Alexander et al. 2011), where non-native species expand upwards from anthropogenic sources in lowlands, progressively dropping out as the environmental filters associated with increasing elevation become stronger (Alexander et al. 2011;Steyn et al. 2017). Historically, this filtering and the relatively limited anthropogenic activity in mountain and high elevation regions have restricted plant invasions, with non-native plant species richness and abundance decreasing with elevation (Becker et al. 2005;Kalwij et al. 2008; Barni et al. 2012;Seipel et al. 2012;Alexander et al. 2016;Guo et al. 2018). However, global climate change, occurring in mountain and high elevation regions, is weakening the strong climate filter (Pyšek et al. 2011;Carboni et al. 2018), resulting in many species shifting ranges towards higher elevations (Lenoir et al. 2008;Chen et al. 2011). This shift has been demonstrable for non-native plant species (Walther et al. 2002;Kalwij et al. 2015;Koide et al. 2017), which are spreading upwards in elevation twice as fast as native species (Dainese et al. 2017) and anthropogenic propagule sources and roadways are critical for these shifts (Kalwij et al. 2008;Kalwij et al. 2015;Dainese et al. 2017;Lembrechts et al. 2017;Lázaro-Lobo and Ervin 2019). Much of the literature, assessing non-native species along elevation gradients, has used frequency or occurrence data to model species distributions or has described community trends (e.g. species richness) across elevation ranges (Kalwij et al. 2008;Lenoir et al. 2008;Kalwij et al. 2015;Dainese et al. 2017;Koide et al. 2017;Lembrechts et al. 2017;Haider et al. 2018), with few studies describing non-native abundance or reproductive trends over time along mountain elevation gradients (Pauchard et al. 2009;Alexander et al. 2016).
The Greater Yellowstone Ecosystem, USA (GYE) is a mountainous region of high conservation value for its endemic native biodiversity and species of concern. Yellowstone National Park is a conservation area in the region and has management goals of protecting endemic and rare plant species, native biodiversity and limiting and preventing the spread of non-native species (Whipple 2001). Yellowstone National Park and the greater GYE region are attractive tourist destinations accessed almost solely via roadways. Recent research has demonstrated that tourist vehicles in national parks are important seed sources for non-native and fast growing ruderal plant species (Yang et al. 2021), which underscores the risk of roadways facilitating invasions into mountainous conservation areas (Johnston and Pickering 2001;Stohlgren et al. 2010).
The invasive Eurasian perennial forb Linaria dalmatica is a species of concern within the region because of its impacts on forage, land values and native plant communities (Robocker 1974;Lacey and Olson 1991). It was first recorded in GYE in 1957(Whipple 2001, can reproduce vegetatively or via seed and has become widespread throughout the region (Adhikari et al. 2020). Research into the conditions and constraints around the regional invasion of L. dalmatica have found contrasting results. Linaria dalmatica's distribution throughout the region has been best described by summer maximum temperatures (McCartney et al. 2019) or February minimum temperatures (Adhikari et al. 2020), while its growth and abundance have been positively associated with precipitation (Blumenthal et al. 2008;Weed and Schwarzländer 2014). Additionally, studies have found no relationship or a negative relationship between L. dalmatica and elevation (Pollnac et al. 2014;Rew et al. 2005), but positive associations with steep slopes and southern exposures (Rew et al. 2005;Blumenthal et al. 2012).
The importance of roadways, climate, elevation and topography for non-native species, such as L. dalmatica, laid the groundwork for our study, which investigated the temporal and spatial abundance patterns of the L. dalmatica along roadways within the mountains of GYE. Specifically, we assessed the environmental constraints shaping L. dalmatica's abundance and propagule pressure along three roads in the mountains of GYE seven times (2008)(2009)(2010)(2011)(2012)(2013)2015) over an eight-year period (2008-2015). Our research questions were: 1) Does L. dalmatica performance (stem density and fruit production) differ along an elevation gradient? 2) How has L. dalmatica stem density changed over the course of the study? Our remaining questions investigated mechanisms driving the observed elevation and temporal trends: 3a) What climate variables best describe L. dalmatica's performance and 3b) how have these climate variables changed over the course of the study and over a longer term ? 4) Do climate (using the variables that best explained stem density and fruit production), elevation, micro-topography (slope angle and aspect) and fruit production, differ between sites where L. dalmatica populations are stable/growing and where they are shrinking?

Study area
In 2008, L. dalmatica presence was surveyed along three predominately tarred roads near Gardiner, MT, USA (45°01'60"N, 110°42'33"W), one within Yellowstone National Park and the other two just north of the park boundary (Suppl. material 1: Fig. S1A). The road transects proceeded from ca. 1700 m above sea level to the highest elevation extent of the specific roads (2900, 2400, 2200 m). During this initial survey, all L. dalmatica populations within 200 m of the roads were identified and, from the identified populations, six study sites, representing an even spread of elevations along each road, were established. The study sites ranged in elevation between 1736 m and 2316 m. Although our sites were near roads, the community vegetation of our sites was predominantly native and ranged from sagebrush-steppe vegetation at lower elevations to conifer forests at higher elevations. The mean annual temperature  of the study sites was 3.9 °C and the mean annual precipitation was 431 mm. The roads extended beyond the upper elevation range limit of L. dalmatica in the region, ensuring dispersal opportunities; see Pollnac et al. (2014) for more site details.

Field and climate measurements
At each site, ten 1 m 2 monitoring plots were placed randomly in patches of L. dalmatica. During the years of the study (2008)(2009)(2010)(2011)(2012)(2013)2015), L. dalmatica stem density (stems m -2 ) and fruit production (seed capsules m -2 ) were recorded at the end of each growing season (late August) in each monitoring plot. Additionally, in each plot, slope angle was taken with an inclinometer and slope aspect was taken with a compass. We had to remove one site from the study because it was located within a riparian zone that became a construction zone for a period during the project, thereby creating unrepresentative outliers.
Climate data (4 km resolution) used in the analyses were obtained from the University of Idaho's Climate Engine, which uses gridded surface meteorological data, based on PRISM and NLDAS2 (Abatzoglou 2013). Due to proximity, our 18 sites fell into nine climate cells (Suppl. material 1: Fig. S1B), which were collated annually and by season: winter (December-February), spring (March-May), summer (June-August) and autumn (September-November). Autumn climate data used in the analyses were from the year prior to data collection. We downloaded daily minimum, maximum and mean temperatures, as well as daily precipitation and evaluated 24 climate variables: annual variables assessed were mean (°C), mean maximum (°C) and mean minimum (°C) temperatures and precipitation (mm); seasonal variables were the same as the annual variables, additionally, because other studies have determined it to be an important factor for species' growth and distribution (Dainese et al. 2014;Pollnac et al. 2014), the coldest minimum daily temperature (°C) for each season was assessed.

Data analysis
To address our first question, we first averaged performance data from the ten subplots at each site and created linear mixed effects models by assessing stem density and fruit production across the elevation gradient, with elevation as a fixed effect and site as a random effect to account for repeated measures through time.
To address our second question, we first assessed how stem density changed at each site over time using a linear mixed effects model with stem density as the response vari-able, year as a fixed effect and plot as a random effect to account for repeated measures through time. In this analysis, we observed sites with shrinking stem densities (n = 9, Suppl. material 1: Fig. S1C, Suppl. material 2: Fig. S2) and sites with stable or growing stem densities through time (n = 8; Suppl. material 1: Fig. S1C, Suppl. material 3: Fig. S3). We then analysed the shrinking/growing population grouped data through time with a linear mixed effects model. In this model, the response was stem density, the fixed effects were shrinking/growing grouping, year (factor) and the interaction between the two and the random effect was site to account for repeated measures through time. Finally, we calculated the yearly stem density estimated marginal means, also known as least-squares mean, for both groups and the Tukey Method was used for post-hoc pairwise comparisons among years for each group and between groups in each year (R-package "emmeans"; Lenth 2020).
To address the first part of our third question, we evaluated the effects of 24 climate variables on L. dalmatica stem density and fruit production. To do this, following similar methodology as previous studies (Barni et al. 2012;Taylor et al. 2014), we first used our 27 environmental variables (24 climate variables, elevation, slope and aspect) to ordinate our sites using principal components analysis (PCA; "rda" in the R-package "vegan"; Oksanen et al. 2019). Linear regression models were then used to assess the strength and direction of the relationship between the PCA axis 1 and each of the environmental variables used to ordinate the plots. We used the climate variables with the strongest relationship with the PCA axis 1 (Table 1) and created linear mixed effects models for both stem density and fruit production ("lmer" in the R-package "lmerTest"; Kuznetsova et al. 2017). For both response variables, additive and interactive climate models were created with uncorrelated climate data (cor. < 0.7; Suppl. material 4: Table S1) as fixed effects. Performance data were averaged across the ten subplots at each site and site was then treated as the random effect to account for repeated measures through time. Road accounted for very little variance as a random effect and models improved without its inclusion. To evaluate for hump-shaped responses, all climate variables were evaluated and compared with and without a second order polynomial transformation. Model creation was followed by manual stepwise backward model selection by removing insignificant variables, from which a final model was chosen by comparing AIC values (ΔAIC > 2, indicating a better fit model). If there was no difference (ΔAIC < 2) between final models, the most parsimonious was selected.
For the second part of our third question, we used separate linear mixed effects models for the years of the study (2008)(2009)(2010)(2011)(2012)(2013)2015) and over a longer period  using climate variables that were significant in the stem density and fruit production models as responses, year and population group as fixed effects and site as the random effect to account for repeated measures through time.
To address our fourth question, we compared the environment of the shrinking population group with the environment of the stable/growing population group using permutational multivariate analysis of variance with a Euclidean distance ("adonis" in the R-package "vegan"; Oksanen et al. 2019), where the environmental ma-trix (24 climate variables, elevation, slope and aspect) was the response, population group the predictor variable and year as strata to take into account repeated measures through time. We also used the previously-mentioned principal components analysis to visually compare the differences in environments between the two groups. We then individually compared the climate, elevation, slope angle and aspect, and fruit production, between the population groups. Climate of the population groups was compared in the mixed effects models addressing the second part of question three. Mean elevation of the two groups was compared using a linear mixed effects model, with population group as a fixed effect and site as a random effect. For microsite characteristics without repeated measures (slope angle and aspect), an analysis of variance was performed to compare mean values between our groups. We assessed L. dalmatica's association with the slope aspect for the shrinking and stable/growing groups by comparing the mean absolute value of the degrees from south (180°). Finally, we compared fruit production (response) between population groups (fixed effect) using a linear mixed effects model with site as a random effect. All comparisons between populations were made using pairwise comparisons of the estimated marginal means using the Tukey Method.
The linear mixed effects models' assumptions of normality and homoscedasticity were assessed by visually inspecting model residuals and fruit production and temperature data were log normally transformed. Significant relationships between predictor (fixed effects) and response variables at the P < 0.05 level were calculated from F statistics, based on Satterthwaite's approximations of degrees of freedom. All statistics were performed in the statistical software R, version 3.6.0 (R Core Team 2019) and graphics were produced using the packages "ggplot2" (Wickham 2016), "cowplot" (Wilke 2019) and "ggmap" (Kahle and Wickham 2013).

Linaria dalmatica response to elevation
Elevation had the strongest relationship with axis 1 of our principal components analysis, which explained 80% of the variance (Table 1)

Linaria dalmatica's temporal response
In nine of 17 sites, stem density decreased over the eight-year study (shrinking populations), while stem densities at eight sites remained stable (6) or demonstrated a growing (2) trend (Suppl. material 2: Fig. S2). Analysis of this trend on the grouped data (shrinking, stable/growing populations) demonstrated that the stem densities differed between groups (F = 8.62, Den df = 15.0, P = 0.010), that year had an effect on stem density (F = 2.80, Den df = 89.0, P = 0.015) and there was an interaction between the population group and year (F = 11.61, Den df = 89.0, P < 0.001). The temporal trend of the shrinking group demonstrated a consistent and substantial stem density decrease over the course of the study (Table 2; Fig. 2). The stable/growing group's increase over the course of the study was consistent, but more gradual and less in magnitude than the decrease of the shrinking group (Table 2; Fig. 2). Comparing stem densities between the shrinking and stable/growing populations through time demonstrated a divergence as the study progressed: in 2008 and 2009, there were no differences between mean stem densities of the two groups; however, in each subsequent year, the difference between the two groups increased (Table 3; Fig. 2). Results for the ten environmental variables that were most correlated with PCA axis 1 shown in descending order of importance. 'Precip' represents precipitation, 'temp' represents temperature, 'max' represents maximum and 'min' represents minimum.  Linaria dalmatica response to climate Summer mean temperature best described L. dalmatica's stem density response to climate, with stem density demonstrating a hump-shaped response with a peak at 12.5 °C (Table 4; Fig. 3A). Linaria dalmatica fruit production decreased with summer maximum temperature (Fig. 3B), but demonstrated a hump-shaped response to winter  precipitation, with peak fruit production occurring at ca. 150 mm of winter precipitation (Table 4; Fig. 3C).
Over the years of the study (2008)(2009)(2010)(2011)(2012)(2013)2015), summer mean and summer mean maximum temperatures increased (P < 0.001, P < 0.001, respectively). Consistent with this short-term trend, between 1980-2016, summer mean temperatures increased by an average of 0.28 °C per decade (P < 0.001) and summer mean maximum temperatures increased by an average of 0.2 °C per decade (Fig. 4A, B). While winter precipitation demonstrated no temporal trend across the years of our study (P = 0.09), it increased between 1980 and 2016 (P = 0.03), with an average increase of 2 mm per decade (Fig. 4C).

Environmental comparison between shrinking and stable/growing sites
Our permutational multivariate analysis of variance indicated that environments (climate, elevation, slope angle and aspect) differed between the sites of the shrinking populations and the sites of the stable/growing populations (F = 39.0, r 2 = 0.25, P = 0.001). Results of our principal components analysis similarly demonstrated a difference in environments between the two population groups (Fig. 5). There was a substantial trend that summer mean temperatures of the shrinking sites were warmer (15.5 ± 0.71 °C) than the stable/growing sites (13.6 ± 0.66 °C; P = 0.06; Fig. 4A). The shrinking sites had significantly warmer summer maximum temperatures (24.4 ± 0.83 °C) than the stable/growing sites (22.0 ± 0.79 °C; P = 0.05; Fig. 4B) and significantly drier winters (60.6 ± 11.7 mm) than the stable/growing sites (98.7 ± 12.4 mm; P = 0.04; Fig. 4C). Shrinking populations were producing fewer fruit (5.02 ± 1.58) than the stable/growing populations (29.05 ± 8.35; F = 17.7, Den df = 15, P < 0.001). Likewise, the sites of the shrinking and stable/growing populations differed in elevation and micro-topography. The mean elevation of the shrinking populations (1900 m) was 183 m lower than the mean elevation of the stable/growing populations (2083 m; F = 5.33, Den df = 15, P = 0.04). Mean slope angle for the shrinking populations (11.2°) was 4.6° less than the mean slope angle of the stable/growing populations (15.9°; F = 16.7, P < 0.001). Across both groups, L. dalmatica was predominantly found on south-facing slopes; however, the mean degrees from south for shrinking populations was 12.2° greater than the stable/growing populations (F = 4.50, P = 0.04). . Linaria dalmatica stem density response to A summer mean temperature and L. dalmatica fruit production response to B summer mean maximum temperature and C winter precipitation. Fruit production regression line is back-transformed predicted values derived from the mixed effects model. Table 4. Linaria dalmatica response to climate. Results of the linear mixed effects models assessing L. dalmatica stem density (stems m -2 ) and fruit production (seed capsules m -2 ) in response to climate variables. All response variables underwent a second order polynomial transformation.  Precip' represents precipitation, 'su mean' represents summer mean temperature, 'su max' represents summer mean maximum temperature.

Discussion
High elevation and mountainous regions have historically presented more obstacles to plant invasions than lowlands (Alexander et al. 2011;Alexander et al. 2016) and evidence of upward shifts by non-native species has been limited (Cannone and Pignatti 2014); non-native species richness and abundance declines with elevation (Becker et al. 2005;Kalwij et al. 2008;Seipel et al. 2012;Alexander et al. 2016;Seipel et al. 2016;Guo et al. 2018;Haider et al. 2018). However, recently, non-native species have started to shift higher along elevation gradients in mountainous regions (Becker et al. 2005;Kalwij et al. 2008;Kalwij et al. 2015;Dainese et al. 2017;Koide et al. 2017;Lembrechts et al. 2017). While consistent with these studies, our findings are novel because we observed a temporal increase in abundance of a non-native plant species  at higher elevations, while previous research modelled shifts in optimal elevations or species distributions using either floristic diversity (i.e. species richness) or occurrence/ frequency data. Abundance and reproduction of non-native species through time has only rarely been assessed across elevation gradients (Pollnac et al. 2014;Seipel et al. 2016) and findings have demonstrated no relationship between elevation and population growth rate (Pollnac et al. 2014) or a decrease in density and propagule production of two non-native species with elevation . Studies have found non-native roadside mountain invasions to be characterised by species with high temperature affinities (Lembrechts et al. 2017) and broad climate ranges (Alexander et al. 2011;Steyn et al. 2017). We demonstrated that, while L. dalmatica is observed along the elevation gradient, it has relatively cool and wet climate optima, which are likely driving the density and reproduction patterns observed across our elevation gradient. The directional ecological filtering hypothesis (Alexander et al. 2011) predicts that species are generalists that have migrated upslope due to their broad climate niches and predicts that successful high elevation invaders are nested subsets of a lowland species pool (Alexander et al. 2011;Marini et al. 2013;Steyn et al. 2017). However, by evaluating species traits, McDougall et al. (2018) demonstrated that perennial non-native species comprised a higher percentage of non-native species invading higher elevations than annual non-natives. Our findings demonstrate specific climate optima for L. dalmatica that is more abundant at high elevations and less at low elevations, i.e. it is leaning upslope (Breshears et al. 2008).
It has been suggested that non-native species that are shifting upwards in elevation are not in equilibrium with their environment and their expansion is not due to climate change, but is due to filling of potential niches (Kalwij et al. 2015). We agree that niche filling is a mechanism driving invasions that are recent and where dispersal or invasion opportunities have previously been limiting. However, we do not believe this is happening with L. dalmatica at our sites in the mountains of the Greater Yellowstone Ecosystem. Linaria dalmatica was first recorded in GYE in 1957 (Whipple 2001), with historic maps suggesting it was introduced as an ornamental at Mammoth Hot Spring in Yellowstone Park, where it was initially observed spreading mainly along roads (L.J. Rew, unpublished data). Mammoth Hot Springs is geographically close to our sites and we believe with over 60 years of substantial tourist vehicle traffic to facilitate its spread, the species has had ample time and opportunity to achieve equilibrium with the environment along our roads of study.
Biotic characteristics of ecosystems affect non-native plant invasions (Dainese et al. 2014;Richardson and Pyšek 2006) and it is possible that the plant communities of our sites affected the observed differences in abundance. However, we do not believe this is occurring because a previous study, using sites on the same roads as ours, found greater levels of invasion by non-native species in sagebrush habitat compared with conifer forests (Pollnac et al. 2012). Our lower elevation sites were predominantly sagebrushsteppe vegetation, while our higher elevation sites were predominantly conifer forests; therefore, if the biotic characteristics of the sites influenced our results, we would expect higher L. dalmatica in the lower elevation sagebrush communities, and not in the higher elevation conifer forests.
The observed upward shifts in elevation by non-native species have so far been closely associated with anthropogenic sources of propagules and roadways (Kalwij et al. 2008;Kalwij et al. 2015;Dainese et al. 2017;Koide et al. 2017;Lembrechts et al. 2017) with surrounding native communities demonstrating higher resistance to invasion (Haider et al. 2010;McDougall et al. 2018). Similarly, intact native communities are currently resistant to L. dalmatica invasions (Blumenthal et al. 2012) and the risk of invasion along our elevation gradient decreases as the distance from the road increases (Rew et al. 2005). However, warming and elevated atmospheric CO 2 have been shown to favour non-native over native species (Liu et al. 2017) and the L. dalmatica invasion will likely benefit from increases in winter precipitation and atmospheric CO 2 (Blumenthal et al. 2008;Blumenthal et al. 2013). Therefore, while the risk of invasion to interior and intact native communities along our roads is currently low, it will likely increase as the climate warms, winter precipitation increases and atmospheric CO 2 levels continue to rise.

Linaria dalmatica regional responses
Interestingly, our L. dalmatica abundance and reproduction findings are consistent with a recent L. dalmatica species distribution model that found elevation and summer maximum temperature to be key drivers of its regional distribution in Montana, Wyoming and Colorado (McCartney et al. 2019). Similar to our findings, McCartney et al. (2019) concluded that L. dalmatica establishment will likely occur below 2800 m, which is above our highest site. However, McCartney et al. (2019) found a higher summer optimal temperature for L. dalmatica and concluded that its occurrence was associated with warm and dry areas. However, they did not address how their results contrasted with previous findings that found a positive association between L. dalmatica and snowfall and precipitation (Blumenthal et al. 2008;Weed and Schwarzländer 2014), though winter precipitation did remain in their model (McCartney et al. 2019). Our findings were also inconsistent with another recent species distribution model that specifically focused on roadside populations, which found February minimum temperature to be the most important variable for L. dalmatica distribution in the state of Montana (Adhikari et al. 2020). Species distribution models use environmental data and occurrence data to determine relationships and can fail to capture detail, especially with regards to species abundance and reproduction (Greiser et al. 2020 and references therein). This is likely the reason for the discrepancy in results between our study using abundance data in GYE versus the two species distribution models that encompassed the state of Montana (Adhikari et al. 2020) or three states (McCartney et al. 2019). At finer scales, our findings contrasted with a field study conducted along roadsides in the mountains of GYE which found no relationship with L. dalmatica population growth rate and elevation (Pollnac et al. 2014). We believe the difference in L. dalmatica responses to elevation between our study and Pollnac et al. (2014) stem from study duration. Our study was a long-term field study (seven sampling years over an eight year period) and had enough data to cut through the inter-annual noise to perceive trends with more power than the shorter term study of Pollnac et al. (2014).
Fruit production has not been evaluated in other regional studies and we observed an increase with elevation and one of the most important climate variables was winter precipitation. Initially, this is ecologically perplexing; however, L. dalmatica growth and abundance responds positively to snowfall and precipitation (Blumenthal et al. 2008;Weed and Schwarzländer 2014) and we believe our findings indicate ideal growing conditions. Winter precipitation falls as snow in this region, melting in spring to provide ample soil recharge for the plant community. In addition, snow maintains a constant temperature at the soil surface, preventing early season mortality due to freezing events, which could benefit stem survival, thereby increasing seed production. Winter precipitation was not a significant variable for our stem density results because it was correlated with summer mean temperature (> 0.70) and was removed during our model selection process.

Micro-topography: slope angle and aspect
Mountainous regions provide a strong environmental filter, but are also topographically complex at a local scale (e.g. slope angle and aspect), which alters climate and growing conditions (Geiger et al. 2009;Ackerly et al. 2010;Graae et al. 2012;Lenoir et al. 2013) and can increase germination and shape species distributions (Rew et al. 2005;Blumenthal et al. 2012;Boehm et al. 2021). Micro-topography can provide footholds, or refugia, for non-native species above their ideal climate conditions, thereby facilitating plant invasions into higher elevations . Linaria dalmatica's preference for steep south-facing slopes is well established (Robocker 1974;Rew et al. 2005;Blumenthal et al. 2012) and our abundant and increasing populations were at higher elevations and on steeper, more directly south-facing slopes. We speculate that microrefugia at higher elevations could be facilitating L. dalmatica's invasion into the mountains of GYE. Simultaneously, our results may provide evidence that microrefugia are facilitating the persistence of L. dalmatica at sites where the climate is moving away from ideal conditions. We have demonstrated that L. dalmatica has a thermal optimum, which is lower than the warm temperatures of our low elevation sites. While populations at the low elevation sites are shrinking, when compared with growing populations, they are also found on terrain (slope aspect and angle) that reduces solar radiation and temperature, which may be aiding their persistence.

Probable response to continued climate change
Our results present evidence that the climate constraints of a mountainous region are weakening with a warming climate (Pyšek et al. 2011;Carboni et al. 2018). In southwest Montana, temperatures have increased by 0.2 °C per decade since the 1950s (Whitlock et al. 2017) and approximately 1.14 °C since L. dalmatica was introduced. Consistent with this, we found that between 1980 and 2016, summer mean temperatures of our sites have risen by 0.28 °C per decade, while summer maximum temperatures have risen by 0.2 °C per decade. Under different climate scenarios, this warming pattern is predicted to continue (Whitlock et al. 2017). As mean summer temperature best predicted L. dalmatica stem density, with a peak around 12.5 °C and the most abundant and stable populations were found at cooler, higher elevations, the regional warming pattern will likely shift L. dalmatica's optimal climate up in elevation, facilitating L. dalmatica's invasion along the roads of the mountains of GYE.
Our results also suggest that the predicted warming pattern will be detrimental to the L. dalmatica trailing edge populations found at lower elevations. The lower elevation L. dalmatica populations may be in an extinction lag phase, where extinction (at these locations) has been delayed, but is imminent . This delay could be due to the presence of microrefugia and its life history as a perennial plant that can reproduce vegetatively, both which can facilitate the persistence of a species in suboptimal conditions (Eriksson 1996;Cotto et al. 2017). The delay could also be the result of an individual asynchrony with the changing climate . However, as the climate continues to change, the lag phase will pass and populations at these sites will likely vanish. This is consistent with Adhikari et al. (2020), who predicted an overall decline of L. dalmatica in Montana under climate change. Interestingly, regional climate predictions suggest precipitation will increase by mid-century (Whitlock et al. 2017). Consistent with this prediction, the winter pre-cipitation of our sites increased between 1980 and 2016. Fruit production demonstrated a hump-shaped response to winter precipitation, with a relatively high optima. Therefore, if the observed trend continues, L. dalmatica will not only become more abundant and stable at higher elevations, it will also produce more seed, which will only further the invasion.

Conclusions
In the mountains of the Greater Yellowstone Ecosystem, high elevations, especially on steep south-facing slopes, appear more suitable for L. dalmatica, with conditions of many low sites being too warm and dry for the species. Under future climate scenarios, as a warming climate weakens climate restrictions, L. dalmatica will likely expand its range along roads into the mountains of GYE, as has been suggested for other nonnative species under climate change (Dainese et al. 2017;Guo et al. 2018). At lower elevations, due to microrefugia and its perennial life history, L. dalmatica will likely persist for a time, although it is unlikely to spread and will continue to decline, becoming locally extinct with a warming climate.
These findings are important for land managers because they show that high elevation L. dalmatica populations are producing seed and potentially acting as source populations. Therefore, we suggest the focus of L. dalmatica management be on monitoring current populations at high elevations, surveying for new populations on suitable landscape positions and containing high elevation roadside populations before the invasion expands into interior native communities under a changing climate.