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

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.

INTRODUCTION

When facing changing environments, organisms can persist by one of three strategies: fleeing (migration), coping (plasticity), or adapting. If migration and plasticity can lead to rapid and reversible changes in the average phenotype of a population, adaptation proceeds through genetic changes and toward phenotypes with the highest fitness in a given environment. The literature describing adaptation in natural populations is vast (Bay et al., 2017; Côté & Reynolds, 2012; Kremer et al., 2012; Olson‐Manning et al., 2012), and the recent rise of next generation sequencing has enabled tremendous progress in our knowledge about the genetic architecture of adaptation at the species level (Barrick & Lenski, 2013; Brown, 2012; Fournier‐Level et al., 2011; Jones et al., 2012).

Long‐term temporal surveys (Visser, 2008) and resurrection studies, where ancestors and descendants are compared under 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 & Kremer, 2012; 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 self‐fertilization 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.worldclim.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.

MATERIALS AND METHODS

Studied population and experimental design

The focus population (F20089 or CO3 according to Jullien et al., 2019) is located in Cape Corsica (42°58.406′N–9°22.015′E). In 1987 and 2009, around 100 pods were collected along three transects running across the population, with at least one meter distance between each pod collected, in order to avoid over‐sampling the progeny of a single individual. Seeds collected in 1987 were stored in a cold room. In 2011, pods collected in 1987 and 2009 were threshed and seeds were replicated through selfing in standardized greenhouse conditions to control for maternal effects and build families of full‐sibs produced by selfing. Seeds for this generation of multiplication were randomly selected from pooled samples of seeds from 1987 and 2009. 64 families collected in 1987 and 96 in 2009 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 for a significant change in flowering time between the sampling years. The model included two fixed effects: sampling year (1987 or 2009) and treatment (short or long vernalization) as well as their interaction. Block (nested in treatment), block × year, and family were random effects. The family effect was nested in years because we were interested in estimating the genetic variance within population each year of collection. The interaction between family and treatment was included in the family effect as a vectorial random effect. The complete model is summarized in Equation (1), where Y denotes the flowering time, μ the average flowering time over the whole sample, and the residuals:

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 , where is the genetic variance as estimated by the family effect and 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 between individual values and the environmental value (estimated as the average phenotype, ) (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 , the expected evolutionary change in the mean phenotypic trait z per generation is equal to , 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 . 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 (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 , 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 bootstrap method (DiCiccio & Efron, 1996) over loci.

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/laugay/GenetHaplo 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 then measured the selection gradient for flowering time as the slope of the regression of the change in frequency of the MLGs between 1987 and 2009 with the genetic value of flowering time (measured as the average flowering time for a given MLG in the short vernalization treatment). We compared this pattern with the predictions from the Robertson‐Price selection gradient. The MLGs found in 2009 but absent in 1987 may have been undetected in 1987 due to low frequency, or may be recent migrants. Their change in frequency between 1987 and 2009 is thus necessarily positive and may not accurately reflect their realized fitness. We therefore reiterated these analyses using a dataset restricted to the MLGs present in 1987 only. 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 each generation of drift by drawing MLG counts from a multinomial distribution parameterized with the effective population size N e estimated from the temporal F ST, and the MLG frequencies in the previous generation. Note that this simulation assumes a generation time of one year and therefore neglects seed dormancy and that the presence of a seed bank would reduce the rate of genetic drift. After 22 simulated generations, we randomly sampled 75 individuals to estimate the frequencies of each MLG and measured the change in MLG frequencies across the 22 generations. This was repeated for 104 replicates in order to draw the distribution of the change in MLG frequency expected by drift alone. To account for the potentially large estimation variance for the F ST (as observed in the simulations performed in Appendix S1: Section S3), we examined the sensitivity of the analysis to the effective population size using a range of values (). Finally, we examined the simulated selection gradient as the relationship between the simulated changes in MLG frequencies through time and the genetic value of flowering time previously measured for each MLG, using a linear model. This provided us with a null distribution of the slopes of the regression between frequency change and flowering time, expected under drift alone. We then tested for the significance of the observed slope against the simulated distribution, by computing the proportion of the simulated slopes that were greater than the observed value.

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.

RESULTS

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; ). 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). The genetic variance for flowering time in the population remained the same in 1987 and 2009, as shown by a LRT between the full model (Equation 1) and a model where family was not nested in year (χ 2 = 6.65; df = 7; p = .47). 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; , Table S2) and the heritability of the sensitivity to vernalization was 0.19 (±0.04) (Table 2). Finally, the multivariate analysis highlighted a strong positive genetic correlation between flowering time (measured in the short vernalization treatment) and the sensitivity to vernalization (in 1987: 0.54 p = .008; in 2009: 0.60 p < .0001), which means that early flowering plants are less sensitive to vernalization cues. Using the flowering time measured in the long vernalization treatment, we observed the same pattern of correlation.

TABLE 1

Effect of sampling year and treatment on flowering time in the cape Corsica population, taking into account the family effect (genetic effect). Effect values on mean flowering time are given for fixed effects and variance components are given for random effects (with standard errors in brackets). The family effect was nested into year (1987 or 2009) and treatment (T1: short vernalization treatment; T2: long vernalization treatment), leading to four variance components. For each component, the degrees of freedom, likelihood ratio (χ 2), and p‐values are reported. None of the interactions considered in the complete model [1] were significant: between year and treatment (LRT χ 2 = 1.8; df = 1; ); between block and year (χ 2 = 0.0006; df = 1; )

Tested effect on flowering time	Mean effect or variance component (SE)	df	χ 2	p
Year	−28.76 a	1	7.3	.007
Treatment	−162.84	1	42.2	8×10^‐11
Block	92.34 (9.61)	1	34.5	4×10^‐9
Family|year × treatment	1987‐T1: 2807.90 (872.97)	10	850.4	2×10^‐16
	1987‐T2: 1793.51 (500.25)
	2009‐T1: 5449.80 (1200.16)
	2009‐T2: 3557.01 (1408.88)
Error	1500 (38.73)	1081

Assuming an average daily temperature of 15°C over the time period considered, the difference of 28.76 degree.days corresponds to two days.

FIGURE 1

Average flowering time per family for the two sampling years and the two vernalization treatments. Short vernalization is in gray and long vernalization in black. The large dots and the horizontal lines stand for the average flowering date for each vernalization treatment, for the years 1987 (dotted lines) or 2009 (dashed lines). Black crossing lines indicate that the reaction norms differ between families, as expected if genotype × environment interactions are significant

TABLE 2

Heritabilities (H 2) and coefficients of genetic variance (CVg) 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

Trait	H 2 (SE)		CVg
	1987	2009	1987	2009
Flowering time	T1: 0.64 (0.06)	T1: 0.77 (0.04)	5.70	8.11
	T2: 0.53 (0.07)	T2: 0.69 (0.07)	5.49	8.03
Sensitivity to vernalization	0.19 (0.04)		18.14
Relative seed production	0.34 (0.03)		30.00

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: = −20.5; LRT comparing this model with a model where the genetic covariance was constrained to be zero: χ 2 = 60.2; df = 1; . The covariance remained significantly different from zero when we used the lines derived from the sampling in 2009: = −18.5; LRT: χ 2 = 12.4; df = 1; . A similar negative relationship was observed among lines derived from each of the two years, which means that the selection gradients predict an evolution toward early flowering under the environmental conditions of the greenhouse (Figure 2).

FIGURE 2

Selection gradients for flowering time. Established as the relationship between the genetic value for flowering time (family average, in degree.days) and the genetic value for relative fitness (family average of the relative number of seeds), for each sampling year and vernalization treatment. Lines stand for the linear regression

Changes in the genetic composition of the population

The analysis of microsatellite data highlighted high levels of genetic diversity for both sampling years, with an increase between 1987 and 2009 only significant for H e (Table S4). This suggests that the increased diversity between 1987 and 2009 reveals more balanced allele frequencies rather than an increase in the average number of alleles. The temporal differentiation measured using the 16 loci was high (F ST = 0.226; 95% confidence interval: 0.182–0.269), which translates into a particularly small effective size (N e = 19 diploid individuals; 95% confidence interval: 15–25). According to equation 16 in Nordborg and Donnelly (1997), we predict that , where N e is the effective size as estimated above. Using mutation rates for dinucleotide microsatellite loci measured in Arabidopsis thaliana ( to ) (Marriage et al., 2009), and assuming an isolated population at equilibrium, we expect that H e should lie between 0.004 and 0.134, which is nearly three times lower than the H e estimated here (Table S4). The observed heterozygosity was particularly low, resulting in large F IS estimates, as expected for a predominantly selfing species. The estimated selfing rate was about 94% in 1987 and rose to 98% in 2009 (statistically significant increase, Table S5). This high selfing rate results in extensive linkage disequilibrium between loci (nearly all pairs of loci are in linkage disequilibrium, Table S4), which makes the analysis of multilocus genotypes particularly relevant.

The analysis of MLG identified 60 different MLGs in this sample of 145 individuals. Out of the 60 MLGs, 48 were fully homozygous at the 16 loci and 12 MLGs displayed some level of heterozygosity (Figure S1). We found no evidence for a link in terms of recombination or segregation between the heterozygous MLGs and any of the fully homozygous MLGs. These heterozygous MLGs were therefore excluded from the following analyses, leaving us with 48 MLGs (58 individuals in 1987 and 75 in 2009). The two predominant MLGs represented more than 50% of the population in 1987 and nearly 20% in 2009. These, as well as three other MLGs, were observed in both years (Figure S2). The absolute changes in homozygous MLGs frequencies through time tended to covary positively with the total number of seeds produced by a plant in the greenhouse (Figure 3a, regression only significant with the sample restricted to the MLGs present in 1987, n = 12 MLGs), which provides support to use it as a proxy to estimate the realized fitness. We therefore used the change in frequency of the 48 MLG (58 individuals in 1987 and 75 in 2009) to build selection gradients for flowering time. Again, we found a gradient with a negative slope (Figure 3b), suggesting that the late flowering MLGs have a reduced realized fitness compared to earlier ones. This confirms the reduced fitness of late flowering genotypes observed in our greenhouse experiment (Figure 2). Yet, the effect of flowering date on the realized fitness was small and only significant when the dataset was restricted to the MLGs present in 1987 and measured in 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 simulation of 22 years of drift with an effective population size of 19 showed that the slope of the observed selection gradient did not deviate significantly from the distribution expected by drift alone (p = .182). Yet, again, when we restricted the dataset to the MLGs that were present in 1987, the observed selection gradient deviated significantly from the distribution expected by drift alone (p = .047), which suggests that the drift‐alone hypothesis could be rejected.

FIGURE 3

Analyses of the “realized fitness,” estimated as the absolute change in frequency of the MLGs through time. MLGs with residual heterozygosity were removed from this analysis. (a) Relationship with the average number of seeds produced by plants of a given MLG in the greenhouse. (b) Selection gradient for flowering time. Each point stands for the average flowering date for a given MLG. The black regression lines are estimated using all points (n = 48; a: slope =  points of frequency per seed p = .094; b: slope = −0.0002 95% confidence interval: −0.0006; 0.0001 p = .179). This includes MLGs that were not observed in 1987 (black dots), for which the change in frequency is necessarily always positive. The dotted lines are the regression lines for the analysis restricted to the MLGs present in 1987 (white dots only; n = 12; a: slope = 0.0002 p = .024; b: slope = −0.0009 95% confidence interval: −0.0017; −0.0002 p = .038). Q–Q plots for the selection gradients are provided in Figure S3

Because selfing reduces the effective recombination, it reduces the number of independent loci. Measuring F ST from linked loci therefore amounts to measuring it from a lower number of markers, and it is known that F ST estimates based on a few loci suffer from a large sampling variance (Weir & Hill, 2002). Alternatively, we could have concatenated the genotypes at the different loci to compute a diploid version of the haplotype‐based F ST (Mehta et al., 2019). Using the changes of frequencies for 48 homozygous MLGs, we estimated a temporal F ST of 0.075, which corresponds to an estimated effective size of 136. However, our simulations (Appendix S1: Section S3) show that these haplotype‐based F ST estimates are strongly downward biased, due to the dependency of F ST with allelic diversity (Alcala & Rosenberg, 2017; Edge & Rosenberg, 2014; Jakobsson et al., 2013) and could therefore overestimate the effective population size. Instead of using this unreliable estimate of 136, we assessed the sensitivity of our neutrality test for MLG frequency changes to the effective population size estimates, using a range of values (). We found that the observed selection gradient can no longer be explained by drift alone if the effective population size exceeds 150 (or even 10 if we consider only the MLGs present in 1987, Figure 4).

FIGURE 4

Test of selection for increasing values of N e. p‐Value, defined as the proportion of simulated datasets where the slope of the selection gradient is steeper than the observed slope, for the simulations of drift alone (a) considering all the homozygous MLGs (n = 48) or (b) considering only the MLGs that were already present in 1987 (n = 12). The dotted line indicates the 0.05 threshold value for significance. The vertical dashed line is the effective size estimated using the temporal F ST and considering the 16 microsatellite loci as independent (N e = 19; p = .182 with n = 48 (a); p = .047 with n = 12 (b))

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 (Table 3): The random block effect only explained 5% of the total variance, whereas the population effect accounted for 34% of variance. The significant year effect showed that the material we collected in 2005 or 2009 in Corsica flowered about five days earlier (78 degree.days, Table 3) compared to the one we collected between 1987 and 1990.

TABLE 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

Tested effect on flowering time	Mean effect or variance component (SE)	df	χ 2	p
Year	−78.00 a	1	9.3	.002
Block	2379 (1029)	1	5.7	.017
Line	14,874 (4423)	1	40.1	2×10^‐10
Error	26,971 (8260)	167
Total variance	44,224

Assuming an average daily temperature of 15°C over the time period considered, the difference of 78.00 degree.days corresponds to five days.

DISCUSSION

Pairing up a resurrection study with population genetic analyses proved highly insightful to understand how flowering time changed through time in M. truncatula and to get insights into the mechanisms involved. Growing plants collected in the Cape Corsica population 22 generations apart in a common garden experiment provided evidence for a diminution of flowering times by about two days (i.e., a reduction between 2 and 4% in flowering time). This study also highlighted the peculiar genetic structure of this highly selfing population, where some multilocus genotypes are persistent through time. This enabled us to measure the fitness of a genotype as its frequency change through time and to establish a multilocus selection gradient. We used this multilocus fitness measure as well as a fitness measure based on individual seed production in the greenhouse to estimate the selection gradient for flowering time. Both gradients predict evolution toward earlier flowering but only the selection gradient using seed production in the greenhouse as a proxy for fitness was significant. It should be kept in mind that the selection gradient could change if the plants were growing in their natural environment, due to potential Genotype x Environment interactions. Simulating evolution across 22 generations showed that the observed change in flowering time can be caused by drift alone, providing the effective size of the population is lower than 150. These analyses suffer from the difficulty to estimate the effective size in a highly selfing population, where effective recombination is severely reduced.

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 selection (Le Rouzic et al., 2015). Overall, due to the reduced effective recombination and potential migration, predominantly selfing populations can strongly deviate from the assumptions of the temporal method to estimate effective size and such estimates should be treated with caution (see figure 3 in Jullien et al., 2019).

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 genotype‐specific 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 approximation due to the possibility that a MLG that was absent in 1987 appeared within the 22 years of time period. A larger sample size in 1987 or additional temporal samples could help improve this analysis. Despite these imprecision, such integrative estimates of fitness are highly valuable because of the difficulty to obtain lifetime measures of fitness in the field (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).

FIGURE 5

Hypotheses for the expected selective pressure on flowering time under climate change. (a) Selective response expected under the hypothesis that the phenotypic optimum for flowering date remains the same. The selective response is expected in the opposite direction compared to the plastic response to increased temperatures. This corresponds to the countergradient hypothesis. (b) Selective response expected under the hypothesis that the phenotypic optimum for flowering date is displaced with climate change and that it becomes advantageous to flower earlier. The selective response is expected in the same direction as the plastic response to increased temperatures. This corresponds to the cogradient hypothesis. (c) Selective response expected under the hypothesis that flowering time is under directional selection

Yet, our temporal survey rejects the countergradient hypothesis, both at the population and at the regional scale. Instead, we found evidence for a genetic change toward earlier flowering, in the same direction as the plastic response to the environmental change (here a rise in temperatures). Such a co‐gradient is expected if climate change has shifted the phenotypic optimum toward earlier flowering dates (Figure 5b). Several hypotheses could explain such a shift and the resulting cogradient. First, in a plant with undetermined flowering such as M. truncatula, reduced frost risk early in the season should favor earlier flowering, because plants that manage to flower early in the season will carry on producing flowers until summer drought becomes limiting (end of May–June). 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). Therefore, a lower number of chilling units received during winter (short vernalization treatment) result in higher heritability of flowering date. This correlation could favor the selective response of flowering date to climate warming because warmer winters will inflate the genetic variance of flowering date. Alternatively, if early flowering genotypes are selected for, or arrive in the population by migration, the evolution of the sensitivity to vernalization might be constrained by the positive genetic correlation with flowering time: Early flowering genes tend to be associated with genes reducing the sensitivity to vernalization cues.

CONCLUSIONS

Because it is difficult to rule out the effect of drift on the observed genetic change in phenology, our results do not entirely answer the question of the adaptive potential in selfing populations raised in the Introduction. Yet, several lines of evidence support the role of selection. First, the observed genetic change is in the direction expected for a response to raising temperatures and reduced rainfalls in the Mediterranean region. Second, the selection gradient measured in the greenhouse suggests that early flowering genotypes produce more seeds. The changes in MLG composition through time provide more equivocal results, but are also compatible with the hypothesis that MLGs with early flowering times had a better reproductive success than later flowering genotypes and replaced them, resulting in the observed genetic change in flowering time. Our simulations of the effect of drift are impacted by uncertainty in effective population size estimations, but the highest effective population size compatible with the observed change caused by drift alone remains relatively low (, Figure 4a). Finally, the shift in flowering date observed in the Cape Corsica population was also detected at the regional scale, which suggests that the set of populations studied could be geographic replicates for this response to the selection of flowering times in M. truncatula in Corsica. Ultimately, only a longer survey of this population combined with a pattern test (Sheets & Mitchell, 2001) could provide a definitive answer to the question of adaptation to climate change through a genetic change in flowering time in this predominantly selfing population. Finally, it is worth pointing out that, in contrast with the theoretical predictions presented in the Introduction, this population displays significant genetic variance for a quantitative trait such as flowering time. As suggested before for M. truncatula (Jullien et al., 2019), it is likely that other evolutionary mechanisms, such as migration, contribute to maintain the adaptive potential of populations in this predominantly selfing species.

CONFLICT OF INTEREST

The authors of this article declare that they have no financial conflict of interest with the content of this article.

AUTHOR CONTRIBUTIONS

Laurène Gay: Conceptualization (equal); data curation (lead); formal analysis (lead); funding acquisition (equal); investigation (lead); methodology (lead); project administration (equal); resources (equal); software (equal); supervision (lead); validation (lead); visualization (lead); writing—original draft (lead); writing—review and editing (lead). Julien Dhinaut: Data curation (supporting); formal analysis (supporting); investigation (supporting); methodology (supporting). Margaux Jullien: Data curation (supporting); formal analysis (supporting); investigation (supporting); methodology (supporting). Renaud Vitalis: Formal analysis (supporting); investigation (supporting); methodology (supporting); software (supporting); validation (supporting); writing—original draft (supporting); writing—review and editing (supporting). Miguel Navascues: Formal analysis (supporting); investigation (supporting); methodology (supporting); project administration (supporting); software (supporting); validation (supporting); writing—original draft (supporting); writing—review and editing (supporting). Vincent Ranwez: Methodology (supporting); software (equal); supervision (supporting). Joëlle Ronfort: Conceptualization (equal); data curation (equal); formal analysis (supporting); funding acquisition (equal); investigation (supporting); methodology (supporting); project administration (equal); resources (equal); supervision (equal); validation (supporting); visualization (supporting); writing—original draft (supporting); writing—review and editing (supporting).

OPEN RESEARCH BADGES

This article has been awarded Open Data and Open Materials Badges. All materials and data are publicly accessible via the Open Science Framework at https://data.inrae.fr/dataset.xhtml?persistentId=doi:10.15454/ZY83BE; https://github.com/laugay/GenetHaplo.

Appendix S1

Click here for additional data file.