Evolution of flowering time in a selfing annual plant: Roles of adaptation and genetic drift

Abstract Resurrection studies are a useful tool to measure how phenotypic traits have changed in populations through time. If these trait modifications correlate with the environmental changes that occurred during the time period, it suggests that the phenotypic changes could be a response to selection. Selfing, through its reduction of effective size, could challenge the ability of a population to adapt to environmental changes. Here, we used a resurrection study to test for adaptation in a selfing population of Medicago truncatula, by comparing the genetic composition and flowering times across 22 generations. We found evidence for evolution toward earlier flowering times by about two days and a peculiar genetic structure, typical of highly selfing populations, where some multilocus genotypes (MLGs) are persistent through time. We used the change in frequency of the MLGs through time as a multilocus fitness measure and built a selection gradient that suggests evolution toward earlier flowering times. Yet, a simulation model revealed that the observed change in flowering time could be explained by drift alone, provided the effective size of the population is small enough (<150). These analyses suffer from the difficulty to estimate the effective size in a highly selfing population, where effective recombination is severely reduced.

common conditions (see box 1 in Franks et al., 2014) or stratified propagule banks (Orsini et al., 2013), are powerful tools to reconstruct the evolutionary dynamics of populations that have faced environmental changes. Yet, observing a genetic change through time is not sufficient to claim that it is adaptation. Testing for selection as opposed to drift is one of the essential criteria for demonstrating adaptive responses, but is often overlooked (e.g., overlooked in 34% of the 44 reviewed studies based on phenotypic variation reviewed by Hansen et al., 2012). Demonstrating the influence of selection on a phenotypic change can be achieved by one of four methods (detailed in table 2 in Hansen et al., 2012;Merilä & Hendry, 2014): reciprocal transplants (Blanquart et al., 2013), Q ST -F ST comparisons (Le Corre Rhoné et al., 2010), genotypic selection estimates (Morrissey et al., 2012;Wilson et al., 2010), or tests of neutrality (pattern or rate tests, Lande, 1977). These methods all rely on measuring quantitative traits (fitness traits or traits supposed to be under selection) but require specific experimental settings. Pattern tests of neutrality rely on comparing evolution across replicates, for example, by comparing phenotypic or allele frequency changes across replicates in experimental populations, or across natural populations, assuming that they are independent replicates of the evolutionary process (same effective size and selective pressure, no migration). Pattern tests can also apply through time if a long sequence of observations is available (Sheets & Mitchell, 2001).
Alternatively, rate tests can be useful to examine the rate of genetic change in a population and compare it to the expectation under a neutral model with a given effective population size (Lande, 1976).
The effective population size (thereafter N e ) is defined as the size of an ideal Wright-Fisher population experiencing the same rate of genetic drift as the population under consideration (Crow & Kimura, 1970). Unlike experimental populations, where N e can be monitored, an accurate estimate of N e is required to perform such neutrality tests in natural populations. Temporal changes in allele frequency at neutral loci can be used to infer the effective size of the population considered (Nei & Tajima, 1981;Waples, 1989).
The ability for a population to adapt to environmental changes depends on several factors such as genetic variability, generation time, population size, or mating patterns, in particular selffertilization rates. In plants, a large fraction (40%) of species do, at least partially, reproduce through selfing (Goodwillie et al., 2005;Igic & Kohn, 2006). Selfing could challenge the process of adaptation because it directly decreases the effective population size N e (reduced number of independent gametes sampled for reproduction (Pollak, 1987); increased homozygosity; reduced efficacy of recombination (Nordborg, 2000); and increased hitchhiking and background selection (Gordo & Charlesworth, 2001;Hedrick, 1980). It is therefore expected that genetic variability is reduced in selfing populations, and empirical measures of diversity from molecular markers strongly support this prediction (Barrett & Husband, 1990;Glémin et al., 2006;Hamrick & Godt, 1996). Furthermore, several theoretical models also predict that selfing reduces quantitative genetic variation within populations (Abu Awad & Roze, 2018; Charlesworth & Charlesworth, 1995;Lande & Porcher, 2015), which has been recently confirmed by a meta-analysis of empirical data (Clo et al., 2019).
We can expect that this depleted genetic variation in predominantly selfing populations will limit their ability to adapt to changing environmental conditions and their long-term persistence and different theoretical models support this prediction (Glémin & Ronfort, 2013;Hartfield & Glémin, 2016;Kamran-Disfani & Agrawal, 2014).
Yet, empirical data examining the response of predominantly selfing populations to environmental changes remain scarce, especially for data showing short-term adaptation in the face of climate change (Qian et al., 2020). In a recent review focusing on evolutionary and plastic responses to climate change in plants, Franks et al. (2014) reported "at least some evidence for evolutionary response to climate change […] in all of these studies," and six of these 31 studies considered selfing populations.
Because there is no consensus between theoretical predictions, empirical, and experimental data, the ability of selfing populations to adapt to environmental changes remains an open question. This calls for further fine-scale population genetics analyses, with a focus on the evolutionary mechanisms involved and on the dynamics of adaptation. Here, we present a temporal survey in the barrel medic (Medicago truncatula) that enabled us to perform a resurrection study. M. truncatula is annual, diploid, predominantly self-fertilizing (>95% selfing, Bonnin et al., 2001;Siol et al., 2008) and has a circum-Mediterranean distribution. Flowering time is a major heritable trait (broad-sense heritability >0.5, Bonnin et al., 1997) that synchronizes the initiation of reproduction with favorable environmental conditions and could play a role in the adaptation to climate change. In M. truncatula, flowering time is highly variable along the distribution range and within some populations (Bonnin et al., 1997). It is mainly controlled by two environmental cues: photoperiod and temperature (Hecht et al., 2005;Pierre et al., 2008). In the Mediterranean region, there has been a significant increase in temperatures between the 80s and nowadays accompanied by a decrease in mean precipitations (http://www. world clim.org/). Most studies about adaptation in M. truncatula have so far relied on large collections of individuals representing the whole species with the aim of detecting selection footprints in the genome linked with flowering time (Burgarella et al., 2016;De Mita et al., 2011) or climatic gradients (Yoder et al., 2014).
However, the complex population structure observed at the species level can make it difficult to understand the selective history of those genes (De Mita et al., 2007). Indeed, natural populations of M. truncatula are composed of a set of highly differentiated genotypes that co-occur at variable frequencies (Bonnin et al., 2001;Loridon et al., 2013;Siol et al., 2008), a genetic structure typical for predominantly selfing species. How does this peculiar genetic composition constrain adaptation to changing environments remains unclear, but preliminary results in M. truncatula have shown that surveying the multilocus genotypic composition through time could reveal a large variance in the relative contributions of these genotypes to the next generations (Siol et al., 2007). Here, we examined the temporal change of flowering time at the population level across 22 generations characterized by changing environmental conditions (temperature and rainfall). We describe the peculiar genetic structure of this highly selfing species and investigate the genetic mechanisms involved in adaptation. In particular, we test for the role of selection as opposed to genetic drift, following four steps. First, we consider the direction of the change in trait value in relation to the environmental change. Second, we estimate the extent of genotypic selection (Morrissey et al., 2012;Wilson et al., 2010) using selection gradients for flowering time based on several fitness estimates (including an estimate of the realized fitness based on changes in frequency of the multilocus genotypes through time). Then, we estimate the effective population size, test the rate of evolution for neutrality by simulating how the frequency of the multilocus genotypes would change under genetic drift alone, and explore the effect of the imprecision in the estimation of effective size. Finally, we examine the change in flowering time during the same time period at the regional scale, using one individual per population across the distribution range of M. truncatula in Corsica. A similar genetic change at the regional scale would lend weight to the hypothesis that the change in flowering time occurred in response to selection.

| Studied population and experimental design
The focus population (F20089 or CO3 according to Jullien et al.,  were successfully multiplied. Out of these, 55 families for each of the two sampling years were randomly chosen in 2012. Seeds from the 110 families were scarified to ease germination and were transferred in Petri dishes with water at room temperature for six hours. We then used two different vernalization treatments (at 5°C during 7 or 14 days) to compare the vernalization requirement between the two years. Five replicates from each vernalization treatment were transferred back to the greenhouse, according to a randomized block design (five blocks and two treatments, adding up to a total of ten replicates per family, 1100 plants in total). Data loggers were placed on each table to monitor temperature and humidity. For each individual, the number of days after germination to form the first flower was recorded. In addition, the total number of seeds produced by each plant throughout its lifetime was measured as a proxy for fitness.

| Temporal changes in flowering time
Individual flowering times were converted to thermal times following Bonhomme (2000). The thermal time was calculated as the sum of the mean daily effective temperatures of each day between sowing and the emergence of the first flower, where the mean daily effective temperature is the day's mean temperature minus the base temperature (T b ). We used T b = 5°C, as reported by Moreau et al. (2007) for the M. truncatula reference line A17.
Plants noted as sick or failing to produce leaves were removed from the datasets (22 individuals removed). Collected measures were tested for normality using quantile-quantile (Q-Q) plots (Nobre & Singer, 2007). All analyses were conducted using R version 2.15.2. We used linear mixed models (lme4 package) to test This maximal model was simplified, using likelihood ratio tests (LRT) to compare the models. In addition, we tested for a significant change in genetic variance between 1987 and 2009 using a LRT between the model [1] and a model where family is not nested into year. Standard errors for variance components were estimated using jackknife resampling. We used the variance components estimated for the random effects to calculate broad-sense heritability as H 2 = V G ∕V P , where V G is the genetic variance as estimated by the family effect and V P is the total phenotypic variance, including block, family, and residual variance. Standard errors for H 2 were estimated through jackknife resampling on families (Sokal & Rohlf, 1995).

| Temporal changes in sensitivity to vernalization
Selection on a trait in an environment can shift both the mean and the plasticity of that trait. Here, we considered the sensitivity to vernalization cues, measured as the slope of the regression line (1) between individual values and the environmental value (estimated as the average phenotype, Y) (Falconer & Mackay, 1996), for each individual i: For each family, the five individuals in each treatment were paired according to their position in the greenhouse (block 1 with block 5, etc.). This coefficient assumes that reaction norms are linear (Gavrilets & Scheiner, 1993;Scheiner, 1993) and this approximation is expected to work well (Chevin et al., 2013). We used a linear mixed model, with sampling year (1987 or 2009) as a fixed effect, a random block effect and its interaction with year, and a family effect (genetic effect) nested into year. As for flowering time, we estimated the broad-sense heritability of the vernalization sensitivity.
A genetic correlation between flowering time and sensitivity to vernalization would affect the response to selection in the context of climate change. We therefore used a bivariate model with the sensitivity to vernalization and the flowering time measured in the short vernalization treatment as two dependent variables to estimate their genetic covariance with a random family effect, including block as a random effect, using AsReml (Gilmore et al., 2009). We ran an independent model for each sampling year. The significance of genetic covariances was tested by comparing the residual deviance of the final model with that of a model with a fixed covariance of zero in a log-likelihood ratio test.

| Selection gradient for flowering date: genetic covariance analysis
In the absence of selection for the trait considered, its observed variation is expected to be independent from fitness. We tested this by measuring the selection gradient, that is, the statistical relationship between a trait and the fitness. Selection gradients were established for each year (and per treatment) following the Robertson-Price identity that states that ΔZ, the expected evolutionary change in the mean phenotypic trait z per generation is equal to Θ a (z, w), the additive genetic covariance of the trait z, and the relative fitness w (Price, 1970;Robertson, 1966): Here, we estimated the broad-sense genetic covariance Θ g . Assuming that the dominance variance is negligible due to the very high levels of homozygosity in selfing populations (Holland et al., 2010), genetic covariance should be a good approximation of the additive genetic covariance (we neglect maternal genetic effects here). As a preliminary step, we checked whether our proxy for fitness, the relative seed production, had significant genetic variance. The relative seed production was measured as the individual seed production standardized by the average seed production of individuals from the same year and treatment. A mixed model was used to analyze the relative seed production, including two random effects for block and family.
Then, we provided there was significant genetic variance for relative seed production in the population each year, and we analyzed it in a bivariate model with flowering time to estimate the genetic covariance with a random family effect, including block as a random effect, using AsReml (Gilmore et al., 2009). Again, the significance of genetic covariances was estimated by comparing the residual deviance of the final model with that of a model with a fixed covariance of zero in a log-likelihood ratio test.

| Genetic analyses
During the multiplication generation in the greenhouse (2011), 200 mg of leaves was sampled from each plant for DNA extraction, using DNeasy Plant Mini Kit (Qiagen). Twenty microsatellite loci were used for genotyping (see the details of amplification reactions and analyses of amplified products in Jullien et al., 2019;Siol et al., 2007). Briefly, samples were prepared by adding 3 μl of diluted PCR products to 16.5 μl of ultrapure water and 0.5 μl of the size marker AMM524. Amplified products were analyzed on an ABI prism 3130 Genetic Analyzer, and genotype reading was performed using GeneMapper Software version 5.

| Single-locus analyses assuming independence among loci
As a preliminary step, the data were filtered to reduce the percentage of missing data (loci or individuals with >10% missing data were removed), and to discard monomorphic loci. After filtering, the dataset comprised 145 individuals (representing 145 families) and 16 loci (64 from the year 1987 and 81 from the year 2009). We measured the genetic diversity of the population each year using the allelic richness N a − rar (Hurlbert, 1971) and the expected heterozygosity H e . In this predominantly selfing population, we expect a strong deviation from Hardy-Weinberg heterozygosity expectations. Thus, for each sampling year, we estimated the inbreeding fixation coefficient F IS and its confidence interval using 5000 bootstraps over loci. Between year differences for N a − rar , H e and F IS across loci were tested using Wilcoxon signed-rank tests. Analyses were performed in R using the packages adegenet (Jombart, 2008) and hierfstat (Goudet, 2005) and the program ADZE for rarefaction analyses (Szpiech et al., 2008). The percentage of pairs of loci showing significant linkage disequilibrium (LD) was calculated using Genepop (Rousset, 2008) with a threshold of 0.05. Finally, we measured the temporal variance in allele frequencies using the F ST estimator by Weir and Cockerham (1984). To estimate the effective population size (N e , measured in number of diploid individuals) from the temporal variance of allele frequencies, we used F ST estimates to account for the correlation of alleles identity within individuals due to selfing (Navascués et al., 2020) and followed the method outlined in Frachon et al. (2017). We measured a confidence interval for N e using an approximate boot-

| Analyses based on multilocus genotypes
We used the program RMES to estimate selfing rates from the distribution of multilocus heterozygosity (David et al., 2007). We tested for a difference in selfing rates between years using a likelihood ratio test between models where the selfing rate was constrained to be constant or not. For each sample (1987 and 2009), we examined the genetic structure by sorting out the number of multilocus genotypes (thereafter called MLG) and measuring their frequency and redundancy through time using GENETHAPLO (available on GitHub at https://github.com/lauga y/Genet Haplo and described in Appendix S1: Section S1). GENETHAPLO takes into account the uncertainty of the assignment of a genotype to a MLG group due to missing data: In case of ambiguity, an individual is randomly assigned to one of the candidate MLG group with a probability proportional to the MLG group size. The approach also considers a genotyping error rate: If two individuals differ by less than the error rate, they are considered to belong to the same MLG. After an initial run with an error rate of zero, we checked the distribution of the distances between MLGs. We found an excess of small distances, which could indicate errors in genotype assignation (Arnaud-Haond & Belkhir, 2007). We corrected this by re-running the program with an error rate of 1/16 (= one mis-read locus). GENETHAPLO also searches for residual heterozygosity (defined as the proportion of heterozygous loci in the multilocus genotype) and evidence for recombination (S1). To identify putative recombination events between MLGs, it uses the genetic distances: a MLG is a recombinant candidate if the sum of its allele differences with two other MLGs ("parental MLGs") equals the number of allele differences between these two parental MLGs.
If a MLG has a high fitness in a given environment, plants carrying this MLG will produce on average a larger progeny and the frequency of the MLG will rise in the following generations. We therefore propose to use the absolute change in frequency of the fully homozygous MLGs through time as an indicator of their "realized fitness." As a preliminary step, we checked whether selection quantified in the greenhouse is likely to mirror the predominant selection between 1987 and 2009 using a linear model to verify whether the change in MLG frequencies covaries positively with and can be predicted by the seed production in the greenhouse. We For each of these models, we verified the normality of the residuals and estimated a confidence interval for the slope using profile likelihood confidence bounds.
In addition, we tested whether the change in frequencies of the MLGs reflects a response to selection or can be expected by drift alone. This was tested by simulating the effect of 22 generations of drift, using an extension to multi-allelic data of the approach described in Frachon et al. (2017) and inspired by Goldringer and Bataillon (2004). Again, only the fully homozygous MLGs were kept for this analysis. We assumed complete selfing during the time interval so the whole genome behaves as a single super-locus. Details about the procedure used to simulate individual MLG frequency trajectories are provided in Appendix S1: Section S2. We simulated

| Regional analysis
Finally, we attempted to disentangle selection and drift by considering other populations located in the same geographic region as the focal population and therefore likely submitted to the same selective pressure due to climatic constraints (pattern test, as described in the Section 1). For this regional analysis, we used 16 populations of M. truncatula across Corsica that were sampled twice, once in the 80s and again in the early 2000s (listed in Table S1). Samples consisted of around 100 pods collected along transects running across the populations. Seeds collected were stored in a cold room. In 2010, one pod randomly selected from each sample (80s and 2000s) was threshed and one plant per population per year was replicated through selfing in standardized greenhouse conditions. This greenhouse generation allowed suppressing potential maternal effects (as in the experiment with the Cape Corsica population) and resulted in 32 families (16 populations × 2 years) of full-sibs produced by selfing. In 2011, seeds from the 32 families were germinated following the same protocol as described earlier for the intrapopulation analysis, but with only one vernalization treatment at 5°C during seven days. Five plants for each family were then transferred to tables in the greenhouse according to a randomized block design (five blocks). We monitored the temperature and humidity and the flowering time for each plant.
Individual flowering times were converted in thermal time, in the same way as it was done for the intrapopulation analysis. Again, we used linear mixed models (lme4 package) to test for the effect of sampling year on flowering time. The model included a single fixed effect for the sampling year (1980s or 2000s). The block effect was included as a random effect, along with its interaction with sampling year. A random population effect was also included and replaced the "family" effect of Equation (1) seen as there was only one family per year in this regional sample. The resulting model was written as: Again, this maximal model was simplified using likelihood ratio tests.

| Changes in flowering time
Visual inspection of the Q-Q plots indicated that the residuals from all the linear models we used were normally distributed. We found that flowering time differed significantly between years: plants sampled in 2009 flowered on average over two days earlier than plants sampled in 1987 (Table 1, Figure 1). This effect remained significant when we analyzed flowering time as a number of days rather than degree.days (results not shown). Longer vernalization also sped flowering up (treatment effect, Table 1). The block effect only explained a low proportion of variance (micro-environment) and the largest variance component was the family effect, for all combinations of years and treatments.
The comparison of a model where family was nested in years only or in years × treatments showed that the family × treatment interaction was significant (χ 2 = 66.1; df = 7; p = 9 × 10 −12 ). It means that the reaction norms for the different genotypes were not parallel (Figure 1), because the genotypes responded differently when exposed for a shorter period to cold temperatures. To account for this genotype × environment interaction, the heritability for flowering time was estimated in each vernalization treatment separately (four components of variance, Table 1). It varied between 0.53 and 0.77 (Table 2) We found no significant year effect on the sensitivity to vernalization (χ 2 = 1.7; df = 1; p = .185). There was no significant difference in the family effect between years (interaction family × year not significant; LRT: χ 2 = 1.2; df = 2; p = .552), but the family effect was highly significant (χ 2 = 32.6; df = 1; p = 1 × 10 −8 , Table S2) and the heritability of the sensitivity to vernalization was 0.19 (±0.04) ( Table 2)

| Selection gradient for flowering date
The relative seed production showed significant genetic variance (family effect, Table S3, heritability of 0.34, Table 2), which enabled us to build multivariate models to examine selection gradients following Equation (2). In 1987, we found a significant genetic covariance between flowering time and relative fitness: Θ a (z, w) = −20.5; (3)     Figure   S1). We found no evidence for a link in terms of recombination or seg-

TA B L E 2
Heritabilities (H 2 ) and coefficients of genetic variance (CV g ) for flowering time in each vernalization treatment (T1: short vernalization; T2: long vernalization) and each sampling year, for sensitivity to vernalization, and for relative seed production the short vernalization treatment (n = 12; Figure 3b). In addition, the negative slope was mostly supported by the decreasing frequency of the two late flowering MLGs that were prevalent in 1987. The simula-  Figure 4).

| Changes in flowering time at the regional level
At the regional level (Equation 3), we found no effect of the interaction between block and sampling year (LRT χ 2 = 0; df = 1; p = 1).
All other effects were significant (

| DISCUSS ION
Pairing up a resurrection study with population genetic analyses

| Can we use effective population size estimates to test whether the genetic change is caused by selection or drift in a predominantly selfing population?
As pointed out in the Introduction, simulating drift is one of the methods to test whether selection has occurred, but it requires knowledge about the effective population size. Using changes in allele frequencies between 1987 and 2009 in a natural population, we estimated a temporal F ST of 22.6%, which corresponds to an effective size of 19 (95% confidence interval: 15-25). This estimate is several orders of magnitude lower than the census population size (>2000 individuals) and lower than expected given the observed levels of diversity (Nordborg & Donnelly, 1997). Similarly, low effective population sizes have been estimated previously in other M. truncatula populations, based on the temporal variance in allele frequencies (Siol et al., 2007), and attributed to the high selfing rate of this species. Yet, the observed levels of polymorphism are often incompatible with such drastically low effective sizes (see figure 3c in Hereford, 2009;Jullien et al., 2019). N e estimates are likely biased and/or imprecise, because some of the assumptions underlying the temporal method are violated, for example, isolation of the populations under scrutiny, absence of selection, and independence of marker loci (Jullien et al., 2019). For example, the quick change in allele frequency caused by a migration event will be misinterpreted as strong drift because temporal methods estimate N e using the pace at which allele frequency changes and therefore underestimate it (Wang & Whitlock, 2003). In addition, strong selfing affects the precision of temporal F ST estimates because the number of independent loci is reduced (Appendix S1: Section S3).
In our focal population, the whole genome behaves practically as a single locus, which limits the precision of our effective size estimates. Unfortunately, we show in Appendix S1: Section S3 that inferring effective size from the variation of MLG frequencies (i.e., considering a single, multi-allelic superlocus) is unlikely to improve the quality of our estimates.
Finally, if selection occurs in a nonrandom mating population, it will exacerbate the Hill-Robertson effect and further reduce the effective size (Comeron et al., 2007). Indeed, selection will create heritable variance in fitness among individuals, thereby locally reducing N e (Barton, 1995;Charlesworth & Willis, 2009;Robertson, 1961). In predominantly selfing species, due to drastically reduced effective recombination (Nordborg, 2000), selection will extend the reduction in diversity caused by the selective sweep to a larger proportion of the genome compared to a random mating population (Caballero & Santiago, 1995;Kamran-Disfani & Agrawal, 2014). With selection, the effective size estimated using the temporal variance in allele frequencies can therefore not be considered as a "neutral" effective size but rather reflects the combined effects of inbreeding and se- If highly selfing organisms strongly deviate from the general assumptions of population genetics models, a major benefit, however, is that the temporal survey of MLGs provides a highly integrative measure of fitness, which is analogous to measures of genotypespecific growth rates in asexual organisms. Our results show that changes in frequencies of MLGs through time are positively correlated with the fitness measured as the seed production in the greenhouse (Figure 3a). This relationship is not significant if we consider all the MLGs found in 2009, but this is not surprising considering the potentially strong environmental variance in the field and the TA B L E 3 Effect of sampling year on flowering time at the regional scale, taking into account the effect of the population of origin of each line. The effect on the mean flowering time is given for the fixed year effect and variance components are given for random effects (with standard errors in brackets). For each component, the degrees of freedom, likelihood ratio (χ 2 ), and p-values are reported  (Shaw et al., 2008), which are generally hindered by pervasive trade-offs between life history traits such as reproduction and survival (Ågren et al., 2013;Anderson et al., 2014).

| What selective pressure could have led to this genetic change in flowering time? Insights from ecophysiology
The evidence that the change in phenology observed in this population across 22 generations is the result of selection as opposed to drift remains equivocal. A further step toward evaluating whether selection is responsible for the genetic change observed is to characterize the potential selective pressure involved. Phenological changes associated with climate change have been reported in a large number of plants (Amano et al., 2010;Cleland et al., 2007;Parmesan & Yohe, 2003;Root et al., 2003). In this context, ecophysiological models of phenology are insightful to understand how climate change can affect traits such as flowering time (Chuine, 2000;Oddou-Muratorio & Davi, 2014). The phenological response to climate change is complex, because the promoting effect of increased temperatures opposes the influence of reduced vernalization (Wilczek et al., 2010). Ecophysiological models generally predict a plastic shift toward earlier flowering times, as long as vernalization is sufficient during winter (Morin et al., 2009). In agreement with these predictions, a meta-analysis exploring the phenological response to climate change in plant populations showed that phenotypic changes are mostly plastic, while evidence for genetic adaptation remains relatively scarce (Merilä & Hendry, 2014, and other references of Evolutionary Applications special issue, January 2014). However, a large part of the intraspecific variation observed in phenology is genetic (Hendry & Day, 2005) and the architecture of the network underlying flowering time variation is well described in some species such as Arabidopsis thaliana (Sasaki et al., 2018;Wilczek et al., 2010). How climate change will affect the genetic values of phenological traits remains uncertain. In a first hypothesis, we may assume that the phenotypic optimum for flowering time is not affected by climate change. We therefore expect a genetic change occurring in the opposite direction than that of the plastic response ( Figure 5a).
This hypothesis resembles counter-gradient variation, which occurs when the genetic influence on a trait along a gradient opposes the environmental influence, resulting in reduced phenotypic variation across the gradient (Levins, 1969). Counter-gradients are widespread along geographic gradients, as shown by the meta-analysis by Conover et al. (2009), who found evidence for counter-gradient in 60 species and for cogradients in 11 species. Therefore, assuming that the same mechanism observed across spatial gradients could occur in temporal gradients, we would expect the genetic response of flowering time to counterbalance the plastic response to climate change. This could be achieved for example with a genetic change increasing the base temperature T b (temperature below which the development is supposed to be nil).
Yet, our temporal survey rejects the countergradient hypothesis, both at the population and at the regional scale. Instead, we found . We can therefore expect that the earliest a plant flowers, the highest its fitness. Second, climate change in the Mediterranean region also tends to reduce precipitations in spring and early summer (Goubanova & Li, 2007;Schröter et al., 2005), thereby shortening the reproductive period.
Severe early summer drought could therefore create a strong selective pressure toward earlier flowering. Such a genetic shift in flowering time in response to extended drought has been reported before in the literature (Franks et al., 2007). In terms of ecophysiology, it can be caused by lower requirements of degree.days, or a reduction of the base temperature T b .
Finally, although it is generally assumed that flowering date should be under stabilizing selection in order to avoid frost or drought when flowering occurs, respectively, too early or too late, a recent meta-analysis found widespread evidence for frequent directional selection toward early flowering (Munguía-Rosas et al., 2011).
Selection estimates considered in this meta-analysis largely ignore the effect of variation in number of flowers and plant size, which could bias the results. Yet, it remains that early flowering could have several advantages, among which an increased time for seed maturation in early reproducing plants and a longer period of growth for the progeny issued from seeds that germinate immediately (as reviewed by Elzinga et al., 2007;Kudo, 2006). Under this scenario of directional selection, we also expect a pattern of cogradient, as observed in the data (Figure 5c).
Besides the evidence for a genetic change in flowering date in M. truncatula in Corsica, we found no evidence for a change in the sensitivity to vernalization, despite genetic variance for this trait in the population (H 2 = 0.19). In the literature, most studies have found at least some genetic variation for plasticity, but corresponding heritabilities were generally low (Scheiner, 1993). Our results also suggest that the sensitivity to vernalization is not independent from flowering date, because the intercept and the slope of the reaction norm to the vernalization treatment are genetically correlated (Gavrilets & Scheiner, 1993

CO N FLI C T O F I NTE R E S T
The authors of this article declare that they have no financial conflict of interest with the content of this article.

This article has been awarded Open Data and Open Materials
Badges. All materials and data are publicly accessible via the Open Haplo.

DATA AVA I L A B I L I T Y S TAT E M E N T
Phenotypic data for the intrapopulation and interpopulation experi- year index of first flowering dates and its response to temperature