Climate change and micro-topography are facilitating the mountain invasion by a non-native perennial plant species
expand article infoChristian D. Larson, Fredric W. Pollnac, Kaylee Schmitz, Lisa J. Rew
‡ Montana State University, Bozeman, United States of America
Open Access


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.


Alien Plant Invasion, Climatic Effects, High-elevation Refugia, Linaria dalmatica, Species Range Shift


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–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 (1980–2016)? 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?

Materials and methods

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 (1980–2016) 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 m2 monitoring plots were placed randomly in patches of L. dalmatica. During the years of the study (2008–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 variable, 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–2013, 2015) and over a longer period (1980–2016) 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 matrix (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) and both stem density (F = 8.09, Den df = 15, P = 0.012; Fig. 1A) and fruit production (F = 10.77, Den df = 15, P = 0.005; Fig. 1B) increased with elevation.

Figure 1.

A relationship of Linaria dalmatica stem density and B fruit production with elevation. Fruit production regression line is back transformed predicted values derived from the mixed effects model.

Table 1.

The primary environmental variables of the study sites. Results from the linear regression models used to assess the strength and direction of the relationship between PCA axis 1 and the environmental variables used to ordinate the 17 study sites.

Variable Coefficient SE P r 2
Elevation 0.03 < 0.001 < 0.001 0.91
Annual precip. 0.03 0.001 < 0.001 0.82
Summer mean temp -1.80 0.11 < 0.001 0.71
Winter precip. 0.10 0.006 < 0.001 0.70
Autumn mean max temp. -1.72 0.10 < 0.001 0.70
Summer mean min temp. -2.10 0.13 < 0.001 0.68
Annual mean max temp. -1.80 0.12 < 0.001 0.67
Summer mean max temp. -1.42 0.09 < 0.001 0.66
Spring precip. 0.06 0.004 < 0.001 0.65
Autumn mean temp. -2.03 0.14 < 0.001 0.65

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).

Figure 2.

The estimated marginal means for stem density of those sites with shrinking (dark grey) and stable/growing populations (light grey) across the years of the study (2008–2013, 2015).

Table 2.

Assessing Linaria dalmatica stem density (stems m-2) among years for the shrinking and stable/growing populations, using estimated marginal means with Tukey pairwise comparisons. Contrast indicates the years being compared, while estimate is the difference between the mean value of the first year compared with the second.

Population group Contrast Estimate SE df T. ratio P
Shrinking 2008–2012 4.1 1.1 89 3.6 0.01
2008–2013 5.0 1.1 89 4.4 < 0.001
2008–2015 6.1 1.1 89 5.3 < 0.001
2009–2010 3.6 1.2 89 3.1 0.044
2009–2011 3.7 1.2 89 3.1 0.036
2009–2012 5.1 1.2 89 4.3 < 0.001
2009–2013 6.0 1.2 89 5.1 < 0.001
2009–2015 7.0 1.2 89 6.0 < 0.001
2010–2015 3.4 1.1 89 3.0 0.054
Stable/growing 2008–2011 -5.8 1.2 89 -4.80 < 0.001
2008–2012 -4.5 1.2 89 -3.7 0.006
2008–2013 -4.7 1.2 89 -3.8 0.004
2008–2015 -4.5 1.2 89 -3.7 0.006
Table 3.

Assessing Linaria dalmatica stem density (stems m-2) between shrinking and stable/growing populations in each year of the study, using estimated marginal means with Tukey pairwise comparison. Estimate is the difference between the mean value of the shrinking group and the stable/growing group.

Contrast Year Estimate SE df T.ratio P
Shrinking : growing 2008 -0.2 2.5 22.6 -0.1 0.933
2009 -2.3 2.5 23.2 -0.9 0.365
2010 -5.7 2.5 22.6 -2.3 0.034
2011 -8.7 2.5 22.6 -3.5 0.002
2012 -8.8 2.5 22.6 -3.5 0.002
2013 -9.9 2.5 22.6 -4.0 < 0.001
2015 -10.8 2.5 22.6 -4.3 < 0.001

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–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).

Figure 3.

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.

Fixed effects Random effects
Response Predictor Estimate SE df T value P Variance
Site Residual
Stem density (Intercept) 9.71 1.22 15.77 7.94 < 0.001 24.04 ± 4.90 9.49 ± 3.08
Poly(su mean)1 -21.33 7.14 103.57 -2.99 0.004
Poly(su mean)2 -10.21 5.18 114.05 -1.97 0.051
Fruit production (Intercept) 2.55 0.36 11.07 7.03 < 0.001 2.13 ± 1.46 0.78 ± 0.88
Poly(wint precip)1 7.83 1.57 113.99 5.00 < 0.001
Poly(wint precip)2 -5.03 1.15 101.80 -4.38 < 0.001
Poly(su max)1 -12.74 1.87 108.41 -6.80 < 0.001
Poly(su max)2 -3.58 1.39 110.50 -2.57 0.012

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, r2 = 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).

Figure 4.

Comparisons of climate trends (1980–2016) between the sites of the shrinking populations (dark grey) and the stable/growing populations (light grey) for A summer mean temperature B summer mean maximum and C winter precipitation. Dashed vertical line denotes the beginning of the study (2008).

Figure 5.

Principal components analysis among the 17 study sites, based on 27 environmental variables (24 climate variables, elevation, slope angle and aspect). PCA axis 1 explained 80% of the variation and PCA axis 2 explained 13% of the variation. Shrinking sites in each year are dark grey, while stable/growing sites in each year are light grey. Ellipses are standard deviations centred on the mean of each group.


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 (Seipel et al. 2016).

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 sagebrush-steppe 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 CO2 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 CO2 (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 CO2 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 (Lembrechts et al. 2018). 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 (Alexander et al. 2018). 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 (Alexander et al. 2018). 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 precipitation 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.


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 non-native 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.


Funding for this project was provided by the National Research Initiative, project number 2009-55320-05033. LJR is supported by the National Institute of Food and Agriculture, U.S. Department of Agriculture Hatch MONB00363. We would like to thank the United States Forest Service and the National Park Service for their cooperation. We would also like to thank Adam Gebauer, Barbara Keith, Alexandre Wing, Kimberley Taylor, Landon Zimmerman, Jordan Meyer-Morey and Mel Cheeseman for assistance in the field. Finally, we would like to thank the reviewers and the Subject Editor, Johannes Kollmann, for their invaluable insight and comments, which made this manuscript possible.


  • Abatzoglou JT (2013) Development of gridded surface meteorological data for ecological applications and modelling. International Journal of Climatology 33: 121–131.
  • Ackerly D, Loarie S, Cornwell W, Weiss S, Hamilton H, Branciforte R, Kraft N (2010) The geography of climate change: implications for conservation biogeography. Diversity and Distributions 16: 476–487.
  • Adhikari A, Rew LJ, Mainali KP, Adhikari S, Maxwell BD (2020) Future distribution of invasive weed species across the major road network in the state of Montana, USA. Regional Environmental Change 20: e60.
  • Alexander JM, Chalmandrier L, Lenoir J, Burgess TI, Essl F, Haider S, Kueffer C, McDougall K, Milbau A, Nuñez MA, Pauchard A, Rabitsch W, Rew LJ, Sanders NJ, Pellissier L (2018) Lags in the response of mountain plant communities to climate change. Global Change Biology 24: 563–579.
  • Alexander JM, Kueffer C, Daehler CC, Edwards PJ, Pauchard A, Seipel T, Consortium M (2011) Assembly of nonnative floras along elevational gradients explained by directional ecological filtering. Proceedings of the National Academy of Sciences, USA 108: 656–661.
  • Alexander JM, Lembrechts JJ, Cavieres LA, Daehler C, Haider S, Kueffer C, Liu G, McDougall K, Milbau A, Pauchard A, Rew LJ, Seipel T (2016) Plant invasions into mountains and alpine ecosystems: current status and future challenges. Alpine Botany 126: 89–103.
  • Barni E, Bacaro G, Falzoi S, Spanna F, Siniscalco C (2012) Establishing climatic constraints shaping the distribution of alien plant species along the elevation gradient in the Alps. Plant Ecology 213: 757–767.
  • Becker T, Dietz H, Billeter R, Buschmann H, Edwards PJ (2005) Altitudinal distribution of alien plant species in the Swiss Alps. Perspectives in Plant Ecology, Evolution and Systematics 7: 173–183.
  • Blumenthal DM, Norton AP, Cox SE, Hardy EM, Liston GE, Kennaway L, Booth DT, Derner JD (2012) Linaria dalmatica invades south-facing slopes and less grazed areas in grazing-tolerant mixed-grass prairie. Biological Invasions 14: 395–404.
  • Blumenthal DM, Resco V, Morgan JA, Williams DG, LeCain DR, Hardy EM, Pendall E, Bladyka E (2013) Invasive forb benefits from water savings by native plants and carbon fertilization under elevated CO2 and warming. New Phytologist 200: 1156–1165.
  • Boehm AR, Hardegree SP, Glenn NF, Reeves PA, Moffet CA, Flerchinger GN (2021) Slope and aspect effects on seedbed microclimate and germination timing of fall-planted seeds. Rangeland Ecology & Management 75: 58–67.
  • Breshears DD, Huxman TE, Adams HD, Zou CB, Davison JE (2008) Vegetation synchronously leans upslope as climate warms. Proceedings of the National Academy of Sciences 105: 11591–11592.
  • Cannone N, Pignatti S (2014) Ecological responses of plant species and communities to climate warming: upward shift or range filling processes? Climatic Change 123: 201–214.
  • Carboni M, Guéguen M, Barros C, Georges D, Boulangeat I, Douzet R, Dullinger S, Klonner G, van Kleunen M, Essl F, Bossdorf O, Haeuser E, Talluto MV, Moser D, Block S, Conti L, Dullinger I, Münkemüller T, Thuiller W (2018) Simulating plant invasion dynamics in mountain ecosystems under global change scenarios. Global Change Biology 24: e289–e302.
  • Chen I-C, Hill JK, Ohlemüller R, Roy DB, Thomas CD (2011) Rapid range shifts of species associated with high levels of climate warming. Science 333: 1024–1026.
  • Cotto O, Wessely J, Georges D, Klonner G, Schmid M, Dullinger S, Thuiller W, Guillaume F (2017) A dynamic eco-evolutionary model predicts slow response of alpine plants to climate warming. Nature Communications 8: 1–9.
  • Dainese M, Aikio S, Hulme PE, Bertolli A, Prosser F, Marini L (2017) Human disturbance and upward expansion of plants in a warming climate. Nature Climate Change 7: 577–580.
  • Dainese M, Kühn I, Bragazza L (2014) Alien plant species distribution in the European Alps: influence of species’ climatic requirements. Biological Invasions 16: 815–831.
  • Eriksson O (1996) Regional dynamics of plants: a review of evidence for remnant, source-sink and metapopulations. Oikos: 248–258.
  • Geiger R, Aron RH, Todhunter P (2009) The climate near the ground. Rowman & Littlefield, Lanham, 587 pp.
  • Graae BJ, De Frenne P, Kolb A, Brunet J, Chabrerie O, Verheyen K, Pepin N, Heinken T, Zobel M, Shevtsova A (2012) On the use of weather data in ecological studies along altitudinal and latitudinal gradients. Oikos 121: 3–19.
  • Greiser C, Hylander K, Meineri E, Luoto M, Ehrlén J (2020) Climate limitation at the cold edge: contrasting perspectives from species distribution modelling and a transplant experiment. Ecography 43: 637–647.
  • Guo Q, Fei S, Shen Z, Iannone III BV, Knott J, Chown SL (2018) A global analysis of elevational distribution of non-native versus native plants. Journal of Biogeography 45: 793–803.
  • Haider S, Alexander J, Dietz H, Trepl L, Edwards PJ, Kueffer C (2010) The role of bioclimatic origin, residence time and habitat context in shaping non-native plant distributions along an altitudinal gradient. Biological Invasions 12: 4003–4018.
  • Haider S, Kueffer C, Bruelheide H, Seipel T, Alexander JM, Rew LJ, Arévalo JR, Cavieres LA, McDougall KL, Milbau A, Naylor BJ, Speziale K, Pauchard A (2018) Mountain roads and non-native species modify elevational patterns of plant diversity. Global Ecology and Biogeography 27: 667–678.
  • Kalwij JM, Robertson MP, van Rensburg BJ (2008) Human activity facilitates altitudinal expansion of exotic plants along a road in montane grassland, South Africa. Applied Vegetation Science 11: 491–498.
  • Kalwij JM, Robertson MP, van Rensburg BJ (2015) Annual monitoring reveals rapid upward movement of exotic plants in a montane ecosystem. Biological Invasions 17: 3517–3529.
  • Koide D, Yoshida K, Daehler CC, Mueller-Dombois D (2017) An upward elevation shift of native and non-native vascular plants over 40 years on the island of Hawai’i. Journal of Vegetation Science 28: 939–950.
  • Lacey JR, Olson BE (1991) Environmental and economic impacts of noxious range weeds. In: James LF, Evans JO, Ralphs MH, Child RD (Eds) Noxious Range Weeds. Westview Press, Boulder, 5–16.
  • Lázaro-Lobo A, Ervin GN (2019) A global examination on the differential impacts of roadsides on native vs. exotic and weedy plant species. Global Ecology and Conservation 17: e00555.
  • Lembrechts JJ, Alexander JM, Cavieres LA, Haider S, Lenoir J, Kueffer C, McDougall K, Naylor BJ, Nuñez MA, Pauchard A, Rew LJ, Nijs I, Milbau A (2017) Mountain roads shift native and non-native plant species’ ranges. Ecography 40: 353–364.
  • Lembrechts JJ, Lenoir J, Nuñez MA, Pauchard A, Geron C, Bussé G, Milbau A, Nijs I (2018) Microclimate variability in alpine ecosystems as stepping stones for non-native plant establishment above their current elevational limit. Ecography 41: 900–909.
  • Lenoir J, Gégout J-C, Marquet P, De Ruffray P, Brisse H (2008) A significant upward shift in plant species optimum elevation during the 20th century. Science 320: 1768–1771.
  • Lenoir J, Graae BJ, Aarrestad PA, Alsos IG, Armbruster WS, Austrheim G, Bergendorff C, Birks HJB, Bråthen KA, Brunet J (2013) Local temperatures inferred from plant communities suggest strong spatial buffering of climate warming across Northern Europe. Global Change Biology 19: 1470–1481.
  • Lenth R (2020) emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.4.5.
  • Liu Y, Odour AMO, Zhang Z, Manea A, Tooth IM, Leishman MR, Xu X, van Kleunen M (2017) Do invasive alien plant species benefit more from global environmental change than native plants? Global Change Biology 23: 3363–3370.
  • Marini L, Bertolli A, Bona E, Federici G, Martini F, Prosser F, Bommarco R (2013) Beta-diversity patterns elucidate mechanisms of alien plant invasion in mountains. Global Ecology and Biogeography 22: 450–460.
  • McCartney KR, Kumar S, Sing SE, Ward SM (2019) Using invaded-range species distribution modeling to estimate the potential distribution of Linaria species and their hybrids in the US northern Rockies. Invasive Plant Science and Management 12: 97–111.
  • McDougall KL, Lembrechts J, Rew LJ, Haider S, Cavieres LA, Kueffer C, Milbau A, Naylor BJ, Nuñez MA, Pauchard A, Seipel T, Speziale KL, Wright GT, Alexander JM (2018) Running off the road: roadside non-native plants invading mountain vegetation. Biological Invasions 20: 3461–3473.
  • Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’hara RB, Simpson GL, Solymos P, Stevens MHH, Szoecs E, Wagner H (2019) Vegan: Community Ecology Package. R package version 2.5-6.
  • Pauchard A, Kueffer C, Dietz H, Daehler CC, Alexander J, Edwards PJ, Arévalo JR, Cavieres LA, Guisan A, Haider S (2009) Ain’t no mountain high enough: plant invasions reaching new elevations. Frontiers in Ecology and the Environment 7: 479–486.
  • Pepin N, Bradley RS, Diaz H, Baraër M, Caceres E, Forsythe N, Fowler H, Greenwood G, Hashmi M, Liu X (2015) Elevation-dependent warming in mountain regions of the world. Nature Climate Change 5: e424.
  • Pollnac F, Seipel T, Repath C, Rew LJ (2012) Plant invasion at landscape and local scales along roadways in the mountainous region of the Greater Yellowstone Ecosystem. Biological Invasions 14: 1753–1763.
  • Pollnac FW, Maxwell BD, Taper ML, Rew LJ (2014) The demography of native and non-native plant species in mountain systems: examples in the Greater Yellowstone Ecosystem. Population Ecology 56: 81–95.
  • Pyšek P, Jarošík V, Pergl J, Wild J (2011) Colonization of high altitudes by alien plants over the last two centuries. Proceedings of the National Academy of Sciences, USA 108: 439–440.
  • R Core Team (2019) R: A language and environment for statistical computing. R foundation for statistical computing. Vienna.
  • Rew LJ, Brummer TJ, Pollnac FW, Larson CD, Taylor KT, Taper ML, Fleming JD, Balbach HE (2018) Hitching a ride: Seed accrual rates on different types of vehicles. Journal of Environmental Management 206: 547–555.
  • Rew LJ, Maxwell BD, Aspinall R (2005) Predicting the occurrence of nonindigenous species using environmental and remotely sensed data. Weed Science 53: 236–241.
  • Richardson DM, Pyšek P (2006) Plant invasions: merging the concepts of species invasiveness and community invasibility. Progress in Physical Geography 30: 409–431.
  • Robocker WC (1974) Life history, ecology and control of Dalmatian toadflax. Washington Agricultural Experiment Station, Pullman, 79 pp.
  • Seipel T, Alexander JM, Edwards PJ, Kueffer C (2016) Range limits and population dynamics of non-native plants spreading along elevation gradients. Perspectives in Plant Ecology, Evolution and Systematics 20: 46–55.
  • Seipel T, Kueffer C, Rew LJ, Daehler CC, Pauchard A, Naylor BJ, Alexander JM, Edwards PJ, Parks CG, Arevalo JR, Cavieres LA, Dietz H, Jakobs G, McDougall K, Otto R, Walsh N (2012) Processes at multiple scales affect richness and similarity of non-native plant species in mountains around the world. Global Ecology and Biogeography 21: 236–246.
  • Steyn C, Greve M, Robertson MP, Kalwij JM, le Roux PC (2017) Alien plant species that invade high elevations are generalists: Support for the directional ecological filtering hypothesis. Journal of Vegetation Science 28: 337–346.
  • Taylor K, Brummer T, Taper ML, Wing A, Rew LJ (2012) Human-mediated long-distance dispersal: an empirical evaluation of seed dispersal by vehicles. Diversity and Distributions 18: 942–951.
  • Walther G-R, Post E, Convey P, Menzel A, Parmesan C, Beebee TJ, Fromentin J-M, Hoegh-Guldberg O, Bairlein F (2002) Ecological responses to recent climate change. Nature 416: e389.
  • Weed AS, Schwarzländer M (2014) Density dependence, precipitation and biological control agent herbivory influence landscape-scale dynamics of the invasive Eurasian plant Linaria dalmatica. Journal of Applied Ecology 51: 825–834.
  • Whitlock C, Cross W, Maxwell B, Silverman N, Wade A (2017) Montana climate assessment. Montana Institute on Ecosystems, Montana State University, and University of Montana, Bozeman and Missoula, Montana.
  • Wilke C (2019) cowplot: streamlined plot theme and plot annotations for ‘ggplot2’. R package version 1.0.0.
  • Yang M, Pickering CM, Xu L, Lin X (2021) Tourist vehicle as a selective mechanism for plant dispersal: Evidence from a national park in the eastern Himalaya. Journal of Environmental Management 285: e112109.