Evolving while invading: rapid adaptive evolution in juvenile development time for a biological control organism colonizing a high-elevation environment.

We report evidence of adaptive evolution in juvenile development time on a decadal timescale for the cinnabar moth Tyria jacobaeae (Lepidoptera: Arctiidae) colonizing new habitats and hosts from the Willamette Valley to the Coast Range and Cascades Mountains in Oregon. Four lines of evidence reveal shorter egg to pupa juvenile development times evolved in the mountains, where cooler temperatures shorten the growing season: (i) field observations showed that the mountain populations have shorter phenological development; (ii) a common garden experiment revealed genetic determination of phenotypic differences in juvenile development time between Willamette Valley and mountain populations correlated with the growing season; (iii) a laboratory experiment rearing offspring from parental crosses within and between Willamette Valley and Cascades populations demonstrated polygenic inheritance, high heritability, and genetic determination of phenotypic differences in development times; and (iv) statistical tests that exclude random processes (founder effect, genetic drift) in favor of natural selection as explanations for observed differences in phenology. These results support the hypothesis that rapid adaptation to the cooler mountain climate occurred in populations established from populations in the warmer valley climate. Our findings should motivate regulators to require evaluation of evolutionary potential of candidate biological control organisms prior to release.

Introduction

Population ecology has played an important role as a guide and explanation for biological control, but the importance of population genetics and adaptive evolution in biological control remains a controversial topic (Force 1967; Remington 2008; Messenger and van den Bosch 1971; Roush 1990; Hopper et al. 1993; Holt and Hochberg 2011; Jervis 1997; Hufbauer 2002; Hufbauer and Roderick 2005; Phillips et al. 2008; Henry et al. 2010). Arguments for and against evolution as a driver of biological control outcomes mirror those arising at the boundaries of ecology and evolution in other contexts, including climate change (Travis and Futuyma 1993; Hoffmann and Sgrò 2011). Contemporary evolution may play a minor role in the outcome of biological control if (i) local extinction is more common than local adaptation as a response to environmental change (Parmesan 2006), (ii) most organisms lack sufficient genetic variation to adapt to changes in abiotic and biotic environments (Parmesan 2006), and (iii) the fundamental niches of introduced biological control agents are conserved, predictable, with little evidence of contemporary adaptive evolution in response to climate (Zalucki and van Klinken 2006) or hosts (Pemberton 2000; Schaffner 2001; van Klinken and Edwards 2002; Zalucki and van Klinken 2006; Thrall et al. 2011; Petitpierre et al. 2012). Conversely, contemporary evolution may play a major role that until recently has gone largely undetected in biological control. Rapid evolutionary change is a well-documented response to rapid environmental change in other contexts (Carroll and Fox 2008), and ecological and evolutionary processes can act and interact on similar timescales (Hairston et al. 2005; Schoener 2011). Populations and species are not genetically or ecologically uniform entities (Burdon et al. 1981; Wajnberg 2004), and evolution (i.e., a change in gene frequencies) can occur in control-organism populations from standing genetic variation (Barrett and Schluter 2007) without mutation because of forces of genetic drift, founder effects, and natural selection (Hufbauer and Roderick 2005).

There have been few reports of adaptive evolution in biological control organisms colonizing new areas despite the ∼100 years over which releases have been made. Changes in gene frequencies in control organisms (Phillips et al. 2008) or their hosts (Burdon et al. 1981) have been reported in a few cases, and a few investigators have screened biological control organisms for genetic variation on which natural selection might act, as in the case of parasitoids (Wajnberg 2004). Rapid life-history evolution has been documented for the Italian biotype of the ragwort flea beetle Longitarsus jacobaeae (Coleoptera: Chrysomelidae), introduced to western Oregon from Italy via California to control ragwort Jacobaea vulgaris (Asteraceae). A Cascade Mountain population derived from a Willamette Valley population evolved shorter duration of summer diapause (time between emergence and first oviposition), yielding a better match between insect phenology and the shorter summer season in the mountain environment (Szűcs et al. 2012). A leaf beetle Diorhabda carinulata (Coloptera: Chrysomelidae) introduced into North America from China for control of an exotic shrub Tamarix spp. evolved shorter critical day length for diapause induction (day length at which 50% of the population enters diapause) in populations located south of the latitude of origin, producing a closer local match between insect and plant phenologies (Bean et al. 2012). The general absence of evidence for evolutionary change is not, in fact, evidence of absence. There is a need for more eco-evolutionary studies to assess whether and how control organisms colonizing new environments adapt to new climates and new hosts (Secord and Kareiva 1996; Simberloff and Stiling 1996; Hopper 2001; Myers 2001; Louda et al. 2003; Cox 2004; Sheppard et al. 2003) and to determine whether rates of evolution can be considered rapid because of their potential to influence an ongoing ecological process like the invasion of new habitats and the acquisition of new hosts (Hairston et al. 2005; Schoener 2011).

Here, we report evidence of rapid adaptive evolution in phenology and juvenile development times for the cinnabar moth Tyria jacobaeae (L.) (Lepidoptera: Arctiidae) derived from the Willamette Valley (87 m) and colonizing new habitats and hosts in the high elevations of the Coast Range (877 m) and Cascade Mountains (1572 m) of Oregon, USA. Our study of the cinnabar moth and a similar study also conducted in Western Oregon of the ragwort flea beetle L. jacobaeae (Szűcs et al. 2012), present a unique opportunity to test whether unrelated organisms introduced for biological control of the same host in the same region show similar patterns of evolution. We used four lines of evidence to test whether shorter juvenile development time evolved after anthropogenic redistribution of the cinnabar moth from low elevations of the Willamette Valley to the high elevations of the Coast Range and Cascade Mountains, where temperatures are cooler and growing seasons are shorter. We compare: (i) field observations on phenology of low- and high-elevation populations, (ii) development times of field-collected individuals reared in a common garden experiment in a greenhouse, (iii) development times of laboratory-reared offspring from parental crosses within and between populations, and (iv) observed differences in phenology with those expected under two null hypotheses of evolution by random processes (founder effects, genetic drift) and the alternative hypothesis of evolution by natural selection. Results show that biological control organisms are capable of rapid adaptive evolution while colonizing and invading new environments. The practical consequence is that evaluating the evolutionary potential of candidate agents may become a necessary component of risk-benefit-cost analysis prior to their release. We conclude by offering some simple ways to assess evolutionary potential of biological control organisms prior to release in novel environments.

Study system

Our study system consists of an interaction between a phytophagous insect species and two of its host plant species in western Oregon: the cinnabar moth (T. jacobaeae) native to Europe, its ancestral host tansy ragwort (Jacobaea vulgaris Gaertn = Senecio jacobaea L., Asteraceae) native to Europe, and one of its acquired host arrow leaf ragwort (Senecio triangularis Hook.) native to North America.

The cinnabar moth was initially introduced from Europe (near Paris, France) to North America (Ft. Bragg, California) in 1959 and then to Oregon (in Scio, Linn County, Oregon on the east side of the Willamette Valley near the Cascade Range and in Valley Junction, Polk County, Oregon in the foothills of the Coast Range) in 1960 (Frick and Holloway 2000; Ritcher 1966). The insect has expanded its geographic range from the Willamette Valley to the Cascades and Coast Range Mountains and its host range from its host (J. vulgaris) native to Europe and targeted for biological control abroad in North America, to new hosts (S. triangularis and Packera pseudaurea (Rydb.) W. A. Weber & Á. Löve var. pseudaurea, Asteraceae) native to North America. Human-assisted (i.e., anthropogenic) redistribution of the T. jacobaeae for ragwort biological control over the period 1970–1985 exposed new habitats and hosts in the Coast Range and Cascade Mountains, where outbreaks of T. jacobaeae occurred from 1985 to 1989 (Coombs et al. 1996).

The climatic tolerances of T. jacobaeae appear to be broad, as inferred from its wide geographic distribution of natural populations and broad thermal tolerance (10–35°C) exhibited in laboratory populations under controlled conditions (Harman et al. 1990). However, local genetic variation and adaptation in climate tolerances have not been previously reported. Genetic differentiation of T. jacobaeae populations introduced into western North America has been inferred from variation in isozyme frequencies (Myers 1978). Heritable variation has been reported for adult emergence times (Richards and Myers 1980), but previous studies report populations of T. jacobaeae living in different climates at low elevations in western North America have similar minimum threshold temperatures and development times for egg to pupal stages (Myers 1979).

The timing of events in the univoltine life cycle of T. jacobaeae in relation to its seasonal environment is as follows (Dempster 1982). Overwintering is in the pupal stage, adults emerge in spring, mate, and adult females oviposit in batches (mean of 40 eggs per batch) on the underside of basal leaves of the host. Larvae develop through five stages before entering the pupal stage. Dispersal occurs in both the adult and late-larval stages, but dispersal distances are generally short (maximum of 300 m for adults)(Harrison et al. 1995; Rudd and McEvoy 1996), and without human assistance, spatial spread rates (combining population growth and dispersal) are slow (41 m year−1)(Hawkes 1967).

The fundamental host range of T. jacobaeae has been estimated by host specificity tests conducted prior to introduction in New Zealand (Miller 1929; Cameron 1935), the United States (Parker 1960), and Canada (Bucher and Harris 1961). The fundamental (=genetically determined) host range appears to be broad, including Jacobaea, Senecio, Packera, Erechtites, and Petasites species; it is constrained by phylogeny, chemistry, and perhaps plant architecture and phenology (Wink and Legal 2001; Bernays et al. 2004). The realized host range is much narrower than the fundamental host range; reports of host use in Oregon environments include natives of Europe (J. vulgaris, S. vulgaris, S. sylvaticus) and North America (S. triangularis, Packera pseudaurea var. pseudaurea) (Diehl and McEvoy 1990; Harris and McEvoy 1995; McEvoy et al. 2008). Cinnabar moth populations introduced in other countries reportedly use plants native to those countries including S. minimus and S. biserratus in New Zealand (Fowler et al. 2000; Paynter et al. 2004).

Our study was carried out along an environmental gradient marked by regions and sample locations labeled Coast, Coast Range, Willamette Valley, and Cascades Mountains. The environmental gradient varies in abiotic characteristics (elevation, distance from the coast, temperature, and precipitation) and biotic interactions (with host plants) (Table 1). Prior studies confirm strong biotic interactions between T. jacobaeae and its ancestral host J. vulgaris (McEvoy et al. 1991) and its acquired host plant species S. triangularis (Diehl and McEvoy 1990; Harris and McEvoy 1995; Fuller 2002) in Oregon measured by high levels of defoliation and defloration in local populations. Ragwort abundance has declined 97% from former levels in Oregon following the introduction of three insects – cinnabar moth T. jacobaeae, ragwort flea beetle Longitarsus jacobaeae (Waterhouse) (Coleoptera: Chrysomelidae), and ragwort seed head fly Botanophila seneciella (Meade) (Diptera: Anthomyiidae) – for biological control (McEvoy et al. 1991). However, the wide distribution and high abundance of the nontarget species, S. triangularis, appears to remain unchanged despite occasionally high levels of repeated defoliation (K. M. Higgs and P. B. McEvoy, unpublished data).

Table 1

Regional abiotic and biotic characteristics for four locations along an environmental gradient in Western Oregon, USA. We report comparisons of Coast, Willamette Valley, and Cascades; comparisons of the Coast Range and Willamette Valley are reported in (Murtaugh et al. 1987)

Region	Coast	Coast Range	Willamette Valley	Cascades
Elevation (m)	31	877	87	1572
Distance from Coast (km)	4	33	60	190
Ave Annual Temp (°C)	10.7	8.13	11.4	4.7
Annual Precip (cm)	262	377	122	226
Host Plant	Jacobaea vulgaris	Senecio triangularis	Jacobaea vulgaris	Senecio triangularis

Methods

Phenology

We compared phenologies of populations and the match between phenology and the length of the growing season using the following conventions. We used physiological time (accumulated degree days above a minimum threshold temperature) as a timescale. We defined phenology as the distribution of life-cycle stages in physiological time, the speed of phenological development as the rate of progress through the stages, and the length of the growing season in physiological time. We used two population statistics to characterize phenology, mean time in stage for each stage and mean stage for each sampling occasion. To illustrate using a hypothetical example, we assume the data for a particular population in Table 2. The sets of observed times in stage for each stage would be:

Table 2

Hypothetical data used to illustrate estimation of two population statistics that characterize phenology, mean time in stage for each stage and mean stage for each sampling occasion

		Number observed
Sample	Time	Stage 1	Stage 2
1	2.6	3	0
2	4.7	5	1
3	4.8	2	4
4	5.1	1	7

X1 = (2.6, 2.6, 2.6, 4.7, 4.7, 4.7, 4.7, 4.7, 4.8, 4.8, 5.1)

X2 = (4.7, 4.8, 4.8, 4.8, 4.8, 5.1, 5.1, 5.1, 5.1, 5.1, 5.1, 5.1)

Mean times in stage for each stage would be = 4.18 and = 4.97.

The sets of stages for each observation time would be:

X2.6 = (1,1,1), X4.7 = (1, 1, 1, 1, 1, 2), X4.8 = (1, 1, 2, 2, 2, 2), X5.1 = (1, 2, 2, 2, 2, 2, 2, 2)

Mean stage for each time would be: = 1.00, = 1.14, = 1.67, = 1.88.

We estimated the minimum threshold temperature using data from Harman et al. (1990). We regressed the inverse of the mean development time in days (=the development rate) on temperature, assuming a linear model, and estimated the slope and X-axis intercept. We used mean development times for each temperature instead of individual development times to ensure that each temperature is weighted equally in the regression. The regression equation for mean speed of development (1/development time) (Y) in relation to temperature (X) is Y = 0.0224 + 0.0026 X, R2 = 0.995, P = 0.00015. The inverse of the slope estimates the degree days required for development from egg to pupa (the thermal constant K = 384.6 degree days), and the X-intercept (8.6°C) estimates the minimum threshold temperature for development. We approximated the environmental temperature wave for each geographic location using 30-year means from the OSU PRISM Group website (http://prism.oregonstate.edu) interpolated for latitude and longitude coordinates of each site.

In 2010, we made weekly observations on the life-stage phenologies (adult, egg through 5th instar) at one Willamette Valley and one Cascades population. Sites were visited 18 and 11 times, respectively. At weekly intervals, we randomly selected 30 plants from a sample of 50 marked plants within the local population and thoroughly searched each plant for eggs and larvae. All cinnabar moth life stages were counted and recorded.

In 2011, we sampled two replicate populations within each region (Baskett Slough and Alsea Highway in the Willamette Valley and Pacific Crest Trail and Wasco Lake in the Cascades) at three times within the year at approximately 300, 400, and 500 degree days. We selected a random sample of plants as in 2010, except we made repeated observations on the same plants on three sampling occasions per population.

We compared phenology between the Cascades and the Willamette Valley using slightly different methods dictated by differences in survey design between years. In 2010, we compared mean times in stage (in degree days) for each stage for complete phenology curves estimated at one Cascades and one Willamette valley population using methods described by Murtaugh et al. (1987), with a null hypothesis of equal mean times and a one-sided alternative hypothesis that mean times in stage are greater (i.e., later) for the Willamette Valley population than for the Cascades population. Transformations did not improve the model fit, and residuals demonstrated some unequal variance owing to the widely varying times in stage when degree days range from 200 to 300. In 2011, we compared mean stage on three sampling occasions (yielding snap shots of phenology) for two replicate populations within each region, Willamette Valley and Cascades. Coding the six stages (Egg, L1, L2, L3, L4, L5) from 1 to 6, the mean stage was calculated across all individuals observed on each plant. Mean stage serves as an index of phenology of the population over all stages and was modeled as a function of accumulated degree days for comparisons of Cascades and Willamette Valley regions. Multiple linear regression was used to estimate fixed effects for regions, populations within regions, an overall slope, and region- and population-specific slopes. These fixed effects were used to calculate a mean slope for each region, and the difference between region-level slopes was tested with a one-sided t-test for an alternative that the mean stage of Cascades populations is larger (i.e., phenological development is more advanced) than the mean stage in Willamette Valley populations at a given physiological time (accumulated degree days). A bootstrap variance was obtained for testing to account plant-to-plant variation in insect phenology within populations, because the individual plant represents the unit of random sampling. The 2011 model required an arcsine square-root transformation after translation of the outcome to a [0,1] interval to more aptly meet the linear regression assumptions of normally-distributed errors with equal variance.

Common garden experiment

We collected eggs in the summer of 2004 from the following hosts and locations: J. vulgaris at Baskett Slough (44.964°N, −123.258°W) in the Willamette Valley, J. vulgaris at Neskowin (45.131°N, −123.964°W) on the Oregon Coast, and S. triangularis at Santiam Pass (44.426°N, −121.850°W) in the Cascade Mountains of Oregon.

We then reared the insects fed ad libitum on J. vulgaris foliage and floral parts in summer of 2004 in a common (greenhouse) environment under conditions of fluctuating temperature and photoperiod and a constant range of humidity (light, 16 h at ∼22°C; dark, 8 h at 12°C; humidity, ∼70–100%). By fixing the environment, we hoped to isolate the component of phenotypic variation owing to differences in genotype. However, by collecting insects in the egg stage from their naturalized-home environment, we did not strictly control (at this phase of our investigation) for possible effects of the parental environment that might be transmitted to offspring (see next section). The sample unit was a cluster of 10 full-sib larvae in a 0.5 × 0.5 × 0.5 m cage covered with nylon ‘Leno weave’ netting covers (open spaces in netting were 0.6 × 1.0 mm) (BioQuip, Rancho Dominguez, CA, USA). We randomized the location of each sample unit on each greenhouse bench. We replicated clusters of insects from each location unevenly because we were constrained by availability of disease-free insects from each environment. Replication was as follows: Coast = 5, Willamette Valley = 2, and Cascades = 4 clusters of 10 larvae each. We started the experiment by introducing 10 first-instar larvae that had just hatched from the egg stage into each cage. We made daily observations on survival (number of individuals surviving per cage) and developmental time (mean number of days to complete development from egg to pupa per cage).

We carried out statistical analysis of variation in development time and survival as follows. We tested whether survival varied among geographic locations (Cascades, Willamette Valley, Coast) (Table 1), using a one-way anova on the number of surviving insects in each experimental unit (cluster in a cage). We regressed mean development times for clusters reared in the laboratory on the accumulated annual degree days available in the naturalized field environment of each population, assuming a linear model, and tested for significant slope to assess (i) whether variation in mean development times among three populations was greater than expected from variation within populations and (ii) whether there was a trend in mean development times among populations related to available degree days in the naturalized field environment.

Crosses within and between populations

We analyzed phenological development in full-sib families of offspring with unique parents sampled in 2010 from populations at two extremes of the environmental gradient (Fig. 1, Table 1), Willamette Valley and Cascades. These were the same populations sampled earlier in 2004 for the common garden experiment.

Figure 1

Environmental temperature curves for locations along the environmental gradient from Coast to Cascades.

We reared the insects individually for a generation in a common environment to control for possible effects of the parental environment and to eliminate a host-specific pathogen Nosema tyriae from laboratory populations, using the methods of Karacetin (2007). We reared T. jacobaeae individually from egg to pupa on foliage and floral parts of J. vulgaris and then overwintered pupae in mulch in a cold frame outdoors at OSU. On 19 January 2011, we transferred pupae to the laboratory for transition to the adult stage.

We randomly crossed emerging adults from within populations to create pure-bred family lines (Cascade F × Cascade M and Valley F × Valley M) and between populations to create hybrid family lines (Cascade F × Valley M and Valley F × Cascade M). Random mating avoided assortative mating (e.g., early with early, late with late), which can bias heritability estimates, although there is no necessary correlation between timing of emergence in adults and the phenotypic trait measured, juvenile development time. We reared offspring from these crosses individually in cups (to minimize sibling resemblance because of a common environment) from egg to pupa on J. vulgaris foliage and flowers in an environmental growth chamber (Hoffman Manufacturing) at a constant, optimal temperature (20°C), long days 18-h L/6-h D and high relative humidity (87% in chamber, nearly 100% in cup with foliage), changing food daily, recording survival, and development time.

We analyzed the results of crossing experiments as follows. First, we screened for evidence of the patterns of inheritance. If offspring from hybrid (between population) crosses were intermediate to offspring from pure-bred (within population) crosses in juvenile development time, then we would conclude development time is a quantitative trait (controlled by multiple genes, each of small effect). If hybrids bore a stronger resemblance to either parent, then we would additionally conclude development time is a sex-linked trait. Second, we screened for evidence of a genetic differentiation between populations in juvenile development time by comparing cumulative distributions of development times from egg to pupal stage, given that individuals pupated. Third, broad-sense (full-sib) heritability for development time was calculated by dividing two times the phenotypic variance component of family by total phenotypic variance (Roff 1997), assuming no correlation between phenotypic values of mates, dominance effects, or effects of a common environment. Standard error of heritability for families with unequal family sizes was calculated as recommended by Roff (1997). We excluded families with very small sample sizes (ni < 3), leaving N = 21 families distributed across four cross types (Cascade F × Cascade M, Valley F × Valley M, Cascade F × Valley M, Valley F × Cascade M) with ni = 3–15 individuals per family and a total of 139 individuals. We used restricted maximum likelihood (REML) to estimate variance components among families (AF) and among progeny within families (AP) when data are unbalanced (Hartley and Rao 1967) for a full-sib design assuming no effects of a common environment (Roff 1997). Note that the variance component AP is just the residual error variance as we do not have enough replication at the progeny level to estimate this variance component separately.

We performed statistical tests of random null hypotheses of founder effects and genetic drift against the alternative hypothesis of natural selection as explanations for observed differences in phenology. To discriminate between founder effect and directional selection, we compared variances of the mean phenology times by plant for mountain populations (Cascades and Coast Range) against the variance from the source population from the Willamette Valley, and we compared the mean phenology times as described earlier. We next applied the rate test of Lande (1977) to distinguish between the null hypothesis of evolution by random genetic drift from the alternative hypothesis of natural selection. The rate test requires information on observed phenotypic variance between populations after t generations and the variance within populations σ2,heritability h2, time since divergence t measured in generations, and effective population size Ne. The ratio of variances to be tested is

The degrees of freedom associated with the F statistic are as calculated by Lande (1977).

Results

The environmental gradient is complex, involving variation in topography, temperature, moisture, and host plant species (Table 1). Temperatures are lower, and growing seasons are shorter at high compared to low elevations; T. jacobaeae's growing season for the Willamette Valley (∼1400 degree days) is >3 × the length of the growing season in the Cascades (∼400 degree days) (Fig. 1).

We found variation in phenology among T. jacobaeae populations along an environmental gradient. The 2010 data reveal differences in phenology between Willamette Valley and Cascades populations. Plotting phenology (as number of individuals in the 5th instar) on a physiological timescale (degree days) illustrates that the differences in juvenile development times between Willamette Valley and Cascades populations are conserved (Fig. 2A), and not simply related to differences in the temperature regime in relation to calendar time that a population experiences. We confirmed significant differences by comparing the mean times in stage for Willamette Valley and Cascades populations using methods described by Murtaugh et al. (1987). We obtained bootstrap estimates of the mean difference in mean times in stage (degree days) between regions of 52.23 (95% CI: 18.35, 86.10) for the egg stage, 27.50 (95% CI: 18.57, 36.42) for stage L1, 55.91 (95% CI: 33.63, 78.18) for stage L2, 66.18 (95% CI: 43.02, 89.33) for stage L3, 55.80 (95% CI: 33.03, 78.57) for stage L4, and 65.24 (95% CI: 48.29, 82.18) for stage L5. The mean differences are substantially larger than zero for all stages, indicating a systematic difference in phenology between a Cascades and a Willamette Valley population. The results of the one-sided t-test indicate significant differences between degree days for all stages, a significant difference across all stages between the two populations (chi-square test, χ2 = 186.36, df = 12, P < 0.0001), and mean times in stage that are less for a Cascades population compared to a Willamette Valley population (Fig. 2B).

Figure 2

(A) Number of individuals in the fifth instar by population. The curves illustrate that 5th instar larvae of Tyria jacobaeae (cinnabar moth) from the Cascades complete their development earlier than those from Willamette Valley. (B) This difference in phenology is confirmed statistically (see P -values for the one-sided test that mean time in stage is less in the Cascades than in the Willamette Valley) and visualized by boxplots of times in stage (degree days) for each instar.

Mean stage at a given physiological time was modeled as a function of degree days, regions, sites within regions (2011 only), and interactions between degree days and regions. The multiple regression analysis resulted in a 2010 model of mean stage of −0.57 + 0.02*Degree Days + 0.33*I (Valley) − 0.01*Degree Days*I(Valley), where I(Valley) is 1 when the outcome comes from the Willamette Valley and 0 otherwise. The 2011 multiple regression model for mean stage is 0.14 + 0.0033*Degree Days + 0.15*I(Valley) + 0.12*I(BS) − 0.04*I(PCT) − 0.0010*Degree Days*I(Valley), where I(BS) and I(PCT) are indicators that the outcome was obtained in the Baskett Slough or Pacific Coast Trail PCT sites, respectively. Mean stage at a given physiological time was higher (i.e., more advanced) in Cascades compared to Willamette Valley regions in both 2010 with one population per region (bootstrap z-test, z-statistic = 8.63, P < 0.0001) and in 2011 with two replicate populations per region (bootstrap z-test, z-statistic = −2.45, P = 0.0071), indicating that phenological development was faster in Cascades populations compared to the Willamette Valley populations. Plots of 2010 mean stage against degree days indicate mean stage near 1 for degree days <200 and mean stage near 6 for degree days over 500, and these extremes serve to anchor the slopes. Scatter plots of mean stage versus physiological time reveal high levels of plant-to-plant variation in the phenologies of insects for Valley and Cascades (Fig. 3).

Figure 3

Scatter plots of mean stage by time in stage (degree days) for two Cascades (Pacific Crest Trail PCT and Wasco Lake) and two Willamette Valley (Alsea Highway and Baskett Slough) populations in 2010 and 2011 with multiple regression lines.

The shorter juvenile development times of the Cascades population match the constraint of a shorter growing season (annual accumulation of ∼500 degree days in 2010) in the Cascades. The longer juvenile development times of the Willamette Valley population do not match the longer growing season there (annual accumulation of ∼1400 degree days in 2010), possibly because phenology is constrained by a lack of genetic variation for longer development times, subject to fluctuating selective pressures owing to a variable abiotic environment, or subject to selection by environmental factors other than temperature.

We found evidence of a strong, positive relationship between mean development time for individuals within full-sib families reared in the greenhouse environment and the available degree days (the amount of heat available for development based on 30-year normal temperatures) in their naturalized-home environment (Fig. 4). The slope of the regression (Y = 21.2886 + 0.0048 X, P = 0.0005) is highly significant and the R2 value indicates that 73% of the variation in development time is explained by its regression on available degree days. We conclude that there is genetic differentiation among three populations (Coast, Willamette Valley, and Cascades) in juvenile development time and that genetically based variation among populations in their development-time phenotype is strongly, positively related to the available degree days in the population's naturalized-home environment.

Figure 4

The development time from egg to pupa measured for populations in the laboratory environment increased linearly with the physiological time available for development in the field environment of the population's naturalized-home environment. Points represent means for survivors from sample units of 10 larvae.

We found no significant variation in the number of surviving larvae among the three geographic populations when reared in a common environment (one-way anova, F = 3.23195, P = 0.10993). Mean survival from egg to the pupal stage was 80%.

Crosses within and between extreme populations

The development times of offspring from hybrid (between population) crosses were intermediate to offspring from pure-bred (within population) crosses in juvenile development times, indicating that development time is a quantitative trait (controlled by multiple loci, each of small effect) (Fig. 5, Table 3). The juvenile development times of hybrids did not bear a stronger resemblance to either parent, indicating sex linkage in this trait is unlikely (Fig. 5, Table 3). We excluded the possibility of maternal effects, a possibility not entirely ruled out by our common garden experiment, by rearing insects in a common environment for a generation prior to conducting the crosses. The contrast of juvenile development times for offspring from two pure-bred (within population) cross types confirms a genetic differentiation of Willamette Valley and Cascades populations; offspring from pure-bred Cascades parental crosses reached the overwintering pupal stage 5 days earlier on average than offspring from pure-bred Willamette Valley parental crosses when reared at a constant temperature of 20°C (Fig. 5, Table 3). That represents a difference of 57 degree days, a number indistinguishable from the difference of 65.24 (95% CI: 48.29, 82.18) degree days accumulated by the 5th instar estimated for Willamette Valley and Cascade populations in the field for 2010. Estimates of heritability (with standard errors) for developmental time were h2 = 1.69(0.10) for families in the Willamette Valley and h2 = 1.84(0.04) for families in the Cascades (in both cases for families with at least three progeny). These estimates are likely inflated by sampling error because of the small number of families or by some combination of unmeasured effects because of correlation between phenotypic values of mates, dominance, and common environment (see Discussion).

Figure 5

Cumulative distribution function for juvenile development times from egg to pupa for offspring from pure-bred crosses (Cascades, Willamette Valley) and hybrid crosses (Cascades F × Valley M, Valley F × Cascades M) reared at a constant temperature (20°C). F = female and M = male. P(x) is the probability of a development time ≤x. Cumulative curves represent pooled data from independent families derived from crosses reported in Table 3.

Table 3

Number of families (n), weighted family means () in days, and standard errors (SE) for juvenile development times of offspring from pure-bred and hybrid lines

Cross	n		SE
Pure-bred lines
Cascades F × Cascades M	8	33.71	0.36
Willamette Valley F × Willamette Valley M	10	38.11	0.59
Hybrid lines
Cascades F × Willamette Valley M	7	35.54	0.30
Willamette Valley F × Cascades M	6	35.86	0.51

Statistical tests on the pattern and rate of evolution of phenology helped distinguish natural selection from founder effects and random genetic drift as possible explanations of the evolution of development time. Both founder effects and directional selection predict a reduction in phenotypic variance. As predicted, a one-sided test of variances of mean degree days by plant (insects on a plant approximate full-sibs because of oviposition of egg masses) within Cascades and Willamette Valley regions showed that both mountain populations are less variable in phenology than the original Willamette Valley population (F test, F47,48 = 0.5727, P = 0.0290 for Willamette Valley versus Cascades; F48,48 = 0.5627, P = 0.0245 for Willamette Valley versus Coast Range). The directional selection hypothesis predicts consistently shorter phenologies with independent translocations of T. jacobaeae from the valley source population to the mountains (Coast Range, Cascades), while founder effects predict no consistent direction to changes in the mean. Contrary to the founder effect hypothesis, we found consistently shorter phenology in mountain populations compared to the valley. We conclude it is unlikely that founder effects explain the observed evolution of phenology, leaving directional selection as the most plausible alternative.

To test the null hypothesis of random genetic drift against the alternative hypothesis of natural selection using the rate test of Lande (1977), we require estimates of phenotypic variance in development time within populations σ2 and between populations for the populations compared, heritability h2, time since divergence t, and effective population size Ne. The observed phenotypic variance between populations after t generations is = 4.28, and the variance within populations is σ2 = 10.47. The effective population size Ne, generally much less than the actual population size Na, was estimated as Ne = 1000. This is an approximation from records kept by the Oregon Department of Agriculture, which indicate the actual release population size averaged 1000 individuals, and 1–5 releases of this size were made annually between 1979 and 1984 in or near the study area. The time since divergence t was estimated as t = 23, the number of generations between 1987 when outbreaks of T. jacobaeae occurred in the Cascades after translocation from the Willamette Valley, and the first year of our observations on phenology in 2010. The estimate of heritability from the Willamette Valley was used, although estimates were similar for the Cascades region. The numerator degrees of freedom are n − 1, where n = 3 is the number of populations. Because the denominator of the F-statistic is not known exactly and the distributions of each of the involved factors are also unknown, the denominator degrees of freedom are conservatively calculated at 10 000 to reduce inflation of the type I error rate (Lande 1977). The F-statistic is calculated as F3,10 000 = 7.24, and the one-sided test yields a P-value <0.0001. Therefore, we conclude that it is extremely unlikely that the observed population differences evolved by genetic drift, leaving directional selection still standing as the most plausible alternative.

Discussion

The results of our study can be summarized as follows: (i) field observations comparing phenotypic variation in the phenologies at each stage for replicate populations in the Willamette Valley and mountains (Coast Range, Cascades) confirmed mountain populations have shorter phenologies corresponding to shorter growing seasons, (ii) a common garden experiment confirmed genetic determination of the phenotypic differences in juvenile development times between Coast, Willamette Valley, and Cascade Mountain populations and further showed that development times decrease linearly with decreasing length in the growing season along the environmental gradient, (iii) a laboratory experiment rearing offspring from parental crosses within and between Valley and Cascades populations further confirmed genetic determination of the phenotype by showing development time is a quantitative trait with high heritability and no evidence of sex linkage, and (iv) a rate test rejected the null hypothesis of evolution by random genetic drift in favor of the alternative hypothesis of natural selection as an explanation for the evolution of development time. In the following paragraphs, we discuss the evidence bearing on the hypothesis of natural selection (phenotypic variation, heritability, and the relative roles of genetic drift, founder effect, and natural selection in the evolution of phenology) and the implications of these findings.

Our analysis revealed the magnitude and direction of phenotypic differences in phenology among populations and confirmed genetic determination of the phenotype. Mountain populations developed faster than Valley populations at each stage from egg to pupa. Magnitudes of differences in phenology between populations varied with the differences in elevation – differences were greater for the comparison of Willamette Valley (elevation 87 m) and Cascades (1572 m), less for the comparison of Willamette Valley and Coast Range (877 m). The magnitudes of phenotypic differences in phenology in field and laboratory, 65.24 (95% CI: 48.29, 82.18) and 57 degree days, respectively, are strikingly similar under natural conditions in the field (where contributions to phenotypic variance from both genotype and environment are expected) and in a common laboratory environment (where the relative contribution owing to genotype is magnified over that owing to environment). Our analysis also confirmed a large standing stock of phenotypic variation in development time on which selection might act (Barrett and Schluter 2007) in all populations that we examined (Coast, Coast Range, Willamette Valley, and Cascades). This is noteworthy because genetic bottlenecks during screening of candidate agents in quarantine and in their initial field release may reduce genetic variation unless the time spent in bottlenecks is short and the rate of increase in release populations is high (Roderick and Navajas 2003).

Our analysis confirmed high and similar levels of broad-sense heritability for development time (egg to pupa) within two populations (Willamette Valley h2 = 1.69 and Cascades h2 = 1.84). Our estimates of broad-sense heritability (which are expected to range between 0 and 1) are clearly inflated, possibly due to sampling error with the small number of families across regions or possibly due to unmeasured effects of assortative mating, dominance, or a common environment. To take the example of bill depth in Darwin's Finches in the Galapagos (Boag and Grant 1978), the value of heritability estimated by intraclass correlation in full-sibs (h2 = 1.42) is 1.73 times the value from offspring–midparent regression (h2 = 0. 82, assumed to be the true value), but a full accounting of the effects of a correlations between phenotypic values of mates, dominance, and a common environment can explain this discrepancy (Falconer and Mackay 1996). Even at half their value, our estimates of heritability are high compared to those reported for life-history traits of invertebrates in general (Mousseau and Roff 1997) but closer in magnitude to heritability (h2 = 0.60) estimated by parent-offspring regression for another phenological trait (heat requirement for diapause) in a species related to the cinnabar moth, that is, in the same family (Hyphantria cunea, Arctiidae) (Morris and Fulton 1970).

The results led us to reject the random null hypotheses of evolution by founder effects or by random genetic drift in favor of the alternative hypothesis of natural selection. The founder effect hypothesis was rejected: While it predicts a reduction in phenotypic variance in phenology consistent with that observed, its does not account for the consistent direction of changes in mean phenology (shorter phenology in mountain populations compared to the valley source population) that we observed. The relatively large release population sizes (1000–5000 individuals) and the relatively low frequency of genotypes controlling short phenology in the parent valley population (Fig. 2) mean such genotypes have a small chance of being fixed by a founder event in the mountains. The hypothesis of genetic drift can be rejected based on the results of the rate test (Lande 1977). We believe that the chief weakness in our analysis under the drift hypothesis, apart from our estimates of heritability as discussed above, is using Ne = 1000 to estimate effective population size Ne; the actual value is unknown and expected to be less than the actual population size (Lande 1977). Any Ne < 320 individuals would yield an F-statistic that is not significant in this case. Were we to decrease our estimates of heritability (by sampling larger numbers of families encompassing a broader spectrum of variation in the regions or by accounting for effects of correlation between phenotypic values of mates, dominance, or a common environment) that would only strengthen the inference that populations did not evolve by genetic drift, leaving directional selection as the most plausible hypothesis? To add further evidence of adaptive evolution to climate, we recommend (i) more replication of populations within and between climatic regions to establish the generality of the pattern of phenotypic variation that we observed, namely that development time decreases with decreases in the growing season, (ii) a reciprocal transplant experiment to decisively distinguish local adaptation from alternative hypotheses (Kawecki and Ebert 2004), and (iii) thermal reaction norms to estimate environmental effects on heritability and describe the pattern of phenotypic expression of a each genotype across a range of temperatures (Angilletta 2009).

What are the implications for management?

The key points to appreciate from our study are that there is genetic variation for ecologically important traits within populations of biological control organisms, and adaptive evolution can and does occur rapidly within populations on ecological timescales. A control organism may acquire new hosts that were previously outside its climatic range through adaptive evolution in local population phenology.

World-wide reviews of risk-benefit-cost analysis for the introduction of classical biological control agents against weeds (Sheppard et al. 2003) and insects (van Lenteren et al. 2006) suggest that there is currently no regulatory framework for evaluating evolutionary potential of biological control organisms prior to their release in new environments. Currently populations, not species, are the level of biological organization permitted for introduction (Hufbauer and Roderick 2005). It is at the population level that evolution by natural selection occurs, and we strongly believe that the population is the place to start when designing appropriate regulations to assess evolutionary potential. Adaptive evolution in biological control organisms is most likely in cases like T. jacobaeae with human-assisted spread, genetic variation, strong suppression of the target host, a broad host range, and maladaptation to the new environment (Thrall et al. 2011). Evolution has a large stochastic component as a result of forces such as random mutation and genetic drift, and no predictive schemes will ever reach 100% accuracy.

Approaches for predicting and describing evolutionary responses of biological control organisms to new environments can be gleaned from reviews of the adaptive potential of organisms in relation to global climate change (Hoffmann and Sgrò 2011). Approaches include (with illustrative examples added from the biological control literature) (i) longitudinal studies testing for genetic changes in populations, (ii) spatial studies across climatic gradients involving transplants or common garden experiments (the present study), (iii) standing quantitative genetic variation estimates within populations (Wajnberg 2004), (iv) quantitative genetic variation estimated through selection experiments (Hopper et al. 1993), (v) experimental evolution in simulated environments (Kraaijeveld 2004), (vi) evidence of loss of function of candidate gene/protein, and (vii) genetic variation in candidate genes for traits pointing to potential for evolution. Foundational research on model biological control systems should combine realistic models of evolution (Perkins 2012) and empirical approaches such as those above to develop a more reliable basis for understanding, predicting, and managing biological control systems. Certainly, evolutionary biology should be part of the training of all biological control scientists (Futuyma 2000).