Rapid evolution of phenology during range expansion with recent climate change

Although climate warming is expected to make habitat beyond species’ current cold range edge suitable for future colonization, this new habitat may present an array of biotic or abiotic conditions not experienced within the current range. Species’ ability to shift their range with climate change may therefore depend on how populations evolve in response to such novel environmental conditions. However, due to the recent nature of thus far observed range expansions, the role of rapid adaptation during climate change migration is only beginning to be understood. Here, we evaluated evolution during the recent native range expansion of the annual plant Dittrichia graveolens, which is spreading northward in Europe from the Mediterranean region. We examined genetically based differentiation between core and edge populations in their phenology, a trait that is likely under selection with shorter growing seasons and greater seasonality at northern latitudes. In parallel common garden experiments at range edges in Switzerland and the Netherlands, we grew plants from Dutch, Swiss, and central and southern French populations. Population genetic analysis following RAD‐sequencing of these populations supported the hypothesized central France origins of the Swiss and Dutch range edge populations. We found that in both common gardens, northern plants flowered up to 4 weeks earlier than southern plants. This differentiation in phenology extended from the core of the range to the Netherlands, a region only reached from central France over approximately the last 50 years. Fitness decreased as plants flowered later, supporting the hypothesized benefits of earlier flowering at the range edge. Our results suggest that native range expanding populations can rapidly adapt to novel environmental conditions in the expanded range, potentially promoting their ability to spread.


Introduction
The redistribution of species in response to climate change (Chen, Hill, Ohlemüller, Roy, & Thomas, 2011;Parmesan, 2006) has focussed ecologists and evolutionary biologists on the processes determining population spread. The persistence of many species will depend on their ability to migrate (Thuiller et al., 2008), and the resulting range shifts will have broad implications for both natural systems and human welfare (Pecl et al., 2017). A central unresolved question in the study of range expansion is whether evolutionary changes that occur in contemporary time facilitate the spread of range shifting populations and thereby contribute to persistence (Hoffmann & Sgrò, 2011;Urban et al., 2016). In particular, as populations spread, they face novel environments, and thus adaptation to local conditionsor the lack thereof -could profoundly affect species' capacity to establish and therefore migrate (Garcia-Ramos & Rodriguez, 2002;Gilbert et al., 2017). Forecasting species' responses to climate change could therefore benefit from a better understanding of the role of adaptive evolution during native range expansions (Urban et al., 2016).
Even if the northward range expansion of climate change migrants is ultimately set in motion by warming conditions (Parmesan & Yohe, 2003), continued expansion may require evolution in response to other environmental variables experienced in the expanded range.
For example, day length, light quality, and seasonal variation in climate all change markedly with latitude (Saikkonen et al., 2012;Taulavuori, Sarala, & Taulavuori, 2010), and may present range expanding populations with conditions not experienced within their historic range. Consequently, the cues used by many plant and animal species to time key events in their life cycle may no longer match the temperature conditions optimal for life cycle transitions (Visser, 2008), and populations may be maladapted to the environment of the expanded range, limiting further range expansion. Phenology in particular, is expected to be under strong selection as populations spread northward because timing reproductive events

Accepted Article
This article is protected by copyright. All rights reserved. combination of genetic drift, gene flow, and natural selection. For example, range edge populations may sometimes be pre-adapted to local conditions beyond the historical range limit by originating from similar environments within the historical native range (Colautti & Lau, 2015;Oduor, Leimu, & van Kleunen, 2016).
Here, we evaluate rapid evolution during the recent native range expansion of Dittrichia graveolens (L.) Greuter (hereafter Dittrichia), an annual plant species in the Asteraceae with a native distribution in the Mediterranean (Brullo & De Marco, 2000). Consistent with climate warming over the last 50 years, Dittrichia has rapidly expanded its range northward over this period (Rameau, 2008), spreading from France into Germany (Brandes, 2009;Jäger, 2017) and then into the Netherlands (Stouthamer, 2007). This northward expansion subjects the plants to a reduction in light availability, especially towards the end of the year (Appendix S1), which may be particularly important because the species flowers in late summer and fall in the core of its range in France (Rouy, 1903). For many plant species, growth is limited by low temperatures, and particularly frost, at the end of the growing season (Larcher & Bauer, 1981). In addition to its northward range expansion, Dittrichia has also expanded its range eastward into Switzerland (Ciardo & Delarze, 2005;Lauber, Wagner, & Gygax, 2012), which is located at a similar latitude and altitude as the historical northern range edge (Fig. 1). The recent nature, large spatial extent, and detailed historical record of the spread of Dittrichia provide a unique opportunity to study the role of rapid adaptation in native range expansions. Moreover, Dittrichia's annual life cycle allows for meaningful evolutionary change to occur within a few decades.
We tested the hypothesis that rapid evolution in response to novel environmental conditions promotes native range expansion with climate change. For Dittrichia, we specifically expected the evolution of earlier flowering time with its northward range expansion. To test the timing and northward shift coincide with climate warming in Europe in the late 20 th century (Lenoir, Gégout, Marquet, Ruffray, & Brisse, 2008). In the 1980s, Dittrichia was reported on roadsides in south-west Germany (Garve & Garve, 2000;Nowack, 1993) and the northwestern Ruhr district (Dettmar & Sukopp, 1991), and began spreading along highways from there around 10 years later (Garve & Garve, 2000;Nowack, 1993;Radkowitsch, 1996). In the Netherlands, it has been expanding its range along highways since approximately 2005 (NDFF, 2015;Sparrius & Van Strien, 2014;Stouthamer, 2007). In addition to its northward range expansion, the species also recently extended its distribution eastward in France (Antonetti, Kessler, Brugel, Barbe, & Tort, 2006), and spread into Switzerland along highways (Lauber et al., 2012). By 2003, it was widely distributed around Lake Geneva (Ciardo & Delarze, 2005). In this study, we examine plant performance at the Dutch and Swiss range edges (Fig. 1).

Seed collection
To establish the common gardens and provide plant material for molecular analysis of population structure, we collected seeds from the core of the range in France, and from the expanding range edges in Switzerland and the Netherlands. We collected seeds from four populations in southern France, representing the Mediterranean conditions in much of the native range, and from nine populations in central France, along a latitudinal gradient towards the historic northern range edge. From the range edge in Switzerland, which is located at a similar latitude and altitude as the populations in central France, we collected seeds from three populations. Finally, we collected seeds from three populations in the Netherlands, representing the most novel northern latitude populations (Fig. 1, Table 1). The Netherlands receives markedly less photosynthetically active radiation as compared to Switzerland, central or southern France (Appendix S1, Fig. S1), meaning that if daylight is limiting, Dutch plants must complete their life cycle earlier in the year to receive comparable

Accepted Article
This article is protected by copyright. All rights reserved. collaboration with the Genetic Diversity Center (GDC), ETH Zurich. Keygene N.V. owns patents and patent applications protecting its Sequence Based Genotyping technologies.

Common garden experiments
To examine the reproductive timing and fitness of plants from all regions of origin under range-edge conditions, we set up parallel common garden experiments in Wageningen, the Netherlands, and Zurich, Switzerland. The gardens were established in 2016 and located at the range edges because these are the locations at which understanding how phenology has evolved and influences fitness is most relevant to future range expansion. To reduce the impact of maternal effects across all regions of origin and produce seeds for use in our main common garden experiments, we first grew plants from seed in a common environment for one (Dutch seeds) or two (Swiss and French seeds) generations of self-fertilization.
Maternal families in our study descend from a single field plant (Appendix S2). Ultimately, we used seeds from three Dutch, three Swiss, six central French and three southern French populations to establish the common gardens (Table 1). Three populations from central France were excluded because they were located very close to another population within the same city, and one population from southern France was excluded due to limited seed availability. For each of the 15 populations, we grew three individuals from four maternal families, comprising 12 plants per population. This design was replicated in the two gardens, using the same seed sources (Appendix S2).
To produce the seedlings that would eventually be planted into these gardens, on June 13 th (2016), seeds were placed on moist germination paper (Zurich) or gamma-sterilized soil (Wageningen) in transparent boxes in growth cabinets set to 28°C with a 16h day / 8h night cycle. After approximately four days, seeds began germinating and were transferred onto soil in seedling trays in a greenhouse. On July 11 th , when the seedlings were strong enough,

Accepted Article
This article is protected by copyright. All rights reserved. they were transplanted into 6L plastic pots (21.2 cm in diameter and 16.8 cm high) filled with a common sterilized sandy loam soil collected from a former agricultural field in the Netherlands (Beneden-Leeuwen; 51.89°N, 5.56°E). In Wageningen, the pots were placed on a tarp in an open field (51.99°N, 5.67°E, 13m); in Zurich, in wooden beds on a roof terrace (47. 38°N, 8.55°E, 460m). We used a randomized block design, with each of three blocks containing one plant from each maternal family. For the duration of the experiment, we provided extra water when plants were at risk of drying out due to hot or dry weather conditions.
To quantify phenology, we conducted a census twice a week on fixed days, and recorded the date each plant first produced flowers (at least one yellow floret visible). The day of first flowering (which is more accurately measured than budding or fruiting date) was our measure of phenology, quantified as days since plants were placed in the common garden.
Plants were harvested on December 1 st and 2 nd , when development had slowed and plants started to senesce. For each plant, we counted the number of fruiting heads (those with pappus-bearing achenes) on the plant. We defined plant fitness as the number of fruiting heads multiplied by seed viability, the methods for which are explained next. To evaluate the possible fate of plants after our harvesting date, we left 11 of the latest-flowering plants in the common garden until early February, but only three of these 11 produced any viable seeds.
After harvesting the experiment, we evaluated seed viability. For each plant, 30 seeds were selected haphazardly from fruiting heads that were collected at the time of harvest. These fruiting heads were open-pollinated, so influence from paternal individuals from different populations is possible. Nonetheless, seed development is strongly linked to the flowering phenology and maternal investment of the mother plant. We first tested the germination

Accepted Article
This article is protected by copyright. All rights reserved.
fraction of the seeds on moist germination paper in transparent boxes in growth cabinets set to 28°C with a 16h day / 8h night cycle. Germinated seeds were counted and removed after 15-16 days. The remaining seeds were covered in a 250mg/L gibberellic acid (GA 3 ) solution and further germinants were counted after 24 hours. We lastly estimated the viability of the still non-germinated seeds at the region of origin level (Netherlands, Switzerland, central France and southern France) for each garden. To do so, we selected one seed at random from 30 individuals from each region, cut the seeds in half and stained them with a 0.25% tetrazolium solution. Viable embryos were counted after 2.5 hours. The final viability rate per individual plant was thus the number of germinated seeds at the plant level plus the projected number of viable, non-germinated seeds (based on 30 seeds from the plant's region of origin and garden), divided by the total number of seeds tested.

Data analysis
(i) Population structure and genetic diversity Bioinformatics analysis of the RAD sequencing data (details provided in Appendix S3) yielded 13,582 single nucleotide polymorphisms (SNPs) across 190 individuals. We evaluated population structure using STRUCTURE v2.3.4 (Pritchard, Stephens, & Donnelly, 2000). To prevent bias in the clustering algorithm, we used the haplotype data and a single individual per maternal family (153 in total) for our analysis, selecting the half-sib with the least amount of missing data. The analysis was conducted with two to eight genetic clusters (K), and replicated 10 times using the admixture model with correlated allele frequencies.
We visualized the output using Structure Harvester (Earl & vonHoldt, 2012) and determined the optimal K using the Evanno et al. (2005) method. In range expanding populations, we may expect many loci not to be in Hardy-Weinberg equilibrium (HWE). We therefore tested the robustness of our results to this assumption by rerunning the STRUCTURE analysis

Accepted Article
This article is protected by copyright. All rights reserved.
using only the SNPs in HWE, and by computing a neighbor-joining tree, which does not depend on a population model (details in Appendix S3).
To infer changes in genetic diversity during range expansion, we examined the genetic diversity within three genetically distinct groups in the core of the range (identified in the population structure analysis) and at the two expanding range edges. To ensure equal numbers of individuals within these five groups, we selected 18 individuals at random per group, distributed uniformly across the available populations. Nucleotide diversity π was calculated for each SNP using VCFtools (Danecek et al., 2011), averaged across the loci and divided by the mean fragment length of 121 bp.
(ii) Differentiation in phenology and adaptation to range edge conditions All statistical analyses were conducted in R v3.4.0 (R Core Team, 2017). We analyzed phenological differences between geographic origins using a mixed-effects Cox proportional hazards model with a Gaussian distribution of random effects (coxme package, Therneau, 2015). In our analyses, we compared four distinct regions of origin: central and southern France in the core of the range, and the two range edges in Switzerland and the Netherlands. Population and maternal family (nested in population) were treated as random effects. We first fit a global model with the day of first flowering as the response variable and region of origin, common garden location, and their interaction as fixed effects. Main and interaction effects were evaluated using Type II partial-likelihood-ratio tests (car package, Fox & Weisberg, 2011). Because there was a significant interaction between origin and garden, we then fit separate models per garden and analyzed pairwise differences between regions of origin using Tukey contrasts (multcomp package, Hothorn, Bretz, & Westfall, 2008).

Accepted Article
This article is protected by copyright. All rights reserved.
To estimate selection on phenology in each garden, we assumed flowering day to be a quantitative trait and standardized it to zero mean and unit variance, separately for each garden. Relative fitness was calculated by dividing each plant's fitness by the average fitness of each garden. We then estimated natural selection on phenology at the garden level using fitness splines (Schluter, 1988), fitting a nonlinear fitness function to all individuals using generalized additive models in R (R code based on Colautti & Barrett, 2013a;package mgcv, Wood, 2011). Additionally, we estimated selection on flowering day separately for individuals from each region of origin, because these regions differed in their range of phenology. Here, we used linear regression as in a classic selection analysis (Lande & Arnold, 1983).
To analyze the effects of origin, garden, and their interaction on fitness, we fit a mixedeffects negative binomial model using maximum likelihood (function glmer.nb in package lme4, Bates, Maechler, Bolker, & Walker, 2015), which is appropriate for overdispersed count data. We tested the main and interaction effects of origin and garden on fitness using likelihood ratio tests, and removed the interaction since it was not statistically significant. We then tested pairwise differences between regions using Tukey contrasts in an additive model with garden and origin as fixed effects.

Spread history
Population genetic analysis revealed three distinct genetic clusters of individuals (STRUCTURE analysis based on all SNPs; Fig. 1

Accepted Article
This article is protected by copyright. All rights reserved. To further explore the importance of early phenology for plant fitness, we estimated phenotypic selection on flowering day. In both common gardens, fitness declined as plants flowered later, approaching zero for plants flowering from early October onwards (Fig. 3).
Although plants from southern France, which flowered last, had the lowest fitness, we also found selection for earlier phenology at the range edges among plants from central France, the putative region of origin of the populations spreading at both range edges; Fig. 3).
Statistically significant selection on earlier flowering time was found in nearly all populations placed at the range edge (Fig. 3).
In both common gardens, plants of Dutch, Swiss, and central French origin had comparable fitness, whereas the fitness of plants from southern France was much lower and often zero (pairwise difference with all other regions: p < 0.001; Fig. 4). Consistent with expectations based on their phenology, the fitness of Dutch plants in the common garden in Wageningen (the Netherlands) was slightly higher than the fitness of central French plants (Fig. 4, left), although there was no statistical support for this result due to the limited number of sampled populations. Overall, plants in the Wageningen garden had almost three times higher fitness than plants in the Zurich garden ( 2 (1) = 169, p < 0.001; Fig. 4).

Accepted Article
This article is protected by copyright. All rights reserved.
populations, which more likely originated from southern France (genetic cluster K1), may have recently evolved earlier phenology as they spread north.
Several lines of evidence support our hypothesis that the differentiation in phenology among central France, Swiss, and Dutch populations is the result of selection for earlier flowering time over the fifty years of range expansion. Since de novo mutations are unlikely to play a large role on this timescale, sufficient genetic variation in phenology would be required within the central France populations, providing the material on which natural selection could act (Barrett & Schluter, 2008). Indeed, the phenology variation among the central France plants was such that the earliest individuals were comparable in timing to those of the Dutch and Swiss populations (Fig. 3). Moreover, flowering time is highly heritable in this species (h 2 = 0.86; Appendix S4). Finally, we found selection for early phenology in both range edge common gardens (Fig. 3). The decline in fitness as plants flowered later matches our expectation that early flowering is beneficial at northern latitudes, and is consistent with the observed differentiation in flowering time between populations from southern and central France within the range core (Fig. 2). Taken together with the short generation time of Dittrichia, the observed genetic variation in flowering time, high heritability of flowering phenology, and selection for earlier flowering allow for the rapid evolution of phenology during the range expansion.

Caveats and open questions
Although our results strongly support the hypothesis that Dittrichia evolved to flower earlier in the year over the course of range expansion, several of our results are surprising in light of the observed differentiation in flowering time between northern and southern populations.
First, plants of a given origin tended to flower one week later in the northern range edge garden in Wageningen than in Zurich (Fig. 2). This may reflect non-adaptive plasticity

Accepted Article
This article is protected by copyright. All rights reserved. (Ghalambor et al., 2015;Ghalambor, McKay, Carroll, & Reznick, 2007), where growing conditions in the Netherlands delay the onset of flowering despite the adaptive benefits of flowering earlier at the northern range edge. Second, the strength of selection on flowering time was similar in both common gardens (Fig. 3), even though stronger selection for earlier phenology might be expected in the north based on the evolution of earlier flowering in these locations. This result may be an artefact of the specific year and location of study. Although 2016 weather conditions in Zurich and Wageningen did not deviate markedly from long term average conditions in these locations (Appendix S1, Fig. S2), it may be that earlier flowering time in the north is selected for by infrequent weather events not experienced during the year of the study. Indeed, for plants with an annual life history and limited seed bank, one year of failed recruitment due to poorly timed reproduction can be catastrophic.
Third, given the evidence for the evolution of earlier average flowering time in the Dutch origin populations, and the implication that this occurred in response to northward range expansion, it is surprising that these plants overall had similar, rather than higher fitness than central French and Swiss plants at the common garden in Wageningen (Fig. 4).
However, the fitness of range edge Dutch populations could have been reduced by other evolutionary forces that change mean fitness during range expansion, including life history trade-offs and the accumulation of deleterious mutations over the course of spread. The latter, known as the expansion load (Peischl, Dupanloup, Kirkpatrick, & Excoffier, 2013) is supported by the lower genetic diversity in Dutch than in range core populations. We therefore believe that the evolution of earlier flowering enhanced the fitness of northern populations relative to the lower fitness they would have otherwise had due to expansion load and other factors. Finally, the overall higher fitness we found in Wageningen compared to Zurich (Fig. 4) is surprising given its more northern location and the later flowering time of the plants in this garden. The most likely explanation for this result is that experimental growing conditions in each garden were influenced by microclimate factors unrepresentative
To evaluate whether late-flowering Dittrichia in our study had a higher reproductive potential than their earlier flowering counterparts, we examined the total number of reproductive structures plants produced, regardless of whether those seeds matured. At the garden in the Netherlands, the late-flowering plants originating from southern France produced 40% fewer heads than the plants from other populations (Appendix S5), indicating that earlier flowering does not seem to come at the expense of total reproduction at the northern range edge. In the Swiss garden, plants from southern France made as many heads as their more northern counterparts (Appendix S5), indicating equal reproductive potential. Although these results suggest that late flowering may be less detrimental to reproductive potential at the more southern latitude of Switzerland, a third common garden in southern France would be necessary to evaluate whether late flowering may in fact be advantageous under the warm climate in the Mediterranean core of the range.

Accepted Article
This article is protected by copyright. All rights reserved.
In combination with recent findings of rapid evolution during the native range expansion of insects (Buckley & Bridle, 2014;Lancaster et al., 2015) and birds (Gunnarsson, Sutherland, Alves, Potts, & Gill, 2011), we conclude that rapid adaptation to novel environments is not only a feature of biological invasions, but may also promote native range expansions with climate change. Though our study focused on adaptation to seasonality, a factor predictably changing with northward expansion, rapid evolution may also promote the expansion of species facing other types of novel environmental conditions. On the biotic side, range expanding populations may experience novel interactions with enemies (Doorduin & Vrieling, 2011) or mutualists. On the abiotic side, resources that were present in the native range may not be available in the expanded range, such as certain soil types or host plants (Buckley & Bridle, 2014).
Our results support growing calls for models predicting the eventual range limits and spread velocity of climate change migrants (e.g., Kearney, Porter, Williams, Ritchie, & Hoffmann, 2009) to incorporate evolutionary change (Urban et al. 2016). Given that both spreading and adapting (cf. Berg et al., 2010) may ultimately prove the most effective response of native species threatened by global warming, models properly accounting for rapid adaptation during range expansion may predict greater species persistence (Bush et al., 2016) than is typical from more traditional modeling approaches (Thomas et al., 2004). Ultimately, understanding the feedbacks between the ecology of spreading populations and their evolution in response to novel environments is key to forecasting species' responses to global change.

Accepted Article
This article is protected by copyright. All rights reserved. with an asterisk (*) were not included in the common garden experiments, because another population was sampled within the same town. In such cases, we excluded roadside populations and kept the population located within the town. Population la Roquebrussanne ( †) was not included in the common garden experiments due to limited seed availability.