Rapid increase in growth and productivity can aid invasions by a non-native tree.

Research on biological invasions has produced detailed theories describing range expansions of introduced populations. However, current knowledge of evolutionary factors associated with invasive range expansions, especially those related to rapid evolution of long-lived organisms, is still rudimentary. Here, we used a system of six 40-year-old invasive pine populations that originated from replicated introduction events to study evolution in productivity, growth, and chemical defence traits. We tested the hypotheses that invasive populations were undergoing rapid phenotypic change as populations spread, that populations exhibit trade-offs between evolution in growth and chemical defences, and that rates of rapid evolution in plant growth and productivity effect rates of invasion. Although all invasions started from replicated pools of genetic material and equal propagule pressure, we found divergence in mean values for the six invasive populations in the six traits measured. Not only were there between-population variations but also invasive populations were also rapidly changing along each invasive population expansion. Two populations displayed greater leaf areas (LAs) and smaller specific LAs (SLAs) during range expansion. Four populations had faster growth rates at the leading edge of the invasion front in comparison with plants at the rear edge. In terms of total plant defences, non-volatile resin increased in plants along one invasion gradient and decreased in a second, total needle phenolics increased in plants along one invasion gradient and total wood phenolics increased in plants along the one invasion gradient and decreased in a second. We found no trade-offs between investments in growth and chemical defence. Also, faster rates of change in growth rate and LA were positively associated with greater dispersal distances of invasive populations, suggesting rapid evolution may increase invasiveness. Understanding the roles of both natural and human-mediated ecological and evolutionary processes in population-level dynamics is key to understanding the ability of non-native species to invade. Published by Oxford University Press on behalf of the Annals of Botany Company.

Introduction

Biological invasions are a leading cause of environmental degradation, are a main focus of concern for conservation practitioners and provide important insights on species responses to climate change (Pyšek ; Caplat ; Moran and Alexander 2014; Kuebbing and Simberloff 2015). Invasive organisms alter ecosystem properties and community dynamics, impacting native species and ecosystem functioning (Wardle ; Yelenik and D'Antonio 2013). Driving the spread and impact of non-native populations are dynamic ecological and evolutionary processes acting at levels ranging from genes to global scale (Zenni 2014; Zenni ; Zenni and Hoban 2015b). Several decades of ecological research have produced detailed frameworks and theories describing range expansions of introduced populations (e.g. Blackburn ; Gurevitch ). More recently, researchers started to disentangle the evolutionary mechanisms driving spread of non-native populations and to incorporate them into invasion theory (Parker ; Prentis ; Sargent and Lodge 2014). However, the role of contemporary evolution in invasive range expansions is still poorly understood (Colautti and Barrett 2013; Zenni ), especially for long-lived organisms such as trees.

Climate has been hypothesized to be a major selection force on adaptation of invasive plant populations, changing plant phenology and size (Colautti and Barrett 2013; Colomer-Ventura ). However, evolution can also occur in growth rates and defence traits (Buswell ). Investments in defence against herbivory is costly for plants, and resource allocation for chemical defences can limit plant growth or reproduction (Coley ; Mithöfer and Boland 2012). The well-established relationship of constant trade-offs among investment in growth, reproduction and defence traits has led to the formulation of the ‘evolution of increased competitive ability’ hypothesis, which states that introduced plants are liberated from natural enemies and, thus, can allocate towards growth and reproduction resources previously required for defence (Blossey and Nötzold 1995). By increasing their resource investment in growth and reproduction, non-native populations may become abundant and widespread (Keane and Crawley 2002). Indeed, many invasives, including pines, escape from natural enemies when introduced to a new range (Liu and Stiling 2006). In their native ranges, several species of Pinus exhibit trade-offs between growth rates and chemical defences; slow-growing species and populations invest more in constitutive defences (Moreira ). For non-native populations, the absence of significant herbivory pressure may reduce plant resource allocation towards defence, possibly increasing fitness of non-native populations. Consequently, invasive populations may evolve increased growth rates and reproduction during the invasion process.

Besides local adaptation, genetic drift and phenotypic plasticity, the human-mediated introduction of pre-selected and adapted genetic lineages may also produce fit organisms in non-native ranges and benefit range expansions (Zenni ). The introduction of highly variable groups of individuals may produce the same effect (Forsman ; Zenni and Simberloff 2013; Forsman and Wennersten 2015). Further, the co-introduction of previously allopatric populations can lead to genetic recombination and novelty, possibly increasing levels of heterozygosity and polymorphism, which could trigger invasions that may not occur based on the original genotypes alone. However, the importance of admixture during invasions has not received much support from the literature (Rius and Darling 2014). While local adaptation, drift and phenotypic plasticity act during the naturalization stage and are part of all natural systems, the latter two mechanisms (introduction of pre-selected and adapted genetic lineages, and introduction of highly variable groups of individuals) act prior and during the introduction stage and may be exclusive to human mediated biological invasions (as opposed to natural range expansions, even those owned to human-mediated climate change). The processes acting at each stage of the invasion process can be different, and most studies have lumped together species at different stages of the introduction–naturalization–invasion continuum, or have looked only at processes occurring at the invasion stage (Blackburn ; Richardson and Pyšek 2012; Moodley ). Thus, understanding the roles of both natural and human-mediated agents prior to and during naturalization is key to understanding the invasion process. For this, researchers need to know the history of the introduction, including number and diversity of propagule sources, and follow the fate of the non-native populations during naturalization and invasion (Burton ; Donaldson ).

For animals, several studies have shown association between evolutionary dynamics and invasive range expansions. Invasive cane toads (Rhinella marina), for instance, evolved longer legs in the 70 years since introduction in Australia, increasing 5-fold the annual rate of expansion of the toad invasion front (Phillips ). Also, for invasive European starlings (Sturnus vulgaris) in South Africa, unfavourable environmental conditions enhanced dispersal, which preserves genetic diversity during range expansion and reduces potentially detrimental founder effects (Berthouly-Salazar ). For lizards (Podarcis muralis) invading in Germany, a study found a trend for an increase in genetic differentiation and a decrease in genetic diversity from the invasion centre to the expanding range edge, suggesting genetic drift as a major factor in the structuring of populations (Schulte ). For plants, virtually all examples of rapid evolution during invasive range expansions come from short-lived species. For example, in North America, local adaptation of purple loosestrife (Lythrum salicaria) in response to climate allowed populations to spread northward (Colautti and Barrett 2013). Further, the great majority of studies on evolution of invasive plants have compared invasive and native populations, or have compared distinct invasive populations. Very few studies have studied populations along invasion gradients.

In order to evaluate the role of human-mediated and natural evolutionary processes (natural selection and genetic drift) in the invasion success of non-native tree species, we designed a study using six fully replicated common garden experiments in southern Brazil where Pinus taeda L. (loblolly pine) was introduced at the same time (1973 and 1975), in the same numbers, from the same seed sources, and has formed invasive populations expanding outward from the plantations (i.e. invasive range expansions). Pinus taeda is a long-lived forest tree that has multi-generational populations, reproduces early (5 years) and yearly, and is wind-dispersed with viable seed dispersal distances of less than 20 m (Zenni ). The common garden experiments were originally planted to serve as forestry provenance trials for silvicultural purposes (Shimizu and Higa 1981). The experiments encompass a north-south transect covering about 850 km or 6° of latitude (Fig. 1). Previous work on the system showed provenance-by-environment interactions where genetic lineages of P. taeda exhibited differential naturalization success depending on the climate of the location into which they were introduced (Zenni ). Furthermore, 25 genes were undergoing significant shifts in allele frequencies along the invasion gradients (Zenni and Hoban 2015b). Although the genes evolving were mostly population specific, many were associated to important plant functions, such as phloem sugar transport, nitrate uptake, and pollen tube growth (Zenni and Hoban 2015b). Taken together, both studies make it evident that P. taeda is undergoing rapid evolution in all six invasive populations (Zenni ).

Figure 1.

Origin of the Pinus taeda seed sources in the USA and location of common gardens and invasive populations in Brazil. (A) Light grey represents the continental USA, dark grey represents P. taeda native range, and dots are the location of the 32 seed sources planted in the six common gardens. Open symbols represent the three genetic provenances. (B) Location of the six common gardens in Brazil. Solid symbols (circles, triangle, diamond and square) represent different climates. CB, Capão Bonito; IR, Irati; RN, Rio Negro; SFP, São Francisco de Paula; SM, Santa Maria; TB, Três Barras.

Here, we aimed to identify phenotypic changes undergone by the six populations of P. taeda in order to see how rapid evolution and genotype-by-environment interactions altered these invasive populations in growth, productivity, and chemical defence traits. We tested the following hypotheses: (i) the observed genotypic changes resulted in phenotypic changes for growth, productivity, and chemical defence traits; (ii) patterns of phenotypic changes match patterns of genetic change across populations; (iii) populations exhibit trade-offs between evolution in growth and chemical defences; and (iv) rates of rapid evolution in plant growth and productivity effect rates of invasion.

Methods

Study system

In 1973 and 1975, six common garden experiments for P. taeda were established in Brazil at the Santa Maria Experimental Farm (53.92°W 29.66°S), São Francisco de Paula National Forest (50.38°W 29.43°S), Três Barras National Forest (50.32°W 26.19°S), Rio Negro Experimental Station (49.76°W 26.05°S), Irati National Forest (50.57°W 25.36°S), and Capão Bonito National Forest (48.51°W 23.88°S). The six locations represent four climates and two of these four climates (Santa Maria and Capão Bonito) are different from the climate in P. taeda native range. Annual precipitations vary from 1212 mm to 2068 mm, whereas mean annual temperatures vary from 15.2 °C to 19.1 °C. Both precipitation and temperature patterns exhibit clinal variations from north to south. Climate types were determined using multivariate clustering analyses with data from Worldclim (see Zenni for details). Furthermore, Irati, Três Barras, Rio Negro and São Francisco de Paula original ecosystem are forest (Araucaria ombrophilous forest), Capão Bonito is a neotropical savanna (cerrado), and Santa Maria is a grassland ecosystem. Previous work showed no effect of original ecosystem type on the genetic changes observed for P. taeda (Zenni ; Zenni and Hoban 2015b).

The common gardens were planted with 29 or 32 seed sources of which 20 were present in all gardens (Zenni ). Each seed source corresponded to a seed lot collected from between 5 and 10 trees in natural stands in the species’ native range (Fig. 1) [see ]. In each Brazilian common garden, seed sources were planted in randomized blocks with four repetitions—a total of 144 trees from each seed source were planted in each common garden in four randomly placed squares of 6×6 trees (Shimizu and Higa 1981). Over the years, each common garden and its surroundings received circumstantial and haphazard management (e.g. some of the common gardens were thinned and plants growing on fire breaks and along roads were cut). The invasive populations themselves received only minor and dispersed management interventions. There were also high mortality rates for some seed sources in the common gardens probably owed to the high density of the plantations. In November and December 2014, all seed sources were still represented by at least 10 trees at any given garden, but the mean number of trees per seed source per site was usually higher. Previous work showed the seed sources form three native range genetic provenances (a provenance is an environmentally explicit genetic delineation of seed sources, and may be defined by a genetic cluster) and that the spatial distribution of these provenances across P. taeda native and invasive ranges correlates with temperature and precipitation patterns (Zenni ). The common gardens are considered parallel replicated introduction pools resulting in identical propagule pressures and residence times for these six locations (Zenni ; Zenni and Hoban 2015b).

Since introduction, the common gardens have produced invasive populations expanding between ca. 100 and 450 m from the common garden. By sampling plants at different distances from each of the common gardens, from the border of the plantations to the leading edge of the invasion fronts, we were able to track changes in genotype and phenotype frequencies in these six naturalized populations over multiple generations encompassing 40 years of population growth. This approach already showed that the six naturalized populations are undergoing contemporary evolution and significant genetic changes are occurring along the invasive range expansions. At the provenance level, genetic changes were associated with climate (Zenni ), whereas at the gene level, genetic changes were mostly population specific (Zenni and Hoban 2015b).

Data collection and trait measurements

The experiment design, description of sample and data collection, and specifications of the laboratory and genotyping methods used for the genetic work have been detailed previously (Zenni ; Zenni and Hoban 2015b). Briefly, 50 plants were haphazardly sampled from each of six invasive P. taeda populations (a total of 300 plants) and genotyped for 94 single nucleotide polymorphisms (SNPs) using Fluidigm® SNPtype Assays. Plants sampled were at least 1.3 m tall and had between 3 and 34 years of age. For this study, the original plantings were not included because traits from trees growing at high-density monocultures may not be comparable with trees naturally growing in heterogeneous habitats. All the SNPs chosen were a subset of the ones used by Eckert and were located in functional genes. Complete genotype data is available at Dryad data repository (Zenni and Hoban 2015, ). The distances between the plants and the edge of the common gardens were measured using the function “Hub distance” in the package MMQGIS for QGIS (Quantum GIS Development Team 2015), and distances were normalized for each population (divided by the maximum distance of spread) to account for variations in spread rates across sites (DNORM). For the phenotypic part of the study, fieldwork was carried out in November 2014 and January 2015. To collect phenotypic data, we visited the same trees previously sampled for genotyping in all six invasive populations. However, between 2012 (original sampling) and 2014 (sampling for this study), some trees died owing to idiosyncratic factors (i.e. tree falls and establishment of firebreaks). The number of surviving individuals was 47 for Capão Bonito, 48 for Irati, 49 for Rio Negro, 49 for Três Barras, 40 for São Francisco de Paula and 50 for Santa Maria, for a total of 283 plants. We used the same individuals from previous studies and did not include new plants because the sampled plants had already been genotyped.

For each plant, we collected an increment core at the base of the tree (∼ 20 cm from the ground) using a 30 cm increment borer (Haglöf Sweden), a handful of fully developed and healthy needles from the tip of the lowest branch, and a ∼15 cm segment of wood from the tip of the lowest branch. We also measured circumference at breast height (CBH; 1.30 m), tree total height, and bark thickness. All fresh wood and needle samples collected were stored in a freezer (−20 °C) immediately upon arrival at the laboratory. Increment cores were air dried, mounted in wooden supports and sandpapered to improve measurement accuracy. Lengths of growth rings were measured using a digital calliper to determine annual growth, and growth rings were counted to estimate the number of years contained in the extracted core (plant age). For all the analyses, we used the plants’ mean annual growth (MAG). We measured the area of 20 fully developed and healthy needles using a table scanner (Epson Perfection™ V700 Photo) and Digimizer v. 4.3.0 (MedCalc Software), and divided the total area by the number of needles to calculate mean leaf area (LA) of each plant. To determine specific LA (SLA), we dried the previously measured needles at 70 °C for 72 h and weighted them (leaf mass; LM) immediately after removing from the drying oven (SLA = LA×LM-1). To estimate inductive and constitutive plant defences in our studied plants, we determined non-volatile wood resin content, needle total phenolic content, and wood total phenolic content. Resin and phenolic extraction and determination followed the procedure described by Moreira . Briefly, resin was extracted from 15 cm long × ∼ 0.5 cm wide pieces of wood by mixing them with a solvent (hexane) in an ultrasonic bath at 45 °C for 20 min followed by a 24 h rest (performed twice to maximize extraction outputs). The solution was filtered in paper filters and left to dry. The residual content was weighted to the nearest 0.0001 g, and divided by the weight of the wood dry mass (oven dried at 80 °C for 24 h) to be expressed as mg g−1. Phenolic content was determined colorimetrically by the Folin–Ciocalteu method in a Bel Photonics 2000UV spectrophotometer (Bel Photonics do Brasil Ltda., Osasco, SP, Brazil) at 760 nm, using tannic acid as a standard (Moreira ). Needles were oven-dried for 72 h at 60 °C, ground in a Wiley mill, and sieved in a 5 mm mesh. For each sample, we used 100 mg of ground needles. The ground samples were transferred to centrifuge tubes, mixed with 5 ml of 70 % acetone and left in the refrigerator at 4 °C for 1 h. The solutions were centrifuged for 10 min at 7800 RPM and 0.3 ml of supernatant was removed and mixed with deionized water to make 1 ml of solution. To this solution, we added 0.1 N NaOH in Na2CO3 and the Folin–Ciocalteu reagent. After 2 h in the refrigerator, we took the spectrophotometer readings.

Statistical analyses

First, we performed six analyses of variance with Bonferroni correction of P values to test for differences in mean trait values of LA, SLA, MAG, non-volatile resin, and total needle and wood phenolic contents among the six locations. Second, we built generalized linear models with Gamma probability distribution and the inverse link function to test the effects of climate (mean annual temperature and annual precipitation) on each of the traits measured. Climate data were obtained from the Worldclim database (Hijmans ).

Next, we used linear mixed-effect models to test how LA, SLA, MAG, non-volatile resin, and total needle and wood phenolic contents varied as a factor of distance from the introduction point (DNORM) and/or were genetically controlled at each location. The genetic clustering results (provenances) from Zenni were used as the genetic makeup of individual plants. The dispersal model assumed for this study is that spread occurs as new generations establish farther from the point of introduction with some degree of back dispersal from the invasion front to the rear edge. Thus, the studied populations represent six unique invasions set out by six fully replicated introduction events. Distance of the invasive plant to the common garden and genetic clustering data were considered fixed effects, whereas plant age was added as a random effect (this was done to remove the effect of plant age from the model). For these analyses, we tested traits separately. For each trait, we built four models: a full model including distance, provenance and age; a distance model including distance as fixed and age as random factor; a provenance model including provenance as fixed and age as random factor; and a null model including only an intercept and age as random factor. We compared the four models for each trait using a likelihood ratio test and models were considered significant if statistically different from the null model at α = 0.05. Coefficients of determination were calculated using likelihood-ratio based pseudo-r2 and represent the variance explained by fixed effects (Nakagawa and Schielzeth 2013). These analyses were done in R 3.2 using the packages “lme4” v. 1.1-7 for mixed-effect models, and “MuMIn” v. 1.14.0 for pseudo-r2.

To test the association between mean annual growth rates and chemical defences, we built three linear models for each naturalized population using non-volatile wood resin content, needle total phenolic content or wood total phenolic content as a dependent variable and MAG as an independent variable.

We also used the slopes of the linear mixed-effect models to test if rate of change in the measured traits resulted in increased invasive potential. For this, we built linear models using the mixed-model slopes for each trait at each location as the independent variable and maximum distance of spread at each location as dependent variable. The slope coefficients are the mean slope for each trait of all plants in the model.

Finally, we performed a genome-wide association analysis to identify genes associated with the measured phenotypes. We tested each trait separately and assumed a codominant genetic model. The 94 unlinked SNPs were added as factors to the models and a Bonferroni correction was done to counteract the multiple comparisons problem. The analysis was done in R 3.2 using the package “SNPassoc” v. 1.9-2. SNPs identified as having associations with traits were further investigated for specific functions in the Dendrome database (http://dendrome.ucdavis.edu/DiversiTree/) and in Genbank.

Results

Although all invasions started from replicated pools of genetic material and equal propagule pressure (Zenni ), we found divergence in mean trait values in the six invasive populations for the six traits measured (Fig. 2 and Table 1): LA (F5,273 = 16.69, P < 0.001), SLA (F5,273 = 7.976, P < 0.001), mean annual growth (MAG; F5,276 = 15.56, P < 0.001), non-volatile resin content (F5,273 = 9.945, P < 0.001), needle total phenolic content (F5,273 = 19.36, P < 0.001) and wood total phenolic content (F5,264 = 16.63, P < 0.001). Climate (mean annual temperature and annual precipitation) explained divergence in LA (P < 0.02), SLA (F3,264 = 12.1, P < 0.001 for MAT), MAG (F3,264 = 22.5, P < 0.05), resin content (P < 0.02), needles phenolic content (P < 0.001) and wood phenolic content (P < 0.001 for AP) (Table 2 and Fig. 3). While LA and SLA were smaller in hotter and wetter locations, MAG was higher in locations with higher annual precipitations.

Figure 2.

Trait variations among the studied Pinus taeda invasive populations. Violin plots and boxplots of (A) leaf area (mm2), (B) specific leaf area (SLA, mm2/g), (C) mean annual growth (MAG, mm), (D) non-volatile wood resin content (mg/g), (E) total needle phenolic content (mg/g) and (F) total wood phenolic content (mg/g) for each invasive population (CB, Capão Bonito; IR, Irati; RN, Rio Negro; SF, São Francisco de Paula; SM, Santa Maria; TB, Três Barras). Violins (grey areas) are kernel probability densities of the data at different values. Boxplots are median (bold black line), quartiles (white rectangles), standard deviation (black lines) and possible outliers (dots).

Table 1.

Mean trait values and standard deviation for each invasive population. LA, leaf area; SLA, specific leaf area; MAG, mean annual growth.

	Capão Bonito	Irati	Rio Negro	Três Barras	São Francisco	Santa Maria
LA	210.96±38.64	264.82±61.45	211.45±44.06	257.65±63.79	191.42±39	206.39±40.73
SLA	7904.98±1812.52	10085.93±2657.68	9572.39±1728.83	9849.64±2351.04	9156.25±1682.43	8535.52±1726.33
MAG	7.43±4.18	12.05±6.92	9.2±3.45	13.23±5.39	13.48±6.08	15.27±4.07
Wood resin	43.4±21.79	49.09±26.59	37.93±16.31	52.57±11.12	33.78±10.14	33.16±10.82
Needle phenolics	9.47±5.38	4.75±2.8	10.06±4.56	8.05±3.55	13.07±6.48	13.78±7.07
Wood phenolics	9.2±4.06	9.25±3.44	13.55±4.68	11.33±3.87	15.49±5.52	15±5.33

Table 2.

Summary results of generalized linear models (Gamma family distribution) for the relationship between plant traits and climate. Null deviance shows how well the response variable is predicted by a model that includes only the intercept, whereas residual deviance shows how well the response variable is predicted by the alternative model. SLA, specific leaf area; MAG, mean annual growth; MAT, mean annual temperature; AP, annual precipitation.

Model	Variable	t	P	dfNull	dfResidual	DevianceNull	DevianceResidual
Leaf area	MAT	2.40	0.017	267	265	15.282	14.614
	AP	3.18	0.002
SLA	MAT	3.72	0.000	267	265	13.692	12.936
	AP	1.86	0.065
MAG	MAT	−2.03	0.043	267	265	97.28	89.144
	AP	−6.02	0.000
Wood resin	MAT	2.48	0.014	267	265	51.766	48.595
	AP	3.96	0.000
Needle phenolics	MAT	−3.73	0.000	267	265	98.674	88.972
	AP	−5.38	0.000
Wood phenolics	MAT	−1.05	0.295	267	265	49.618	44.821
	AP	−5.26	0.000

Figure 3.

Relationship between leaf area (mm2; A and B), specific leaf area (SLA, mm2/g; C and D), mean annual growth (MAG, mm/year; E and F), resin (mg/g; G and H), needle phenolic content (mg/g; I and J) and wood phenolic content (mg/g; K and L) and climate (mean annual temperature and annual precipitation) of the six studied locations. Solid lines represent local polynomial regression fitting and grey shades are 95% confidence intervals. Panels D and K show no relationship between variables.

Not only were there between-population variations but also the invasive populations were rapidly changing along each invasive range expansion (Fig. 4 and Table 3). We tested rapid evolution along the six invasive range expansions using linear mixed-effect models with distance from plantation and provenance as fixed effects and plant age as random effect. Two populations (Rio Negro and Três Barras) showed increases in LAs (Fig. 4 and Table 3, pseudo-r2 = 0.08 and 0.26, P = 0.045 and P < 0.001, respectively) and decreases in SLA (Fig. 4 and Table 3, pseudo-r2 = 0.1 and 0.2, P = 0.02 and 0.005, respectively) during range expansion. Four populations (Capão Bonito, Rio Negro, Três Barras and Santa Maria) showed faster growth rates at the leading edge of the invasion front in comparison with plants at the rear edge (Fig. 4 and Table 3, pseudo-r2 = 0.1, 0.1, 0.3, 0.1, and P = 0.03, 0.03, <0.001, 0.02, respectively). In two of these cases (Irati and Três Barras), MAG increased 2-fold in 40 years (Fig. 4). None of the populations showed decreases in MAG. In terms of constitutive and inductive plant defences, non-volatile resin content increased in plants along one invasion gradient (Rio Negro, pseudo-r2 = 0.1, P = 0.03) and decreased in a second (São Francisco, pseudo-r2 = 0.1, P = 0.04), total needle phenolic content increased in plants along one invasion gradient (Capão Bonito, pseudo-r2 = 0.1, P = 0.02) and total wood phenolic content increased in plants along the Capão Bonito invasion gradient (pseudo-r2 = 0.1, P = 0.03) and decreased in a second (Três Barras, pseudo-r2 = 0.1, P = 0.04).

Figure 4.

Relationship between productivity, growth, and defence traits and dispersal distance (D/Dmax, see Methods section for details). Scatterplots of leaf area (mm2), specific leaf area (SLA, mm2/g), mean annual growth (MAG, mm/year), non-volatile wood resin content (mg/g) and total needle phenolic content (mg/g) across the invasion gradient (normalized spread distance) for each invasive population. Lines represent local polynomial regression fit and grey shades are the 95% confidence intervals.

Table 3.

Summary results of likelihood ratio test among linear mixed-effect models. Pseudo-r2 represents the variance explained by fixed effects. Degrees of freedom (df) are reported for the likelihood ratio tests. SLA, specific leaf area; MAG, mean annual growth.

		Capão Bonito			Irati			Rio Negro			Três Barras			São Francisco de Paula			Santa Maria
		df	Pseudo-r2	P	df	Pseudo-r2	P	df	Pseudo-r2	P	df	Pseudo-r2	P	df	Pseudo-r2	P	df	Pseudo-r2	P
Leaf area	Full	5	0.012	0.551	5	0.017	0.811	5	0.179	0.016	5	0.283	< 0.001	5	0.050	0.198	5	0.024	0.447
	Invasion	4	0.006	0.615	4	0.003	0.702	4	0.080	0.045	4	0.262	< 0.001	4	0.049	0.196	4	0.016	0.368
	Provenance	4	0.003	1.000	4	0.016	< 0.001	4	0.074	1.000	4	0.005	1.000	4	0.002	1.000	4	0.012	1.000
SLA	Full	5	0.084	0.108	5	0.110	0.268	5	0.103	0.024	5	0.303	0.003	5	0.056	0.232	5	0.004	0.965
	Invasion	4	0.040	0.169	4	0.030	0.184	4	0.103	0.023	4	0.165	0.005	4	0.039	0.235	4	0.000	0.969
	Provenance	4	0.032	1.000	4	0.088	< 0.001	4	0.002	1.000	4	0.136	1.000	4	0.017	1.000	4	0.004	< 0.001
MAG	Full	5	0.191	0.010	5	0.158	0.113	5	0.098	0.032	5	0.338	< 0.001	5	0.049	0.205	5	0.092	0.026
	Invasion	4	0.096	0.032	4	0.082	0.070	4	0.098	0.030	4	0.335	< 0.001	4	0.049	0.199	4	0.092	0.017
	Provenance	4	0.068	1.000	4	0.093	< 0.001	4	0.002	1.000	4	0.002	1.000	4	0.002	1.000	4	0.016	1.000
Wood resin	Full	5	0.020	0.473	5	0.034	0.211	5	0.135	0.014	5	0.117	0.179	5	0.146	0.037	5	0.004	0.971
	Invasion	4	0.007	0.557	4	0.034	0.203	4	0.100	0.025	4	0.041	0.166	4	0.114	0.040	4	0.000	0.896
	Provenance	4	0.009	< 0.001	4	0.001	1.000	4	0.020	1.000	4	0.081	< 0.001	4	0.033	1.000	4	0.004	< 0.001
Leaf phenolics	Full	5	0.110	0.021	5	0.021	0.527	5	0.008	0.549	5	0.001	0.850	5	0.034	0.760	5	0.036	0.694
	Invasion	4	0.109	0.021	4	0.007	0.581	4	0.007	0.562	4	0.001	0.844	4	0.002	0.770	4	0.008	0.521
	Provenance	4	0.000	1.000	4	0.015	< 0.001	4	0.000	1.000	4	0.000	1.000	4	0.032	< 0.001	4	0.033	< 0.001
Wood phenolics	Full	5	0.116	0.048	5	0.041	0.243	5	0.010	0.663	5	0.147	0.012	5	0.023	0.468	5	0.005	0.616
	Invasion	4	0.094	0.034	4	0.040	0.257	4	0.005	0.634	4	0.131	0.012	4	0.019	0.453	4	0.005	0.615
	Provenance	4	0.036	1.000	4	0.005	1.000	4	0.006	< 0.001	4	0.019	1.000	4	0.008	1.000	4	0.000	1.000

Some instances revealed that evolutionary forces did not produce phenotypic changes along the invasive range expansions for the traits measured. Some trait variations were explained by the genetic provenance of the plants (i.e. genetically conserved) (Table 3). In Irati, LA (pseudo-r2 = 0.02, P < 0.001), SLA (pseudo-r2 = 0.09, P < 0.001), MAG (pseudo-r2 = 0.09, P < 0.001) and total phenolic (pseudo-r2 = 0.02, P < 0.001) values were explained by the plants genetic provenance. Ancestry also explained SLA in Santa Maria (pseudo-r2 = 0.004, P < 0.001), wood resin in Capão Bonito, Três Barras and Santa Maria (pseudo-r2 = 0.009, pseudo-r2 = 0.08, pseudo-r2 = 0.004, all P < 0.001), needle total phenolic content in São Francisco and Santa Maria (pseudo-r2 = 0.02 and 0.03, and P < 0.001 for both) and wood total phenolic content in Rio Negro (pseudo-r2 = 0.006, P < 0.001).

Our analyses of potential trade-offs between growth and chemical defences yielded no significant relationships (P > 0.05) among MAG and wood resin, wood total phenolics, or needle total phenolics in all but two models (Fig. 5). In Irati, concentrations of total phenolics in needles were lower in plants with higher MAG (r2 = 0.12, P = 0.02), but a single outlier (IRT6) with very high phenolic content drove this result. Removing the outlier from the model yielded a non-significant relationship (P = 0.09). In Capão Bonito, contrary to all expectations, concentrations of total phenolics in needles and MAG were positively associated (r2 = 0.22, P < 0.001).

Figure 5.

Relationship between chemical defences and mean annual growth in plants across the six studied naturalized populations. (A) Relationships between resin contents and mean annual growth, (B) relationships between needle total phenolic contents and mean annual growth and (C) relationships between wood total phenolic contents and mean annual growth. Dots of different shapes represent different naturalized populations. None but one of the relationships were statistically significant (needle phenolics and MAG in Capão Bonito; see text for details).

When testing how rate of evolutionary responses affected invasive ranges expansions, we found that rate of change of increased MAG and rate of change of increased LA were positively associated to the total spread distances of the invasive populations (r2 = 0.7 and P = 0.04 for both models, Fig. 6). However, although plant MAG and LA were not correlated (P = 0.07) and both were positively related to total spread distance, the full model was not significant (P = 0.2). This was probably due to the small sample size (n = 6 populations).

Figure 6.

Relationship between evolution in leaf area and growth and invasive range expansion. Regression line and 95% confidence intervals (grey shade) of the relationship between the rate of evolution in (A) leaf area and (B) mean annual growth (MAG) for the six invasive populations. Evolution is represented by the mean slope of the mixed models between normalized distance of spread and trait value (Fig. 4).

In a genome-wide association analysis, we detected one gene associated with LA (SNP_210094-Pita) and two genes associated with MAG (SNP_215552-Pita and SNP_225711-Pita). SNP_210094-Pita is related to phosphorus assimilation by the roots (GenBank Acc: BE496394) and phosphorous has been experimentally shown to affect LA in forest trees (Herbert and Fownes 1995). SNP_215552-Pita is related to root development (GenBank Acc: CO198215), drought resistance (GenBank Acc: CF392949) and fungal resistance (GenBank Acc: DR093744). SNP_225711-Pita is associated with wood formation (GenBank Acc: BE496394).

Discussion

Our results show strong evidences for rapid phenotypic change in all six invasive populations, providing support for our first hypothesis. However, contrary to our second hypothesis, patterns of phenotypic changes did not match patterns of genetic change across populations, and each population showed a unique pattern of phenotypic change along the invasion gradient. Surprisingly, also contrary to our initial expectations, populations did not exhibit trade-offs between evolution in growth and chemical defences. In one instance (Capão Bonito), there was a positive association between MAG and leaf phenolic content. Finally, supporting our last hypothesis, we found positive associations between rates of rapid change in plant growth (MAG) and productivity (LA) and the rates of population spread. The mismatches observed between phenotypic changes, genetic makeup of the plants and climate patterns suggest part of the variation found may be caused by phenotypic plasticity or genetic drift. Although populations are changing genetically along the invasion gradients (Zenni ), the traits measured might be responding to the environment without a underlying genetic change. Also, the genetic and phenotypic changes might not be the result of selection, but the result of genetic drift and founder effects (van Kleunen and Fischer 2008; Monty ).

A previous study in this system showed that climate was one of the factors determining the patterns of genetic change across the six populations (Zenni ). Accordingly, climate (mean annual temperature and annual precipitation) also partially explained patterns of variation in LA, SLA, MAG, wood resin and leaf and wood phenolic contents among populations (Fig. 3). The strongest association was between MAG and annual precipitation, where growth rates were higher in wetter locations. Although statistically significant, mean annual temperature did not strongly affect the growth rates of plants. Relationships between P. taeda phenotypic patterns and climatic factors had also been reported in the species native range (e.g. Eckert ). For other systems, climate has also been shown to be an important factor effecting population growth, spread, or retraction (Chen ; Pfeifer-Meister ; Kerr ). In our study, the fact that all populations started from replicated introduction pools, and genotypic and phenotypic responses to each climate were observed during range expansions in less than 40 years, suggests climate may be an important driver of evolution for invasive populations. Under new climates, either owing to human-meditated introductions to a new range or to climate change, invasive populations may be the ones able to adapt quickly to the new environment.

Climate not only can drive ecological patterns of plant chemical defences but also the evolution of constitutive chemical defences in pines may be driven by it (Moreira ). Our results seem to support this notion (Fig. 3), although we measured total phenolic content, which includes both constitutive and inducible defences. In the current study, we found differences in chemical defence traits both among and within populations (Table 2 and Figs 2 and 4). The among population differences may be due to ecological responses of populations to each environment (climate and habitat), whereas the changes along each invasion gradient may be showing evolutionary responses of each population to each environment. Considering that the among population differences are stronger than the within population changes, it is not possible to determine if adaptive evolution is occurring in chemical defences in the six studied populations or if differences are due to genetic drift caused by founder effects (Monty ). However, it is possible evolution in the studied defence traits is occurring, but the effects were not detected. Another study, disentangling constitutive and inducible defences would be necessary to address the genotypic (constitutive), phenotypic (constitutive and inducible) or epigenetic components of defence responses.

The production of chemical defences is costly, and increased resource allocation to defences tends to result in disproportionally high decreases in growth rates for plants (Mithöfer and Boland 2012). Consequently, if non-native plants escape from natural enemies during introduction and are able to afford neglecting their chemical defences to favour growth, reproduction, and dispersal traits, they may become invasive (Callaway and Ridenour 2004; Felker-Quinn ). However, in our study, we generally found no relationship between investments in plant growth, productivity and chemical defences. Plants showed increased or kept existing resource use strategies, high growth rates and levels of chemical defences. The quantities of resin and total phenolics we found were similar to quantities reported for the species by other studies in other regions (Moreira ), indicating that P. taeda plants in the six studied populations had levels of chemical defences similar to the species' average. Although investments in chemical defences varied across populations, we did not find a consistent trend towards increased or decreased investments along the invasion gradients as we did for MAG. On one hand, resin and phenolics total contents might be a plastic response of populations and, given the lack of important natural enemies in the introduced ranges, there are no selection pressures for increased chemical defences. Or populations show a slower pace of change for chemical defences than for growth rates given the complexity of the metabolic processes involved in producing chemical defences (Mithöfer and Boland 2012). On the other hand, there may be no limiting resources requiring a trade-off between chemical defences and growth and increase growth rates are simply a result of favourable climatic conditions for P. taeda (e.g. Nuñez and Medley 2011).

Although the invasion literature advocates that high SLAs and fast growth have positive associations with woody plant invasiveness (Grotkopp and Rejmánek 2007; Zenni and Simberloff 2013), it has also been suggested that natural selection can reduce a population invasive potential (Lankau ). Our data suggest that rapid evolution, either via natural selection or genetic drift, can result in populations with increased invasive potential in comparison to the original introduction pool. In most of our cases, plants at the leading edges of the invasion fronts grew faster and were more productive, on an average, than plants at the rear edges (Fig. 4). These effects were mediated by climate, suggesting that in conditions adverse for the species, evolutionary changes may not be so pronounced. However, in none of the cases did populations showed slower growth rates or lower productivities at the leading edges. We found two populations (Rio Negro and Três Barras) that produced lower SLA along the invasion gradients, besides being actively spreading and invading. Four populations increased growth rates along the invasion gradients. Although our results seem to contradict the established literature (Tecco ; Porté ), we believe that the decrease in SLA combined with increases in LA may be explained as a response of P. taeda plants to the forest ecosystems of Rio Negro and Três Barras. MAG results were in agreement with the established literature that fast growth is a key invasion trait (Lamarque ).

Also supporting the idea that rapid evolutionary change in growth rates and productivity can increase the invasive potential of populations (Lamarque ; Monty ; Sargent and Lodge 2014), we found strong relationships between MAG and LA and total spread of populations. Populations with greater rate of change towards increased growth rates and LA spread farther in the 40 years period encompassed by this study (Fig. 6). However, it is also possible that greater spread could lead to higher growth if plants are experiencing less competition at the leading edge of the invasion front, or if growth and spread were both correlated to a third factor. Nevertheless, the strong effect rapid evolution can have on populations’ invasiveness suggests that invasions must be dealt with sooner rather than latter.

Several of the statistical models including only genetic provenance as a factor were significantly different than the null models, but overall they had low explanatory power (Table 2). It is possible that genetic provenance is a coarse scale to test fine evolutionary changes in populations. The fine-scaled genome-wide association analyses, where trait values were tested for each one of the 94 genes included in this study, we found three genes associated with the traits measured. Given that we used putatively functional markers, we expected that more of the genes would be associated with traits. A possible explanation for the few significant associations found is that most traits result from the expression of many genes, instead of only a few genes of large effect. Also, because we measured plants growing in natural environments and not in a controlled homogenous setting (e.g. greenhouse), phenotypic plasticity may be masking the role of specific genes in trait expression. Taken together with previous findings on the study system that several genes were undergoing rapid changes in genotype frequencies (Zenni and Hoban 2015b), the results support our claim that evolution is occurring in this system and changes were not caused only by phenotypic plasticity.

Because our study took place in a natural setting, and not a controlled environment, it is subject to several caveats. It is possible some of the variations in traits values could be explained by the age and size of the trees, as well as by micro-environmental factors that differ for each individual tree, and that may vary along each of the invasion gradients. We tried to address the effect of tree size and age in the statistical models by having tree age as a random effect in the mixed models and by using only the variance explained by the fixed effects in subsequent analyses and interpretations of the model results. We were unable to disentangle the relative importance of adaptive evolution, genetic drift and phenotypic plasticity for the traits measured. However, given the structure of our statistical models, it is likely that plastic responses only decreased the chances of us finding significant relationships between trait evolution and invasion spread. We still found a number of significant associations. Further, although we acknowledge the existence and the potential importance of micro-environmental factors in affecting trait values on each plant (adaptive and plastic), we think they make our findings even more remarkable. We showed that despite all potential micro-environmental variations in each of six naturalized population spanning an 840 km transect, four climate zones, and four decades of natural processes, rapid evolution in several functional traits could still be detected along invasion gradients.

Conclusions

In summary, our study provides strong support for the role of rapid evolution on the success of non-native range expansions of an invasive species. We also show that evolutionary changes are in part associated with climate, but are mostly fine-scaled and context-specific, even though invasive populations started from fully replicated introduction events. The results of this study highlight the unique nature of the ecological and evolutionary dynamics of each population, suggesting that predicting invasion trajectories may be daunting. Further, we demonstrate that the potential for rapid evolution can affect the capacity of populations to cope with new environments. Taken together, the capacity for adaptation to different conditions, evolution of increased growth rates, and the lack of trade-offs between growth and defences make P. taeda an invasive species that requires management and control before its spread reach large areas.

Sources of Funding

Our work was funded by CNPq-Brazil (313926/2014-0).

Contributions by the Authors

R.D.Z. conceived the study, collected data, measured resin contents, did statistical analyses and wrote the manuscript. W.L.C. collected data, measured SLAs and mean annual growths. G.S. measured phenolic contents. All authors discussed the results and commented on the manuscript.

Conflict of Interest Statement

None declared.