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.

At each site, ten 1 m² 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.

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

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.

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

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.

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² = 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.

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.

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.

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.