Genetically Based Trait Differentiation but Lack of Trade-offs between Stress Tolerance and Performance in Introduced Canada Thistle.

Trade-offs between performance and tolerance of abiotic and biotic stress have been proposed to explain both the success of invasive species and frequently observed size differences between native and introduced populations. Canada thistle seeds collected from across the introduced North American and the native European range were grown in benign and stressful conditions (nutrient stress, shading, simulated herbivory, drought, and mowing), to evaluate whether native and introduced individuals differ in performance or stress tolerance. An additional experiment assessed the strength of maternal effects by comparing plants derived from field-collected seeds with those derived from clones grown in the glasshouse. Introduced populations tended to be larger in size, but no trade-off of stress tolerance with performance was detected; introduced populations had either superior performance or equivalent trait values and survivorship in the treatment common gardens. We also detected evidence of parallel latitudinal clines of some traits in both the native and introduced ranges and associations with climate variables in some treatments, consistent with recent climate adaptation within the introduced range. Our results are consistent with rapid adaptation of introduced populations, but, contrary to predictions, the evolution of invasive traits did not come at the cost of reduced stress tolerance.

Introduction

Invasive species can have profound ecological and economic impacts (Pimentel et al., 2000; Pyšek and Richardson, 2010), so there is considerable benefit to studying the drivers of invasion to inform control and biosecurity programs to mitigate their influence. Assessing trait evolution in introduced species is important because functional traits can influence demographic parameters essential for invasion success. Consequently, there has been considerable interest in identifying trait differences between native and invaded ranges, to try and understand if and when evolutionary changes might causally drive invasion (Thébaud and Simberloff, 2001; Hinz and Schwarzlaender, 2004; Bossdorf et al., 2005; Felker-Quinn et al., 2013; Hodgins et al., 2018).

Although many invasive species are typically found within similar climatic envelopes in the introduced range as the native range (Petitpierre et al., 2012; Atwater et al., 2017), the biotic environment is expected to show substantial differences (Keane and Crawley, 2002; Colautti et al., 2004). Because of this, Blossey and Notzold (1995) hypothesized that escape from natural enemies in the introduced range should drive allocation of resources away from defense to growth or reproduction, resulting in the evolution of invasive traits (evolution of increased competitive ability [EICA]). Specifically, they predicted that introduced populations should evolve to produce more biomass than native populations and that specialist herbivores will experience improved performance on introduced populations relative to native ones. Many common garden experiments comparing native and introduced populations have sought to test the EICA hypothesis (for reviews see Bossdorf et al., 2005; Colautti et al., 2009; Orians and Ward, 2010; Felker-Quinn et al., 2013; Hodgins et al., 2018). Evolutionary shifts in reproduction, growth, defense, and competitive ability have been identified in many cases (Felker-Quinn et al., 2013; Colautti and Lau, 2015), although only size seems to frequently increase during invasion, resulting in limited general support for the EICA hypothesis (Blumenthal and Hufbauer, 2007; Felker-Quinn et al., 2013).

Enhanced resources in the introduced range may also facilitate the evolution of faster growth and greater size, which may come at a cost to abiotic stress tolerance (Grime, 1977; He et al., 2010). Indeed, increased resource availability can enhance invasion (Davis et al., 2000; Blumenthal, 2006; Richardson and Pyšek, 2006) and species abundant in resource-rich habitats often exhibit invasive traits like rapid growth rates (Blumenthal, 2006; Leishman et al., 2007; Rejmanek and Richardson, 2007; Blumenthal et al., 2009). Consequently, as an extension of this hypothesis, we predict that in benign common gardens invasive populations should have improved performance compared with native populations, but under resource-limited conditions this advantage should be lost. Several species in the Asteraceae have enhanced size and reproduction in introduced compared with native populations in control common gardens, while abiotic stressors (mainly drought) have reduced and even reversed this performance advantage (He et al., 2010; Hodgins and Rieseberg, 2011; Turner et al., 2014; Dlugosch et al., 2015; but see Turner et al., 2015). Biotic shifts can also contribute to increased resource abundance in the introduced range facilitating the evolution of invasiveness. This appears to be the case in invasive yellow star thistle where the loss of native species in California grasslands has reduced competition for water resources, facilitating the evolution of larger plants with reduced drought tolerance (Dlugosch et al., 2015).

Studies of rapid trait evolution in invasive species run in opposition to the classic view of evolution as a slow process. However, evolutionary change during invasion can result from both adaptive and non-adaptive processes (Dlugosch and Parker, 2008; Colautti et al., 2009; Colautti and Lau, 2015). Furthermore, shifts in resource allocation are not the only explanation for trait divergence between the native and introduced ranges. For instance, inter- or intraspecific hybridization in the introduced range could lead to enhanced performance (Rieseberg et al., 2007; Schierenbeck and Ellstrand, 2008). Low densities experienced during rapid range expansion can lead to selection for increased investment in reproduction (Burton et al., 2010), potentially contributing to variation within the invaded range as well as differences between the ranges. Similarly, adaptation to heterogeneous environments (e.g., local adaptation to climate) within the native and/or introduced ranges can potentially produce patterns of trait divergence consistent with EICA (Colautti et al., 2009; Colautti and Lau, 2015). Therefore, a broad population sampling, and the inclusion of climate-associated variables in the analysis, are important in assessing likely drivers of trait differentiation between the ranges using common gardens.

Cirsium arvense (L.) Scop. (Cardueae, Asteraceae; also known as Canada thistle) is native to temperate Eurasia (Tiley, 2010) but has become a noxious invasive plant on all remaining continents except Antarctica (Holm et al., 1977; Tiley, 2010; Guggisberg et al., 2012). Canada thistle particularly thrives in open, disturbed habitats under mild climatic conditions (Moore, 1975; Donald, 1994; Tiley, 2010). It is a diploid (2n = 34), dioecious (obligate outcrossing) perennial that spreads both sexually via seeds and vegetatively by creeping roots (Hamdoun, 1972; Moore, 1975; Lloyd and Myall, 1976; Kay, 1985; Donald, 1994; Heimann and Cussans, 1996; Tiley, 2010). This suite of life-history traits, along with the high genetic diversity recorded in both native and invaded ranges (Jump et al., 2002; Hettwer and Gerowitt, 2004; Solé et al., 2004; Slotta et al., 2006, 2010; Guggisberg et al., 2012), may explain why attempts to regulate its spread have so far been largely unsuccessful (Tiley, 2010; Cripps et al., 2011b). The differential response of genotypes to management activities and stress factors (Hodgson, 1964; Hunter and Smith, 1972; Moore, 1975; Guggisberg et al., 2013), however, also indicates that the reasons for its efficient spread around the world are likely manifold.

To test the trade-off hypotheses outlined above, and to look for evidence of rapid adaptation to local habitats in introduced Canada thistle, we asked the following specific questions: (1) is there evidence of genetic differentiation of quantitative traits between native and introduced populations? (2) Do introduced populations have enhanced performance in benign conditions, and is this performance advantage eliminated or reversed under stressful environments? (3) Are there parallel patterns of adaptation to climate within each range? (4) Are trait differences between ranges likely due to maternal effects or are they rather caused by genetic differentiation? To address these questions, we collected seeds from across the introduced North American and the native European range, grew them in benign (control) and stressful conditions (nutrient stress, shading, simulated herbivory, drought, and mowing) in the glasshouse, and measured a number of performance-related traits.

Results

Climate Analysis

We conducted a principal component analysis of the bioclimatic variables downloaded from WorldClim for our sampled populations (Supplemental Tables 1 and 3). Principal component 1 (PC1) explained 56% of the variation in climate (henceforth called CLIMPC1) and was correlated with most of the bioclimatic variables (Supplemental Figures 1 and 2). In particular, temperature seasonality, temperature annual range, and mean diurnal range were strongly positively correlated with CLIMPC1, while minimum temperature of coldest month; mean temperature of coldest quarter; annual precipitation; and precipitation of driest month, driest quarter, and coldest quarter were negatively correlated. CLIMPC1 was also strongly associated with range (F1,39 = 84.29, p < 0.001), with introduced North American populations having higher CLIMPC1 values (mean ± SE: native = −2.59 +/−0.40, introduced = 2.72 ± 0.41), and therefore experiencing greater temperature seasonality and lower temperatures compared with native European populations. PC2 (now termed CLIMPC2) was strongly positively correlated with precipitation of warmest quarter and negatively associated with annual mean temperature and mean temperature of driest quarter. For this principal component, there was no significant difference between the ranges (F1,39 = 0.017, p = 0.90. Mean ± SE: native = −0.043 ± 0.46, introduced = 0.045 ± 0.48).

Germination Rate Differences between the Ranges

We found that germination rate was significantly different between the ranges (F1,40 = 5.03, p < 0.05). Higher germination rate was identified in the native range versus the introduced range (mean ± SE: native = 0.68 ± 0.030, introduced = 0.57± 0.033). However, it was also significantly associated with CLIMPC1 (F1,39 = 4.40, p < 0.05), but not CLIMPC2 or latitude (p > 0.05). Therefore, germination rates were lower in seeds sourced from populations experiencing greater temperature seasonality and lower temperatures, which reflect conditions in the introduced North American range compared with Europe. Once range and CLIMPC1 were both included in the model, the significance of both effects was lost (range: F1,38 = 0.31, p = 0.58. CLIMPC1: F1,38 = 0.50, p = 0.49).

Genetic Differentiation of Reproductive and Vegetative Traits in the Control Common Garden

Using population means for reproductive and vegetative traits, we eliminated traits that were strongly associated (Spearman's ρ2 >0.7) in the control treatment, either retaining traits that were most biologically significant (e.g., fitness-related traits such as flower head number or biomass traits) or randomly choosing traits if they were similar (e.g., leaf number at time point two versus three) when examining pairwise correlations. Following this procedure, we eliminated 20 traits that were strongly associated with other traits in the control treatment but retained 13 traits (Supplemental Figure 3) that were used for the subsequent analysis of the control treatment. We conducted a principal component analysis for this reduced trait set (Supplemental Figure 4). The first two principal components explained 24% and 21%, respectively, of trait variation. PC1 (hereafter called TRAITcPC1) was most strongly positively correlated with above-ground biomass, final stem diameter, and leaf area at time point one and was not significantly different between the ranges (F1,35 = 0.48, p = 0.50). PC2 (hereafter called TRAITcPC2) was most strongly positively correlated with height at time point one, leaf number difference after 3 weeks, and specific leaf area (SLA), and was negatively correlated with final leaf number (Supplemental Figure 4), and was significantly different between the ranges (F1,35 = 4.34, p < 0.05. Mean ± SE: native = 0.60 ± 0.39; mean introduced = −0.51 ± 0.36). In addition, variation appeared to be reduced in both dimensions in the introduced range relative to the native range (Supplemental Figure 4).

The multivariate ANOVA (MANOVA) of vegetative and reproductive traits in the control common garden indicated a significant difference between the native and introduced ranges (p < 0.05; Supplemental Table 4). The difference between the ranges in phenotypes was significant when CLIMPC1 or CLIMPC2 were included as covariates in the analysis (and when they were excluded), and neither climate variable was significantly associated with the multivariate phenotype (Supplemental Table 4). However, for the MANOVA that included latitude, both range and latitude were significant (p < 0.05 in both cases; Supplemental Table 4).

The univariate analyses of the control common garden revealed that greater final (time point t4) stem diameter and leaf number in the introduced range largely drove the differences between the ranges in the MANOVA (Supplemental Table 5). Both of these traits were also positively correlated with latitude in the native and introduced ranges (Table 1; Figure 2). Populations in the introduced range produced larger stems and more leaves at equivalent latitudes in both ranges (Supplemental Table 6). Flower head number declined with increasing latitude in both ranges, although both ranges produced a similar number of flower heads. An increase in below-ground biomass was also observed in the native range with increasing latitude, although this pattern was not evident in the introduced range.

Table 1

Univariate Analysis of Population Mean Trait Responses of Canada Thistle to Range, Latitude, and Their Interaction in the Control Common Garden Using General Linear Models.

	Range	Latitude	Range:Latitude
Flowering day	0.241,34(ns)	0.0621,34(ns)
SLA	3.381,34(ns)	0.591,34(ns)
Leavest1-t0	0.291,34(ns)	1.33,34(ns)
Heightt1	2.361,34(ns)	0.021,34(ns)
Max leaf areat3	0.281,34(ns)	2.931,34(ns)
Shootst4	1.111,34(ns)	1.131,34(ns)
Heightt4	0.081,34(ns)	2.841,34(ns)
Flowers	0.011,34(ns)	11.171,34∗∗
Above biomass	0.021,34(ns)	3.921,34#
Below biomass	4.561,34∗	0.041,34(ns)	4.421,34∗
Max leaf areat1	3.921,34#	0.191,34(ns)	4.001,34#
Stem diametert4	6.851,34∗	5.861,34∗
Leavest4	8.301,34∗∗	7.001,34∗

Type III F-values with degrees of freedom as subscript and symbols specifying significance of effect are reported. Trait descriptions are given in Supplemental Table 2.

Ns, p > 0.1; #p < 0.1; ∗p < 0.05, ∗∗p < 0.01; ∗∗∗p < 0.001.

Figure 2

Trait Divergence in Response to Range and Latitude in Canada Thistle.

Population mean trait responses for final stem diameter (stem diameter t4) (A), final leaf number (leavest4) (B), number of flower heads (flowers) (C), and below-ground biomass (below biomass) (D) in the control common garden with latitude in the introduced North American (peach) and native European (turquoise) ranges of Canada thistle, with model predictions and 95% shaded confidence intervals from stepwise reduced models. Range effect was significant for stem diameter and leaf number (p > 0.05) (see Table 1).

Genetic Differentiation within and between Ranges of Traits in the Control and Treatment Common Gardens

We conducted a principal component analysis of vegetative traits using population means across all three treatments where final vegetative trait measurements were recorded on most individuals (control, nutrient, herbivory). We excluded the time point measured prior to the beginning of the stress (t0). We first eliminated traits that were strongly associated (Spearman's ρ2 > 0.7) and retained six (out of 25) traits (Supplemental Figure 5) for subsequent analyses. We conducted a principal component analysis for this reduced trait set (Supplemental Figure 6). The first two principal components explained 34% and 31%, respectively, of the variation in the traits (hereafter termed TRAITCNHPC1 and TRAITCNHPC2, respectively). The TRAITCNHPC1 was most strongly positively correlated with below-ground biomass, final stem diameter, and final shoot number but negatively correlated with the change in leaf number. We did not detect a significant range by treatment interaction (p = 0.48), but, once the interaction was removed, there was a marginally significant difference between the ranges (F1,119 = 3.76, p = 0.055. Mean± SE: native = 0.35 ± 0.23; mean introduced = −0.30 ± 0.24) and a highly significant treatment effect (F2,119=93.32, p < 0.001). TRAITCNHPC2 was most strongly negatively correlated with final leaf number and leaf difference between time point two and the end of the experiment. There was no significant interaction between treatment and range (p = 0.59), but, when the interaction was removed, we found a significant difference between the ranges (F1,119 = 8.11, p < 0.01. Mean± SE: native = 0.58 ± 0.27; mean introduced = −0.54 ± 0.28) and a highly significant treatment effect (F2,119 = 37.55, p < 0.001).

The analysis of vegetative traits in the control and treatment (herbivory and nutrient stress) common gardens using MANOVA identified a significant effect of range as well as treatment (p < 0.001; Supplemental Table 7). No significant interaction between range and treatment, which could be an indicator of a potential evolutionary trade-off associated with abiotic or biotic stress response, was identified.

Univariate analysis of vegetative traits in the control and treatment (herbivory and nutrient stress) revealed that final leaf number, final stem diameter, and the change in leaf number from time point two to the end of the experiment had a significant range effect, and no interaction with treatment was evident (Figure 3; Table 2). In these cases, the introduced populations had significantly larger trait means than the native populations. This pattern was maintained even when initial leaf number was used as a covariate to control for differences in size at the start of the stress treatments (Supplemental Table 9). Biomass and stem diameter in the herbivory and nutrient stress were significantly lower than the control, indicating a significant fitness impact of the stressors (Supplemental Table 10). Leaf and shoot production were also significantly lower than the control in the nutrient stress. However, this was not the case in simulated herbivory where half of each new leaf was removed, perhaps indicating that the simulated herbivory was compensated for by increased production of leaves and shoots.

Figure 3

Trait Divergence in Response to Treatment and Range in Canada Thistle.

Population least square means trait responses for change in leaf number between time point two to the end of experiment (leavest4-t2) (A), final leaf number (leavest4) (B), and final stem diameter (mm) (stem diametert4) (C) to range and treatment in three common gardens between introduced North American (peach) and native European (turquoise) ranges of Canada thistle, along with 95% confidence intervals. Range and treatment were significant for all traits (p < 0.05) (see Table 2).

Table 2

Univariate Analysis of Population Mean Trait Responses of Canada Thistle to Range and Treatment (Control, Nutrient, and Herbivory Treatments).

	Range	Treatment
SLA	0.031,119(ns)	25.522,119∗∗∗
Shootst4	0.011,119(ns)	19.852,119∗∗∗
Leavest4-t2	7.561,119∗∗∗	48.682,119∗∗∗
Leavest4	10.741,119∗∗	32.032,119∗∗∗
Below biomass	0.271,119(ns)	37.742,119∗∗∗
Stem diametert4	18.141,119∗∗∗	171.852,119∗∗∗

Interactions terms were non-significant and therefore removed. Type III F-values with degrees of freedom as subscript and symbols specifying significance of effect are reported. Trait descriptions are given in Supplemental Table 2.

Ns, p > 0.1; #p < 0.1; ∗p < 0.05, ∗∗p < 0.01; ∗∗∗p < 0.001.

To determine if the apparent trait differences between the ranges reflected patterns in climate adaptation within each range, we ran two separate MANOVAs where we included the first two climate PCs (CLIMPC1 or CLIMPC2) in the analysis. With the inclusion of these climate PCs in the MANOVA, range remained significant (Supplemental Tables 7 and 8). Furthermore, we did not identify a significant range by treatment interaction or a three-way interaction between range, treatment, and either CLIMPC1 or CLIMPC2. The same pattern was also found when leaf number prior to the start of the stress was used as a covariate to control for initial size differences among the ranges (Supplemental Table 8). These patterns do not support the hypothesis of an evolutionary trade-off in response to stress, even when potential differences in local climates are accounted for. However, we did identify interactions between the climate PCs and range, indicating an association between traits and climate variables that differed between the introduced and native range.

As the MANOVAs provided significance for all three main effects and the interaction between range and the climate PCs, we conducted univariate analyses to further dissect the cause of the significant effects. For most traits, associations with CLIMPC1 were not evident (Table 3; Supplemental Figure 7). For final stem diameter (stem diameter t4), the univariate analysis that included CLIMPC1 produced non-significant effects for range and CLIMPC1, although when CLIMPC1 was removed range was highly significant (see Table 2). We identified two traits with a significant CLIMPC1 by range interaction. For final leaf number (leavest4), there was a significantly positive association with CLIMPC1 in the introduced North American range, but no significant pattern in the native European range (Supplemental Table 11). Consequently, leaf number increased in populations with historically greater temperature seasonality and lower temperatures but only in North America. The inclusion of CLIMPC1 also resulted in no significant overall effect of range, and a lack of significant difference between the ranges for equivalent climates where the climate variable overlapped between the ranges (Supplemental Table 11). Finally, the change in leaf number between time point two and the end of experiment (leavest4-t2) was significantly negatively associated with CLIMPC1 in the native range, but this pattern was not significant in the introduced North American range. Here, the change in leaf number declined with temperature seasonality and colder historic temperatures. However, introduced populations had greater leaf production overall and for equivalent overlapping climates (Supplemental Table 11).

Table 3

Univariate Analysis of Population Mean Trait Responses of Canada Thistle to Range, Treatment (Control, Herbivory, Nutrient), CLIMPC1, and Their Interactions.

	Range	Treatment	CLIMPC1	Range: CLIMPC1
SLA	2.771,118#	26.072,118∗∗∗	3.581,118#
Shootst4	1.581,118(ns)	20.032,118∗∗∗	2.101,118(ns)
Leavest4-t2	22.031,117∗∗∗	55.552,117∗∗∗	0.121,117(ns)	11.191,117∗∗
Leavest4	2.841,117#	33.922,117∗∗∗	6.371,117∗	7.781,117∗∗
Below biomass	1.131, 118(ns)	38.271, 118∗∗∗	2.681, 118(ns)
Stemt4	1.701,118(ns)	172.961, 118 ∗∗∗	1.77 1,118(ns)

Trait descriptions are given in Supplemental Table 2. Type III F-values with degrees of freedom as subscript and symbols specifying significance of effect are reported.

ns p > 0.1; #p < 0.1; ∗p < 0.05, ∗∗p < 0.01; ∗∗∗p < 0.001.

Similar to the patterns with CLIMPC1, CLIMPC2 had few significant associations with traits in the univariate analysis (Supplemental Table 12). For SLA, final leaf number (leavest4), and change in leaf number between time point two and the end (leavest4-t2) of experiment, we identified a significant range by CLIMPC2 interaction. SLA and the change in leaf number were significantly associated with CLIMPC2 (positively correlated with precipitation of warmest quarter and negatively associated with mean diurnal range and mean temperature of driest quarter) but only in the native range (Supplemental Table 13). Specifically, CLIMPC2 increased as SLA declined and the change in leaf number increased. By contrast, final leaf number was positively associated with CLIMPC2 in both ranges. We identified greater final leaf number and change in leaf number in the introduced range for equivalent CLIMPC2 values, but only at lower CLIMPC2 values (minimum European CLIMPC2).

Survival in the Drought, Mowing, and Shade Common Gardens

We conducted a survival analysis in three separate stress treatments (drought, mowing, shade) to assess if individuals from native and introduced populations showed differences in mortality rates over time. We found no significant difference between the native and introduced ranges in survival for the drought (χ21 = 2.88), mowing (χ21 = 0.02), and shade (χ21 = 0.13) common gardens using a Cox proportional hazards regression (Supplemental Table 14; Supplemental Figure 8).

Maternal Effects

Using the MANOVA we found a significant experiment by treatment interaction, but no interaction involving range and experiment, suggesting there was no evidence that maternal effects were driving patterns among the ranges. However, we also did not identify a main effect of range, perhaps because the maternal effects experiment only used a small number of populations, limiting power. We identified an experiment by treatment interaction (Supplemental Table 15). Univariate analysis of the three traits used in the MANOVA identified the same treatment by experiment interaction in one (change in leaf number) of the three traits (Table 4; Supplemental Figure 9). For the change in leaf number, the maternal effects and main experiment were equivalent in the control, but, in the nutrient stress, the change in leaf number was much reduced in the maternal effects garden resulting in a larger difference between the nutrient and the control for the glass-house-derived clones. For above- and below-ground biomass, there was a significant treatment and experiment effect, but no interaction of treatment with experiment. The below- and above-ground biomass was lower in the maternal effects experiment and in the nutrient stress. However, the small population sample size limited our capacity to readily detect differences in the greenhouse-derived clones versus the field-collected seed.

Table 4

Population Mean Trait Responses to Range, Treatment and Experiment to Identify Maternal Effects on Traits.

	Range	Treatment	Experiment	Experiment:Treatment
Above biomass	0.011,28 (ns)	113.421,28∗∗∗	77.311,28∗∗∗
Leavest2-t1	1.141,27	2.571,27	0.121,27	6.751,27∗
Below biomass	2.611,28(ns)	10.431,28∗∗	21.661,28∗∗∗

The interactions were tested in a stepwise manner and removed if they were not significant, starting at the highest order interaction. We reported type III F-values with degrees of freedom as subscript and symbols specifying significance of effect. Trait descriptions are given in Supplemental Table 2.

ns, p > 0.1; #p < 0.1; ∗p < 0.05, ∗∗p < 0.01; ∗∗∗p < 0.001.

Discussion

Our results show that plants derived from the invaded range of Canada thistle are significantly larger in terms of leaf and shoot production than those from the native range when grown in a common garden. In addition, several phenotypic traits were found to co-vary with latitude or climate in a similar way across both ranges, suggestive of parallel local adaptation. However, we failed to find evidence to support our hypothesis that evolutionary trade-offs between performance and tolerance to abiotic or biotic stress contribute to invasion success in the introduced range, as introduced populations tended to perform equivalently or better than native populations across treatments. This is despite the fact that we tested a wide range of stress treatments in controlled common gardens.

Evidence for Genetically Based Trait Differentiation of Traits between Ranges

Many previous studies have reported that plants from populations in the introduced range of a species are larger, grow faster, and are more fecund than those from conspecific populations in the native range (reviewed in Felker-Quinn et al., 2013; Colautti and Lau, 2015). While some such reports may not be reliable because of a failure to account for sources of within-range variation (e.g., latitude) or maternal effects, the overall pattern remains (Colautti et al., 2009). Our results are consistent with this previous work; after controlling for latitude and climate variation, we found evidence of significant phenotypic differentiation between native and introduced populations. Therefore, the trait differences that we observed between the ranges are unlikely to be due solely to local adaptation to differential climate conditions within each range. While the differentiation between the ranges was mainly driven by significant increases in stem diameter and leaf number in populations from the introduced range, it is noteworthy that, in most cases, for performance-related traits (e.g., biomass and flower head production), means in the introduced populations exceeded those in the native populations. Lastly, our experiments showed that trait variation was not significantly affected by maternal environment, although our maternal effects experiment was likely underpowered. Thus, we can conclude that the differentiation in quantitative traits reported here is likely genetically based.

Lack of Evolutionary Trade-Offs in Introduced Populations

The most prominent explanations for the widely reported pattern of increased growth rate and reproductive output of introduced relative to conspecific native populations involve evolutionary trade-offs with tolerance to biotic (Blossey and Notzold, 1995) or abiotic stress (Grime, 1977). While strong and significant treatment effects were observed in our experiments, no significant trade-offs were identified. In fact, populations from the invaded range significantly out-performed native populations in both the herbivory and nutrient stress treatments; trait means for introduced populations were higher than for native populations in most cases, and significantly so for final leaf number, final stem diameter, and change in leaf number over time. No significant differences in survivorship were seen between introduced and native populations under mowing, drought, and shading, although there was a trend toward greater survivorship of native populations under mowing stress. If improved performance of introduced populations in benign environments evolved through increased resource allocation to growth and reproduction at the expense of stress tolerance, we would predict reduced performance of introduced compared with native populations in at least some of the stressful environments. However, our data do not support this hypothesis.

Evolutionary trade-offs with tolerance to abiotic and/or biotic stresses have previously been reported for other Compositae weeds, possibly contributing to invasion success. For example, trade-offs with performance under drought stress have been demonstrated in spotted knapweed (He et al., 2010), common ragweed (Hodgins and Rieseberg, 2011), weedy sunflowers (Mayrose et al., 2011), and star thistle (Dlugosch et al., 2015). On the other hand, invasion success of diffuse knapweed has been attributed to an escape from the trade-offs with drought stress that were observed for plants from the native range (Turner et al., 2014). Trade-offs with biotic stress, consistent with the EICA hypothesis, have also been reported in some Compositae weeds, including weedy sunflowers (Mayrose et al., 2011), dandelion (González-Teuber et al., 2017), star thistle (Montesinos et al., 2019), and common ragweed (Fukano and Yahara, 2012; Sun and Roderick, 2019). However, other studies of Compositae weeds have failed to fully support the EICA hypothesis (e.g., Meyer and Hull-Sanders, 2008; Turner et al., 2014; van Boheemen et al., 2019a).

So why do we not detect abiotic or biotic trade-offs with tolerance to either abiotic or biotic stress in Canada thistle? One possibility is that there is a trade-off with survivorship under mowing stress, but our experiment lacked the power necessary to demonstrate this. Another possibility is that we studied the wrong stresses. The EICA hypothesis predicts that introduced plant species will be released from native specialist predators, but here we simulated herbivory, which might not accurately represent the effects of specialist herbivores (Blossey and Notzold, 1995). For example, leaf trichome/spine densities have been shown to reduce the performance of a specialized herbivore on Canada thistle (Cripps et al., 2015), but we did not assay variation in leaf trichome/spine densities in our experiment. Alternatively, it might be that tolerance to different stresses in Canada thistle is largely uncoupled from performance, similar to that reported for introduced populations of diffuse knapweed (Turner et al., 2014). Such decoupling is often reported in cultivated germplasm (e.g., Fereres et al., 1986; Ali et al., 2007; Gao et al., 2019) and is prized by plant breeders. Mechanistically, this can occur in different ways. Some stress resistance traits may be able to evolve independently of general resource use and growth potential or yield (Chapin et al., 1993; Munns, 2011; Sadras and Richards, 2014). Alternatively, decoupling can be associated with inducible resistance responses (Sato et al., 2016), which are less costly than constitutive resistance mechanisms. This explanation is supported by transcriptomic studies, which show that native and introduced populations of Canada thistle differ with regard to inducible R-gene defense and sensitivity in response to abiotic stresses (Guggisberg et al., 2013).

Signs of (Parallel) Local Adaptation

While the focus of our study was on genetic differences between the native and invaded range, we also found extensive, genetically based phenotypic variation within each range. Some traits, such as final leaf number, stem diameter, and number of flower heads, displayed parallel latitudinal clines in both ranges, which represents strong evidence of local adaptation. Several other traits (e.g., below-ground biomass) were correlated with latitude in one of the ranges only. Final leaf number was the only trait that was significantly correlated with climate in both ranges, although two other traits (SLA and change in leaf number) were found to be associated with climate in the native range. The finding of significant trait evolution in the invaded range is not unusual (reviewed in Hodgins et al., 2018), and offers additional evidence that phenotypic evolution in introduced populations can be surprisingly fast, even in perennials such as Canada thistle, perhaps reflecting the high levels of genetic variation present in the introduced range of this outcrossing species (Guggisberg et al., 2012).

Comparisons with Previous Studies on Canada Thistle

Early studies on Canada thistle reported contradictory results regarding the enemy-release hypothesis. On the one hand, comparative field surveys between Europe (native range) and New Zealand (introduced range) failed to detect significant changes in plant performance between the ranges, despite reduced herbivory in New Zealand (Cripps et al., 2010). On the other hand, a common garden experiment aimed at comparing the fitness of European (native) and North American (introduced) populations of Cirsium arvense in the absence of natural enemies in the home range suggested that introduced genotypes always produced more above- and below-ground biomass than their native counterparts, irrespective of the nutrient regime, which is consistent with our finding of increased performance of introduced populations for some traits (Abela-HofbPauerová and Münzbergová, 2011). In the former, plants were not grown in a common garden setup, while the latter may not be representative, for only two populations were investigated per range, only one reciprocal transplant was investigated, and data were not corrected for climatic differences. As a result, neither study is sufficient to support or refute EICA.

Cripps et al. (2019) detected significant genetic variation in response to complete defoliation after three rounds of mowing. We did so as well in terms of survival, but failed to detect significant differences in mortality across the two ranges and therefore cannot ascertain that Canada thistle has evolved along a fitness trade-off. Since Canada thistle repeatedly reallocates its resources between shoot and root, meaning that shoots alternatively serve as source for root elongation or as sink for stem and leaf production (Donald, 1994; Leathwick and Bourdôt, 2012; Verwijst et al., 2018), the timing and strength (i.e., frequency) of mowing are of utmost importance if differential responses between native and introduced populations are to be expected. Irrespective of these considerations, the fact that only half of the individuals died over the course of our experiment indicates that sustained efforts need to be deployed to sufficiently deplete sinks and hence completely eradicate established Cirsium arvense.

Aside from the fact that Canada thistle seems highly tolerant to defoliation, the circumstances where enemy release can be expected and how strong this effect has to be to become advantageous should be considered. Cirsium arvense encountered a large guild of putative enemies upon its introduction to North America, because this continent sustains 92 indigenous Cirsium species with their suite of Cardueae-specialized herbivores (Hill and Kotanen, 2010; Cripps et al., 2011a, 2011b). However, such a hypothesis does not hold for other regions, such as New Zealand, where no native thistles were initially present, and even a transient enemy release may suffice for a successful establishment, until (unfruitful) biological controls are deliberately released (Nunes and Kotanen, 2018). A recent study by Nunes and Kotanen (2017) further indicates that below-ground herbivory has a stronger (negative) impact on plant performance than above-ground herbivory and should therefore be examined more deeply in the future.

Relative Importance of Non-adaptive Trait Evolution in Introduced Range

During invasion, neutral processes can drive patterns of trait differentiation within and between ranges. Founder events and population bottlenecks associated with introductions are expected to reduce variation present in introduced populations and facilitate divergence, at both genetic and phenotypic levels. However, for highly quantitative traits, introduction bottlenecks are expected to lead to more limited change relative to molecular markers (Hodgins et al., 2018) as quantitative genetic variation should be less affected by the loss of rare alleles occurring during bottlenecks (Lewontin, 1965). Furthermore, common garden data of quantitative traits often identify similar levels of differentiation in the native range as the introduced range (Colautti and Lau, 2015) and there is frequent evidence for adaptive trait divergence occurring during invasions (Bock et al., 2015). To control for non-adaptive evolution of traits, ancestry using molecular markers can be assessed to partition the amount of trait variation due to non-adaptive evolution between the native and introduced populations (Schrieber et al., 2017; van Boheemen et al., 2019b; McGoey et al., 2020). Although a study of invasion history using microsatellite markers has been performed in this species (Guggisberg et al., 2012), the same populations were not examined between studies, so this analysis was not possible. However, the molecular data have revealed the presence of an introduction bottleneck in North American Canada thistle, although high levels of genetic variation are maintained within and among introduction populations. This likely reflects the outcrossing mating system of the species and multiple introductions from both western and eastern Europe (Guggisberg et al., 2012). Consequently, it is likely that quantitative trait variation in North America is not strongly limited by invasion history and patterns of divergence between ranges, and trait differentiation with respect to latitude and climate variables reflects to some degree adaptation to local conditions encountered during range expansion. However, future analysis of reciprocal transplant or common garden experiments incorporating population ancestry would be important in further assessing this hypothesis.

Future Directions

Canada thistle can reproduce through seed to maintain high levels of genetic diversity, but its primary mode of reproduction is through vegetative growth, since seedlings are poor competitors, which rarely establish in existing pastures characterized by low-light conditions at early life stages (Moore, 1975; Donald, 1994; Heimann and Cussans, 1996; Leathwick and Bourdôt, 2012). Accordingly, future studies should concentrate on asexual processes (in particular shoot versus root elongation and resource partitioning/mobilization between shoot and root), to deepen our understanding of below-ground interactions. Since root production depends on the photosynthetic opportunity of shoots, thereby determining shoot population through over-wintering root biomass in the following season (Donald, 1994; Leathwick and Bourdôt, 2012), only multi-year experiments should further be undertaken, to adequately assess the long-term effects of stressors on the survival of this perennial plant.

Methods

Seed Collections

Seeds were collected from 22 populations in Europe (i.e., native range) and 20 populations in North America (i.e., introduced range) during summer 2008 (Figure 1 and Supplemental Table 1). Plant material from Europe was quarantined before being shipped to Canada under Canadian Food Inspection Agency (CFIA) import permit no. P-2008-03604.

Figure 1

Sampling Locations of Canada Thistle in the Introduced North American (Peach) and Native European (Turquoise) Ranges.

Common Garden Experiment

In July 2009, 20 seeds from five families (i.e., mother plants) from each population were scarified with 100-grit sandpaper, and sown in petri dishes, on damp filter paper soaked with 1% Plant Preservative Mixture (Plant Cell Technology, Washington, DC) to prevent fungal growth during germination. Petri dishes were then placed in a germination chamber set for 30°C, 80% humidity, and 16 h daylight. Radicle emergence was checked on a daily basis and germination success was recorded after 2 weeks. Newly-emerged seedlings were transferred onto a 1:1 sand/soil mixture when root growth reached 1–2 cm in length, and moved to a growth chamber set for 20°C –22°C, 40% humidity, and 16 h daylight. After 1–2 weeks’ acclimation in the growth chamber, when one to two pairs of true leaves were out, young plantlets were transferred into 9-cm pots filled with the same sand/soil mixture (except plants assigned to drought stress, which were transferred into 100% potting soil), and moved to a flood bench at the UBC Horticulture Greenhouse.

Six individuals from each family (i.e., 30 plants per population from five different mother plants) were randomly assigned to each of six growth conditions: (1) a control, non-stressful condition, (2) a nutrient stress, (3) a light stress, (4) a drought stress, (5) a simulated herbivory, and (6) a mowing treatment. All plants, with the exception of those assigned to the nutrient stress treatment, received 1.5 mL of Osmocote 13-13-13 slow-release fertilizer. The light stress was simulated by growing the plants in a shade box (approximately 4.5 × 3.5 × 1.8 m) made of PVC (polyvinyl chloride) pipes and covered by 121 Lee green filters (Andover, UK) and neutral density shade cloths, to mimic the spectral quality of light that is transmitted through the canopy of neighboring plants (Bonser and Aarssen, 2003). The green filter reduced light transmittance by 73%, and the shade cloth further reduced transmittance to 92% of the original value (light intensities were reduced from an average of 873.3 to 66.1 μmol·m−2·s−1 based on an average of three measurements taken at noon on sunny days). The drought stress was induced by discontinuing the watering of the plants by elevating the trays over the flood bench. We then measured the time until wilting or death. Herbivory was simulated by clipping off half of each new fully expanded leaf and by spraying 1 mM methyl jasmonate on each plant on a weekly basis. For the mowing experiment, above-ground biomass was cut off at 2 cm above ground on a monthly basis during a timeframe of 8 months, as long as any given plant was able to regenerate. Plants were randomly distributed within each treatment block on the bench, for a total of approximately 200 plants per treatment (20 populations times five individuals per population for each range). Within each block and treatment combination, individuals were randomly assigned to a tray, whose position on the flood bench was randomized every fortnight. All traits measured during the experiment are summarized in Supplemental Table 2.

Maternal Effects Experiment

Five plants from four native and four introduced populations (i.e., 20 genotypes per range), which were grown for 4 months in summer 2009 under controlled conditions, were cloned in fall 2009, to assess the extent of putative maternal effects. Briefly, clones were randomly assigned to two growth conditions (control, and nutrient stress as described above) and measured for the same traits summarized in Supplemental Table 2.

Statistical Analysis

Climate and Trait Principal Component Analysis

We conducted all statistical analyses in R v3.4.3 (R Core Team, 2017). We used the 19 WorldClim variables (Supplemental Table 3) and added a geographic dimension to the data by including longitude and latitude, with the latter potentially affecting season length and photoperiod. To explore associations between climatic variables in the sampled populations, we applied a principal component analysis using prcomp. We summarized, using the first two principal components, both trait (TRAITPC1 and TRAITPC2) and climate variation (CLIMPC1 and CLIMPC2) among populations using this approach.

Multivariate Analysis of Trait Differences

To explore patterns of climate-associated trait divergence as well as divergence among ranges, potentially indicative of local adaptation, we tested population mean trait responses to range (North America and Europe), climate variables, and their interaction in multi- and univariate models (MANOVA/MANCOVA or general linear models). We increased the power of the multivariate analysis (Scheiner, 2001) by removing highly correlated traits (Spearman's ρ > 0.7) and calculated the approximate F-statistics and Wilks' λ (multivariate F-value) to measure the strength of the associations. The reproductive traits flowering day and flower head number were only recorded in the control plot, as flowering was highly reduced in the other treatments. Therefore, we analyzed all traits including the reproductive traits from the control treatment only in a model with a main effect of range, and then a separate model including range as well as CLIMPC1 or CLIMPC2 and the interactions. As phenology is strongly affected by local growing season length in many flowering plant species, we also analyzed the control data using latitude as a covariate instead of CLIMPC1 or CLIMPC2. We subsequently removed those reproductive traits and analyzed all treatments (control, herbivory, nutrient), including the effects of treatment, range, and CLIMPC1 or CLIMPC2 as well as all interactions. We also conducted an additional analysis using initial leaf number prior to the onset of the treatments as a covariate to control for any size differences among the ranges that could confound our ability to address the impact of the stress on growth. We improved normality and reduced heteroscedasticity of the data by square root or log transforming traits where appropriate.

Univariate Analysis of Trait Differences

We conducted univariate tests of trait population means using general linear models, removing non-significant interactions in the same way as above. We calculated least-squares means and SEs with the package lsmeans and Tukey's honestly significant difference (HSD) post hoc tests. To explore differences in climate- or latitude-associated trait clines between ranges as revealed in the analyses of covariance (ANCOVAs), we tested for significant two-way interactions between range and the covariates, which is indicative of nonparallel trait-covariate slopes among the ranges. When the interactions with the covariates were significant, to further dissect the extent of trait divergence and its dependence on the covariates, we compared ANCOVA model estimates of traits in the introduced ranges with the native estimates at the minimum and maximum observed covariate values in each range, where the covariate values overlapped. We explored the highest-order significant interactions, using v2 tests with Holm p-value correction using the PHIA package.

We analyzed germination rate using a general linear model using an arcsine square root transformation of mean population germination rate and range as a main effect as well as with and without the covariates in the same manner as above. We did not include these data in the control MANOVA but analyzed them separately since germination rates of field-collected seed could be strongly influenced by the field environment.

Survival Analysis

We conducted a survival analysis for the drought, mowing, and shade experiment. Survival was assessed using Cox proportional hazards regression, using the coxme function in the coxme R package and the survminer package was used to visualize the results. Range was considered a fixed effect and population was included as a random term. A log likelihood ratio comparing the full and reduce model was used to test the significance of range.

Maternal Effects Analysis

For the maternal effects analysis, we implemented a similar approach to above. Population means were taken for each treatment and experiment (maternal effects or main experiment). The same maternal lines in the nutrient and control treatments from the main experiment were used for the maternal effects analysis. We then randomly removed a trait if it was highly correlated with another in pairwise comparisons (Spearman's ρ > 0.7). We conducted a MANOVA on the reduced trait set followed by univariate analysis to determine if the greenhouse produced clones and seed-grown maternal plants had similar trait patterns among ranges and treatments. We had three main effects in the model: treatment (control and nutrient), range (introduced and native) and experiment (main and maternal), and all interactions. We removed interactions that were non-significant starting with the highest-order interaction. Maternal effects that affected range and treatment patterns could be reflected in interaction of those effects with experiment. We calculated the approximate F-statistics and Wilks' λ (multivariate F-value) to measure the strength of the associations. For the same traits, we conducted a univariate analysis using general linear models, removing non-significant interactions in the same way as above. We calculated least-squares means and SEs with the package lsmeans and Tukey's HSD post hoc tests.

Funding

This work was supported by grants (PBZHP3-123301 and PA00P3_134180) from the to A.G., and from the Awards ( 327475 and 353026) to L.H.R.

Author Contributions

A.G., K.A.H., and L.H.R. designed the research. A.G. and K.N. conducted the experiment. K.A.H. and A.G. conducted the analysis. K.A.H., A.G., and L.H.R. wrote the manuscript.