Open Access

Assortative mating and gene flow generate clinal phenological variation in trees

BMC Evolutionary Biology201212:79

DOI: 10.1186/1471-2148-12-79

Received: 5 December 2011

Accepted: 24 May 2012

Published: 8 June 2012



On-going climate change is shifting the timing of bud burst (TBB) of broad leaf and conifer trees in temperate areas, raising concerns about the abilities of natural populations to respond to these shifts. The level of expected evolutionary change depends on the level and distribution of genetic variation of TBB. While numerous experimental studies have highlighted the role of divergent selection in promoting clinal TBB differentiation, we explored whether the observed patterns of variation could be generated by the joint effects of assortative mating for TBB and gene flow among natural populations. We tested this hypothesis using an in silico approach based on quantitative genetic models.


Our simulations showed that genetic clines can develop even without divergent selection. Assortative mating in association with environmental gradients substantially shifted the mean genetic values of populations. Owing to assortative mating, immigrant alleles were screened for proximal or distant populations depending on the strength of the environmental cline. Furthermore, we confirmed that assortative mating increases the additive genetic variance within populations. However, we observed also a rapid decline of the additive genetic variance caused by restricted gene flow between neighboring populations resulting from preferential matings between phenologically-matching phenotypes.


We provided evidence that the patterns of genetic variation of phenological traits observed in forest trees can be generated solely by the effects of assortative mating and gene flow. We anticipate that predicted temperature increases due to climate change will further enhance genetic differentiation across the landscape. These trends are likely to be reinforced or counteracted by natural selection if phenological traits are correlated to fitness.


Apical bud phenology of temperate trees has been intensively studied in recent years owing to predicted shifts in the timing of bud development as a result of climate changes [1]. Monitoring of leaf unfolding in various species across their distributions has shown that global warming will trigger earlier flushing [24]. These observations have raised concerns about the capacity of tree populations to cope with changes in the timing of bud burst (TBB), which is related to the fitness of trees in two ways: (i) it establishes the length of the growing season and is a major determinant of growth [5], (ii) it determines the timing of flowering, so is related to fecundity [6]. The adaptive response of TBB to global warming is dependent on the level and distribution of genetic variation within a species; the more variation, the larger the predicted genetic shift in TBB. Numerous investigations involving common garden experiments have demonstrated that TBB exhibits large intra- and inter-population differences, as shown by high population differentiation (Q ST ) associated with high heritability values [7]. Additional genetic investigations indicated that juvenile-mature correlation in TBB is high and genotype-environment interactions are low [8]. Finally, genetic dissection by quantitative trait loci (QTLs) mapping has shown that many QTLs contribute to TBB, but these QTLs show stable expression over years and sites [9].

Regardless of species, TBB follows strong geographic clinal patterns of variation, either altitudinal, latitudinal or longitudinal. Phenotypic clines revealed by in situ observations of TBB show congruent patterns across species: bud burst in southern latitudes or lower altitudes occurs earlier than in northern latitudes or higher altitudes [1012], because TBB is triggered by heat sum [13]. Genetic clines can be assessed in common garden experiments where TBB is observed under the same environmental conditions for all populations and are illustrated by the linear relationships between TBB of different populations and geographic variables. Interestingly, genetic clines vary across species and exhibit co-gradient variation or counter-gradient variation with geographic variables and associated phenotypic clines [14, 15]. Co-gradient variation corresponds to clines of both phenotypic variation and genetic variation in a species that co-vary in the same way with the environmental gradient. Counter-gradient variation occurs when phenotypic and genetic clines vary in opposite directions. In the case of oak, genetic and phenotypic clines exhibit co-gradient variation; e.g. populations from southern latitudes flush earlier than populations from northern latitudes, when assessed under the same conditions in common gardens [16, 17]. In the case of beech, genetic clines are opposite to phenotypic clines and exhibit counter-gradient variation: provenances from northern latitudes flush earlier than populations from southern latitudes [18, 19].

Clinal variations, either co- or counter-gradient, have usually been interpreted as consequences of divergent selection among populations by either abiotic or biotic selection pressures. For example, late-flushing trees will not suffer the detrimental effects of late frosts [20] or may avoid damage by defoliating insects [21, 22]. However, few studies have considered the impacts of other evolutionary factors, such as gene flow in combination with the peculiar features of bud burst, in shaping the genetic variation of TBB. Indeed, because trees mate assortatively by flowering time [23, 24], and because TBB is tightly linked to the timing of flowering, assortative mating is likely to shape the variation of TBB. Furthermore, under assortative mating, immigrant pollen will introduce genes likely to generate new allelic combinations for TBB, owing to the existence of environmental clines.

A number of theoretical studies have dissected the effects of assortative mating on the evolution of quantitative traits under polygenic inheritance, beginning with the early investigations by Fisher (1918) [25] and Wright (1921) [26]. All predicted that assortative mating will increase genetic variation as a result of the build up of genetic covariations among loci [25, 2729]. Others demonstrated the amplifying role of assortative mating on natural selection [24, 30], as well as its contribution to allopatric speciation [31, 32]. Finally, more recent studies aimed at predicting the effects of assortative mating on the genetic covariance of different traits [3335]. No prior investigations, however, have considered the effects of assortative mating on a trait in multiple populations interconnected by extensive gene flow in the presence of environmental gradients. We tested whether interactions between gene flow and assortative mating under such circumstances could generate the distribution of genetic variation that is observed in common garden experiments, even in the absence of divergent selection. Our main hypothesis was that assortative mating, by filtering incoming alleles among interbreeding populations, will change the genetic composition and the genetic values of the phenological trait in recipient populations and hence generate population differentiation. We mainly focused on the maintenance of high within- and between-population genetic variation and on the build-up of genetic clines. There exists no available analytical theoretical prediction of genetic variation and differentiation taking into account assortative mating. We therefore used a simulation approach allowing us to monitor in silico the evolution of TBB under contrasting levels of assortative mating and environmental clines.


Components of population subdivision

Our main objective was to track components of genetic variation in phenology-related traits in a subdivided population that would mimic extant ecological settings. We were primarily interested in assessing the within- and between-population genetic variances (V W and V B ) as well as the differentiation among populations as measured by Q ST , which are standard genetic measurements used in quantitative genetics.
Q ST = V B V B + 2 V W
where V W is the within-population genetic variance, and V B is the between-population genetic variance. As suggested by recent QTL studies [9, 36], we assumed that phenological traits were controlled by multiple QTLs with only additive effects. Previous theoretical studies have also shown that the genetic variances V B and V W of multilocus traits can be substantially inflated by allelic covariations among loci [37].
V = i σ i 2 + i j i Co v ij

where σ i 2 is the genic variance of locus i and Co v ij is the covariance between allelic effects at locus i and j. V stands for V B or V W with appropriate σ i 2 and Co v ij expressed either at within- or between-population levels.

These covariations build up as a result of within- or between-gametic disequilibrium generated by different evolutionary forces and are scaled by the parameters θ W and θ B .
θ = i = 1 n j i n Co v ij i = 1 n σ i 2
Le Corre and Kremer (2003) [37] and Kremer and Le Corre (2011) [38] showed how the θ values contributed to the final differentiation of the trait together with the genetic differentiation that also arises at the QTLs controlling the trait ( G S T q
Q ST = ( 1 + θ B ) G S T q ( θ B θ W ) G S T q + 1 + θ W

A major finding of previous theoretical work was that divergent selection generates important between-population disequilibria that becomes a major driver of population differentiation (Q ST ) and has only a minor impact on differentiation at QTLs ( G S T q In the absence of selection and under random mating, θ W and θ B should be 0 and Q ST equal to G S T q We will explore in these simulations how assortative mating will shape the distribution of genetic variability by monitoring the different components of Q ST (e.g. V W , θ W , θ B , V B , and G S T q under different evolutionary scenarios.

Models and simulations

We used the Metapop simulation engine to assess evolutionary changes along successive generations in a subdivided population. Essential steps of the evolutionary processes included in the software - mutation, gene flow, selection, demographic growth - have been described in earlier papers [37, 3941]. We will only address here the changes introduced to account for assortative mating and phenotypic clines of phenological traits.

Phenotypic subdivision of phenological traits

Populations are positioned on a two-dimensional grid (Figure 1) that mimics in a discrete way real situations showing continuous environmental variations. Each population is composed of N individuals. The overall phenotypic value Z ij of individual i from population j is composed of three components: the additive part G ij of the genes contributing to the trait, the environmental component E j and a random local environmental deviation ij .
Z ij = G ij + E j + ij
Figure 1

Spatial settings of populations and environmental effects. Fifty-five populations of 500 individuals each were spread homogeneously on a 5×11 grid along 11 latitudinal positions. E(Y) represents the environmental effect at a given latitude Y and is scaled by k E (see equation 8). No selection was introduced: stabilizing selection was canceled with ω2 = 109 and all populations shared a phenotypic optimum Z opt  = 0.

And the within-population phenotypic value is
Z ij = G ij + ij
G ij is the genetic value resulting from the sum of additive effects of alleles present at n QTLs controlling the trait.
G ij = l = 1 n ( α 1 + α 2 ) l

α values are drawn at loci from the distribution N ( 0 , W l σ A 0 2 / 2 ) where W l is the level of contribution of the lth locus considered and σ A 0 2 the initial variance of allelic effects based on estimated values of additive variance in experimental plantations. More details on the method are available in [38].

E j represents the influence of environmental conditions at the location of population j. E j is of the same magnitude for all individuals of population j located at latitude Y . In our study case, E accounts for the effect of temperature on TBB demonstrated in forest trees [11]; indeed, flushing dates of broadleaves and conifers are tightly dependent on the heat sum [13] and exhibit continuous variation with latitude, resulting in environmental clines of E values. This is the rationale of assigning the same E j value to all trees of population j. The linear variation of E j along latitude, which corresponds to an environmental cline, results in the phenotypic cline as observed in natura (Figure 2). The steepness of the environmental cline is scaled by k E , a standardized measure of the between-environment variance relative to the within-population phenotypic variation. We considered different levels of steepness of the environmental cline by taking different values of k E .
k E = σ E 2 ( σ G 0 2 + σ 2 )
Figure 2

An example of environmental and genetic clines for time of bud burst in oaks (data of Alberto et al., 2011 [12]). The time of bud burst (TBB) was recorded in sessile oak (Quercus petraea) stands located along two valleys on the northern side of the Pyrénées mountains. In situ observations (green dots on the graph) showed that trees located at higher elevations flushed much later then trees located at lower altitudes, as a result of strong correlations between TBB and heat sum [4]. This pattern of variation, the phenotypic cline, is clearly linear. Open-pollinated seeds were collected in each stand and were experimentally raised in a common garden at low altitude, and TBB was monitored (blue points). The TBB was plotted as a function of the altitudes where the seeds were collected. A linear pattern of variation corresponds to a genetic cline. This example illustrates a co-gradient pattern of variation, because the slopes of the phenotypic and genetic clines share the same sign. Counter-gradient variation corresponds to cases where the two clines vary in opposite directions.

σ G 0 2 being the total genetic variance observed within the initial population. Hence k E is constant over the generations through the evolutionary process. Given that E follows a linear relationship with latitude, we can assign environmental values E j according to
E j = k E × ( σ G 0 2 + σ 2 ) σ Y 2 × Y j

Finally, ij is a random local environmental deviation following the distribution N ( 0 , σ )

Sequence of evolutionary processes in Metapop

Metapop implements evolutionary processes over successive generations in a subdivided population. Within each generation, processes are simulated along four steps within a main loop, depicted in Additional file 1: Figure S1. First, fitnesses of reproducing individuals are computed according to stabilizing and divergent selection. The level of stabilizing selection is scaled by the parameter ω2 from Turelli’s relation [42] while the strength of divergent selection is scaled by σ Z opt 2, where Z opt of a given population is the phenotypic value for which trees have the highest fitness in that population. Second, from the populations growth settings and seed migration matrix, the number of individuals of each population contributing to the future generation is computed. Third, mates are chosen based on the constraints due to assortative mating scaled by the correlation between Z i and Z j , the overall phenotypic values of individuals i and j.
ρ = cov ( Z i , Z j ) σ Z i × σ Z j
Following 10, the differences in phenotypic values of two mating parents are drawn from the distribution N ( 0 , σ δ with
σ δ = σ Z i 2 ρ 2 σ Z i 2

Fertilization occurs by drawing male and female gametes conditionally to ρ, fitness of the parents and seeds migration matrix. A proportion of male gametes, based on the pollen migration matrix, is drawn from other populations to account for pollen flow. Finally, mutation is also considered.

Monitoring of gene flow

We now consider how the interaction between gene flow and assortative mating may modify the genetic values in natural populations. Because assortative mating will filter immigrant alleles so that they can mate with trees of recipient populations, we compare the genetic values of immigrant alleles to local alleles to explore whether gene flow will modify the mean genetic value of populations.

In each generation, matings take place between trees of the same population, but a fraction m p of matings involves pollen from other populations. We can subdivide the genetic value of the offspring into two components:
G t + 1 = ( 1 m p ) ( 1 2 G t + 1 2 G t ) + m p ( 1 2 G t + 1 2 G t )
where G t and G t stand respectively for the mean genetic values of the female and male parents, and G t stands for the mean genetic value of the male parents providing immigrant alleles at generation t. ( 1 m p ) ( 1 2 G t + 1 2 G t ) represents the component of the genetic value due to intra-population matings and m p ( 1 2 G t + 1 2 G t ) the component of the genetic value due to inter-population matings involving external incoming alleles. Each generation, G t = G t = G t When assortative mating occurs within populations, mating parents share similar phenotypic values, and because they belong to the same population, they also share the same environmental values. However, because male parents from the outside populations should share the same phenotypic value as the female parent, their genetic values are likely to be different from those of the female parents owing to the environmental gradient. Within a population, the mean phenotypic value of the male parents corresponding to the immigrant alleles is equal to
Z t = G t + E
and the mean phenotypic value of the female parents is equal to
Z t = G t + E
Because the phenotypic values of both parents should be similar owing to assortative mating, the mean genetic value of the male parents is
G t G t + E E
As a result, each generation the genetic value of the population is expected to shift by about Δ=Gt + 1G t , which can be expressed in
Δ 1 2 m p ( E E )

More generally, matings that occur within populations can be subdivided in two different kinds: (1) matings between individuals sharing similar genetic values, which would correspond to positive assortative mating and (2) matings between individuals likely to have different genetic values resulting from gene flow. In the extreme case, these matings may result from negative assortative mating. The shift of the genetic value is therefore driven by the level of effective gene flow m p and the difference in environmental values between the recipient and donors populations. Consequently, we monitored the effective pollen flow during the simulations by tracking its spatial origin.

Simulations settings

We simulated the evolution of 55 populations of 500 individuals each spread homogeneously on a 5 × 11 grid depicted in Figure 1. We did not consider overlapping generations and the number of individuals per population was kept constant over successive generations. A fictive gradient of latitudes was set from latitude −0.5 to latitude + 0.5 in steps of 0.1. Three levels of environmental clines were considered along the latitudinal gradient: k E  = 1, k E  = 2 and k E  = 3.

Recent observations in oak populations suggested that assortative mating for TBB is substantial [6]. Indeed, the flowering time in oak may extend over several weeks within a population, but the receptive period of female flowers lasts only a few days at the individual level. We consequently investigated two strengths of assortative mating, encompassing the suspected range of variation, using ρ = 0.3 and ρ = 0.8 to model moderate and strong assortative mating, respectively. Random mating was considered as well with ρ = 0.

We used Wright’s island migration model to generate gene flow among populations located on the grid system, and considered two levels of gene flow: Nm = 5.1 and Nm = 10.2. These values fit the range of variation of F ST values (2.4% to 4.7%) observed in natural oak populations [7]. Pollen and seed migration rates (m p and m s ) were then inferred from Nm values and introduced in the simulations, assuming further that m p  = 100m s (Table 1). In addition to the island model, we also designed gene flow via the stepping stone model using pollen and seed migration rates corresponding to Nm = 5.1.
Table 1

Initial simulation settings

heritability h2


selfing rate s


nb. of populations d


nb. of ind. per pop. N ind


pollen migration rates m p

0.02, 0.04

seed migration rates m s

0.0002 , 0.0004

nb. of QTL n


mutation rate μ


nb. of latitude levels Y


interval of latitudes Y

0 . 5 , + 0 . 5

steepness of environmental cline E scaled by k E

1, 2, 3

variance of Z opt across latitudinal levels σ Z opt 2

0, 1

intensity of stabilizing selection ω2

109, 5

assortative mating intensity ρ

0; 0.3; 0.8

Assuming that the starting populations were in mutation-migration-drift equilibrium, initial allelic frequencies in different populations were drawn from a Dirichlet distribution [38]. We assumed that phenological traits were controlled by 10 QTLs. Additive values of alleles were chosen at random from a Gaussian distribution whose initial variance was adjusted to fit the heritability values observed in extant progeny plantations, 0.83 from [43]. Mutations at each QTL occurred across generations at a rate of μ = 10−5 per generation. The local environmental deviation was drawn at random from the distribution N ( 0 , 1 )

We considered eight different evolutionary scenarios by combining unique slopes of environmental clines, levels of assortative mating, migration models, and levels of gene flow (Table 2). Because our investigations were focused on the impact of gene flow and assortative mating on the evolution of TBB, we purposely excluded selection in the simulations. We consequently canceled stabilizing selection within all populations by setting all ω2 values to 109, and we set all Z opt values to 0. However, as a control, we added one scenario including selection (ω2 = 5 and σ Z opt 2 = 1, corresponding to strong stabilizing selection and moderate divergent selection. This scenario did not consider assortative mating and was designed to compare the steepness of the genetic clines observed in the eight studied cases with a selective scenario. For each evolutionary scenario based on combinations of these settings (Table 2), we performed 50 independent replicated simulations over 1000 generations.
Table 2

Evolutionary scenarios


ρ = 0

ρ = 0.3

ρ = 0.8

k E  = 1




k E  = 2



×,× s m

k E  = 3



*identical scenarios; because under random mating, phenotypic values of individuals have no impact on our simulation outcomes, variations in the environmental component have no influence when ρ = 0. ×s and × m stand respectively for scenarios simulated under the stepping-stone migration model and with a higher migration rate (Nm = 10.2) under the island migration model.


Within population genetic variance

Assortative mating substantially increased allelic covariances during the first generations (Figure 3). After reaching maximum values, covariances decreased very rapidly and evolved to asymptotic levels. These patterns were more pronounced when assortative mating was strong and were only slightly modified by the magnitude of the environmental cline. Under strong assortative mating, covariances accounted for more than 1.5 of the genic variances relative to the total genetic variance, while under moderate assortative mating, the maximum value was only 0.28. Under steeper clines, the maximum values of θ W were slightly higher, 1.5 vs 1.4, and its change over generations was slightly delayed. Overall θ W values were always larger under assortative mating than under random mating.
Figure 3

Variations in within-population allelic covariation ( θ W ) and genetic variance ( V W ) under different evolutionary scenarios. θ W and V W were monitored under three different strengths of assortative mating and two levels of environmental cline. All simulations were conducted under the island migration model with moderate gene flow (Nm = 5.1). The red line indicates strong assortative mating (ρ = 0.8), the blue line moderate assortative mating (ρ = 0.3), and the black line random mating (ρ = 0). Each line represents the mean of 50 independent replicates for each evolutionary scenario.

The variations in θ W had striking consequences on the genetic variances (equation 2). Indeed, under assortative mating, genetic variances increased rapidly during the early generations, then they very rapidly dropped below even the level of genetic variance reached under random mating. As for covariances, there was a strong effect of the level of assortative mating and only a minor effect of the environmental cline. The decrease in genetic variance due to assortative mating could be dramatic after 400 generations. Furthermore, the final heritability for the trait was divided by a factor 2.5 at generation 500. As expected without selection in large populations, genetic variance was maintained under random mating and extensive gene flow.

Between population genetic variance

Assortative mating had a strong effect on allelic covariances at the between-population level; θ B increased during the early generations and was maintained at higher values through the 1000 generations, in contrast to θ W values. There was a stronger impact when environmental clines were steeper. For example, under strong assortative mating, the maximum value of θ B was 2.7 when k E  = 2vs 2.5 when k E  = 1. The initial phase of increase lasted longer under moderate assortative mating than under strong assortative mating: 500 generations vs 230 generations when k E  =  1 (Figure 4).
Figure 4

Variations in between-population allelic covariation ( θ B ), between-population variation ( V B ), and differentiation of TBB ( Q ST ). These measurements were monitored under three different strengths of assortative mating and two levels of environmental cline. All simulations were conducted under the island migration model with moderate gene flow (Nm = 5.1). The red line indicates strong assortative mating (ρ = 0.8), the blue line moderate assortative mating (ρ = 0.3), and the black line random mating (ρ = 0). Each line represents the mean of 50 independent replicates for each evolutionary scenario.

Between-population variances of allelic frequencies at selected loci increased steadily over generations. They increased more rapidly under strong assortative mating, while no substantial differences were observed between random mating and moderate assortative mating. By generation 1000, differentiation at selected loci had reached 0.16, which could be compared with differentiation under random mating (0.03), which was very close to differentiation at neutral markers (0.024) (data not shown). Overall, between-population genetic variances exhibited strong differences between moderate and strong assortative mating and also between low and strong environmental clines (Figure 4).

Trait differentiation and genetic clines

Because assortative mating had strong consequences on within- and between-population genetic variances, it ultimately contributed to population differentiation of the trait. There were striking differences in the levels of differentiation observed under random and assortative mating. Q ST values steadily increased under assortative mating and reached up to 0.7 when k E  = 2. There was only a slight effect of the steepness of the environmental cline on the level of differentiation: Q ST  = 0.7 when k E  = 2vs 0.62 when k E  = 1.

This effect was due to the trade-off between variations in V B and V W in equation 1. The steepness of the environmental cline increased V W (Figure 3) and had a decreasing effect on Q ST , but at the same time, it also increased V B , increasing Q ST (Figure 4). As a result, Q ST showed similar values at both levels of environmental cline. These results suggested that assortative mating differentiated populations and shifted their mean genetic values. We consequently examined the spatial distribution of mean genetic values across the landscape; indeed, a cline of genetic values built up during the early generations following a south-north gradient (Figure 5). The steepness of the genetic cline was stronger under assortative mating and under steep environmental clines resulting in a co-gradient variation with the environmental cline. The temporal dynamics of the cline could be illustrated by the changes in the genetic value of the population located at the extreme northern latitude (Figure 6). This value reached a peak between generation 200 and 400, depending on the steepness of the environmental cline and the level of assortative mating. No genetic cline developed under random mating.
Figure 5

Variations in mean population genetic values at different latitudes and at different generations. The value for each latitude is the average of the five mean genetic values for the populations concerned. Latitudinal means were computed and reported for two levels of environmental cline and three different strengths of assortative mating. All simulations were conducted under the island migration model with moderate gene flow (Nm = 5.1). The red line indicates strong assortative mating (ρ=0.8), the blue line moderate assortative mating (ρ = 0.3), and the black line random mating (ρ = 0). The dashed line depicts the mean genetic value obtained under divergent selection modeled with ω2 = 5 and σ Z opt 2 = 5 and without assortative mating. Each line represents the mean of 50 independent replicates for each evolutionary scenario.
Figure 6

Evolution of the mean genetic value of a population located at the extreme north of the landscape. The mean genetic value of a population located at latitude + 0.5 (dotted circle in Figure 1) was monitored under two different levels of environmental cline and three different strengths of assortative mating. All simulations were conducted under the island migration model with moderate gene flow (Nm=5.1). The red line indicates strong assortative mating (ρ=0.8), the blue line moderate assortative mating (ρ=0.3), and the black line random mating (ρ = 0). Each line represents the mean of 50 independent replicates for each evolutionary scenario.

We also explored the clinal patterns resulting from a more extreme environmental cline, a higher migration rate, and the stepping-stone migration model (Figure 7). Surprisingly the resulting genetic cline was less pronounced under k E  = 3 than under k E  = 2. When k E  = 3, the environmental variance among populations was 3-fold larger than the within-phenotypic variance (equation 8). Consequently, phenological matches between trees from different populations were limited, thus increasing the filtering of incoming genes to proximal populations (Figures 8 and 9). Similarly, when the pollen dispersal distance was a priori reduced to the most proximal populations, as in the stepping-stone migration model, a very shallow genetic cline built up (Figure 7). In this latter case, when Nm = 5.1, ρ = 0.8, and k E  = 2, only populations at extreme latitudes became genetically differentiated. Despite this very shallow cline, Q ST approached 0.45 at generation 1000 under the stepping-stone migration model; under the same simulations parameters, Q ST values reached 0.7 under the island migration model. Finally, when pollen migration rates increased (Nm = 10.2vs Nm = 5.1), no significant change was observed in the slopes of the clines. However, additional investigations indicated that lower migration rates decreased the slopes of the genetic clines and induced higher Q ST values, owing to an important drift effect [37] (Additional file 2: Figure S2 and Additional file 3: Figure S3). Overall large stochastic variations were associated with the genetic parameters that were monitored during the evolutionary scenarios (data not shown). We illustrate these variations only for Q ST and V W (Figure 10). The trend among generations, i.e., the form of the curve, was the same among the replicates.
Figure 7

Variations in mean population genetic values at different latitudes under multiple scenarios. The value for each latitude is the average of the five mean genetic values for the populations concerned at generation 300. All scenarios (except the selection scenario, dashed line) were conducted under strong assortative mating (ρ = 0.8). Red line: steep environmental cline (k E  = 2), island migration model, moderate gene flow (Nm = 5.1). Purple line: very steep environmental cline (k E  = 3), island migration model, moderate gene flow (Nm = 5.1). Brown line: steep environmental cline (k E  = 2), island migration model, extensive gene flow (Nm=10.2). Red line with open circles: steep environmental cline (k E  = 2), stepping stone migration model, high gene flow (Nm=5.1). Dashed line: random mating, divergent selection ( σ Z opt 2 = 1, strong stabilizing selection (ω2 = 5), without assortative mating. Each line represents the mean of 50 independent replicates for each evolutionary scenario.
Figure 8

Amount of immigrant alleles received by a northern population. Absolute number of immigrant alleles into a population located at the extreme northern latitude (+0.5 dotted circle in Figure 1). Numbers on the y-axis are cumulative counts of alleles from generation 16 to 20. Counts of alleles were monitored at three strengths of assortative mating and three levels of the environmental cline. The red line indicates strong assortative mating (ρ = 0.8), the blue line moderate assortative mating (ρ = 0.3), the black line random mating (ρ = 0), and the purple line strong assortative mating under an extreme environmental cline (ρ = 0.8, k E  = 3). Lines are mean values of 50 replicates for each evolutionary scenario.
Figure 9

Amount of southern immigrant alleles received by a northern population over generations. Absolute number of immigrant alleles into a population located at the extreme northern latitude (+0.5 dotted circle in Figure 1) and coming from southern latitudes (-0.5 to -0.1). Only gene flow between populations is represented here. Numbers on the y-axis are counts of alleles at a given generation (x-axis). Counts of alleles were monitored at three strengths of assortative mating and three levels of environmental cline. All simulations were conducted under the island migration model with moderate gene flow (Nm = 5.1). The red lines indicate strong assortative mating (ρ = 0.8), the blue line indicates moderate assortative mating (ρ = 0.3), the black line random mating (ρ  = 0), and the purple line strong assortative mating under an extreme environmental cline (ρ = 0.8, k E  = 3). Lines are mean values of 50 replicates for each evolutionary scenario.
Figure 10

Stochastic variations in Q ST and V W among different simulations within a given scenario. Upper and lower bounds of the 50 simulations conducted per scenario. ρ was set to 0.8 in all cases. k E is the scaling factor of the environmental cline. Plain lines indicate mean values of the 50 simulations for each scenario and dotted lines represent the two simulations that gave the extreme results.

Pollen filtering by assortative mating

We monitored the incoming pollen composition in a population located at the extreme northern latitude. By doing so, we expected to predict the shift in genetic values that contributed to the development of the genetic cline under the island migration model (equation 16). Figure 8 clearly shows that assortative mating filtered incoming alleles by geographic origin. Very rapidly, there was a preferential screening of incoming alleles from neighboring populations in the case of assortative mating, and the trend was more pronounced when the environmental cline grew steeper. The discrepancy between distant and proximal alleles was more pronounced with strong assortative mating. Furthermore, the level of filtering changed over generations. More alleles arrived from distant populations during the first 40 generations, especially when strong assortative mating was occurring (Figure 9). These distant alleles would shift the genetic values of populations as predicted by Δ.


Our simulations demonstrated that genetic clines could be established in the absence of divergent selection. We showed that the combination of assortative mating and pre-existing environmental clines resulted in population genetic differentiation along the environmental cline. We also confirmed that assortative mating increased the within-population genetic variances in the early stages of the evolutionary scenarios. However, assortative mating was also responsible for the severe decline in genetic variation in later generations.

These patterns resulted in a positive covariance between genetic and environmental population values and corresponded to what has been called co-gradient variation [14, 15]. We discuss here how such covariations may build up under assortative mating in the case of phenological traits in trees. Given the pre-existence of environmental clines, genetic clines are generated by the combined effects of assortative mating and gene flow. In particular, we examine how the interplay between assortative mating and gene flow will actually produce the genetic cline we observed. According to equation 15, the larger the physical distance between the mates associated by gene flow, the more different their genetic values. As a consequence, a larger shift in the mean genetic value should be expected at extreme latitudes in our grid settings (Figure 1). In what follows, we illustrate this trend by providing values for the shift Δ obtained at the extreme northern latitude under the strongest assortative mating intensity and across the steepest environmental cline.

We can subdivide the evolutionary process into three main phases, illustrated in Figures 4, 5, 6, 7 and 8.
  1. (1)

    In the very early generations (05), the mean genetic value is 0 for all populations, there is no within-population allelic covariance, and alleles are randomly spread over the landscape. During this period, assortative mating will generate phenotypes with extreme genetic values in each population. Hence the genetic variance within populations increases as predicted by previous analytical models [24, 34] and numerical simulations [28, 32, 44]. Gene flow during the early generations preferentially imports alleles from neighboring populations (Figure 8), owing to the fact that populations at this stage are genetically undifferentiated over the whole grid and parents exhibiting similar phenotypes are more likely to be in neighboring populations. As a result, the shift Δremains limited: 0.0798 at the allelic level for northern populations.

  2. (2)

    From generation 5 to about 30, because the increase in within-population genetic variance has now produced phenotypes with more extreme values, gene flow tends now to import alleles from more distant populations (Figure 9). The fraction of imported alleles enriches the population gene pool and further facilitates an increase in genetic covariances θ W . The genetic variance between populations continues to increase steadily. During the second phase, the Δvalue tends to be larger (0.14) as a result of more divergent alleles imported by distant gene flow. A similar effect that symmetrically decreases the Δvalue of incoming gene flow within southern populations is expected to take place at the same time. As a result, the mean genetic values of the population shift strongly, leading to the progressive formation of the genetic cline.

  3. (3)

    After generation 30, most of the alleles have been spatially redistributed by gene flow constrained by assortative mating at the landscape level. Allelic covariations within populations have been exhausted and the genetic variance has now reached its maximum. Assortative mating within populations tends now to become a selective factor favoring phenotypes following the shift of the mean genetic values. Furthermore, gene flow again becomes strongly restricted to neighboring populations that share fewer divergent alleles than distant populations. Restricted gene flow therefore reinforces the decrease in the genetic variance. Overall, phase 3 is characterized by a continuous decrease in genetic variance and the reaching of an asymptotic mean genetic value in populations; the genetic cline is establishing. We further advocate that restricted gene flow, together with within-population assortative mating, now constrains effective population sizes, accelerating the decrease in genetic variance due to drift. A similar decrease was observed by Devaux and Lande [32] in a single population, despite a high mutation rate. Jorjani et al. also noticed a decreasing effect of negative assortative mating on the evolution of the genetic variance within a single population [44].


These three phases were observed for all of the simulation settings we used. The lengths of the two first phases extended over longer periods, populations differentiated more rapidly, and genetic clines were shaped faster under strong assortative mating. By dissecting the evolutionary process, we showed that the screening of immigrant alleles due to assortative mating triggers shifts in the genetic values of populations (Figures 6 and 9). Indeed, when assortative mating allows for long-distance filtered pollen flow, the shifts in the genetic values of recipient populations are strongly enhanced. Because moderate assortative mating generates less extreme genotypes over generations, distant gene flow is promoted less and the mean expected shift in the mean genetic values of populations remains limited. Consequently, under moderate assortative mating, the final steepness of genetic clines is less dependent on the steepness of environmental clines (Figures 5 and 6).

Increasing the slope of the environmental cline generated more genetic variance and higher genetic differentiation as well. According to equation 15, each generation steeper environmental clines increase the expected divergence between mates from distinct populations. However, the divergence is constrained by the necessary overlap of parental flowering times. If long distance pollen flow is restricted by large phenological differences among populations, then assortative mating will favor matings between proximal populations, and the shift in genetic values will be limited. In our simulations, the latter case occurred when we explored very large k E values (k E  = 3).

A similar outcome was observed under the stepping-stone migration model. In this case, populations do not differentiate except at the northern and southern edges of the landscape (Figure 7). This result is only partly explained by the absence of distant gene flow. Indeed, according to the expression of Δand considering the features of the stepping-stone migration model, limited Δ values are expected owing to pollen flow from adjacent latitudes. However, incoming alleles from neighboring northern populations balance with incoming alleles from neighboring southern populations. As a consequence, the shift induced within populations by southern gene flow is systematically canceled by the one caused by northern flow, resulting in a null contribution to the Δ values. Finally, because under the stepping-stone migration model, incoming gene flow is latitudinally unbalanced at the northern and southern margins of the grid, the genetic values of populations can be shifted by assortative mating at these latitudes. These results suggest that the spatial configuration of the populations in combination with the migration model may also contribute to the building of the genetic cline. Any combination that increases an asymmetry in gene flow between northern and southern populations will enhance the genetic cline, while symmetry will tend to even out the effects of northern and southern gene flow.

To summarize, the construction of a genetic cline as a result of the combined effects of gene flow and assortative mating can only be met under certain circumstances when there is a balance between the intra-population and between-population phenotypic variance (k E varying between 1 and 3), when long distance pollen flow is possible, and when the patterns of incoming pollen flow at population level are unbalanced regarding the environmental cline. Interestingly these criteria are met under realistic situations. Taking oaks as an example, flushing dates may vary over 5 weeks from southwestern to central France [4], while the same range of variation may be observed between early and late flushing trees in a given forest stand. Viable pollen has also been shown to be dispersed over such distances [45].


Our simulations showed that interaction between assortative mating and gene flow across environmental clines may shape the genetic variability of phenologically-related traits and induce cogradient variation without any divergent selection. We also demonstrated that the extent of genetic variability resulting from assortative mating was related to the patterns of incoming pollen flow at the population level. Because phenotypic clines have been very widely reported in forest trees [2, 11, 18], we suspect that assortative mating and gene flow could actually be responsible for the co-gradient variation observed in some species in common garden experiments [12, 17]. However, most tree species actually exhibit counter-gradient variation [46, 47], suggesting that other evolutionary forces, such as divergent selection, actually counteract the combined effects of assortative mating and gene flow. In a subsequent paper, we will explore how selection interacts with assortative mating and gene flow to generate counter-gradient variation. Finally, our simulations also indicated that very large levels of genetic variation should also be expected within populations, generated by genetic covariances in allelic effects due to assortative mating as predicted by other theories or simulations [24, 25, 32]. Experimental data from progeny tests of forest trees indeed show that heritability values of phenologically-related traits can exceed 0.5, much larger than other phenotypic traits generally assessed in experimental plantations [43]. Furthermore, our simulations predict that the steep increase in genetic variation will be temporary and will be followed by a rapid decrease. Once all covariation has been exhausted, assortative mating will act as a selective force by constraining the synchronicity of male and female flowering periods. Given the large genetic variation still existing in extant forest stands, we suspect that the time of decrease has not yet been reached in natural populations, owing to the long generation times of trees. Finally, our simulations should be prolonged under more realistic ecological settings, including different patterns of gene flow and selection on multiple traits. Both authors read and approved the final manuscript.

Author’s contributions

JPS and AK designed the study and identified the different evolutionary scenario to be tested. JPS adapted the Metapop simulation engine by introducing assortative mating and environmental clines. JPS and AK wrote the paper. Both authors read and approved the final manuscript.



This research was supported by the EU project MOTIVE (project number: 226544). We acknowledge the help of Valérie Le Corre and Frédéric Raspail during the adaptation of Metapop to the study of phenological traits. We thank two anonymous reviewers for constructive comments.

Authors’ Affiliations

Univ. Bordeaux, BIOGECO, UMR 1202


  1. Bertin RI: Plant phenology and distribution in relation to recent climate change. J Torrey Bot Soc. 2008, 135: 126-146. 10.3159/07-RP-035R.1.View Article
  2. Menzel A, Sparks TH, Estrella N, Koch E, Aasa A, Ahas R, Alm-Kubler K, Bissolli P: European phenological response to climate change matches the warming pattern. Global Change Biol. 2006, 12: 1969-1976. 10.1111/j.1365-2486.2006.01193.x.View Article
  3. Nordli O, Wielgolaski F, Bakken A, Hjeltnes S, Mage F, Sivle A, Skre O: Regional trends for bud burst and flowering of woody plants in Norway as related to climate change. Int J Biometeorology. 2008, 52: 625-639. 10.1007/s00484-008-0156-5.View Article
  4. Vitasse Y, François C, Delpierre N, Dufrêne E, Kremer A, Chuinee I, Delzon S: Assessing the effects of climate change on the phenology of European temperate trees. Agric Forest Meteorology. 2011, 151 (7): 969-980. 10.1016/j.agrformet.2011.03.003.View Article
  5. Bennie J, Kubin E, Wiltshire A, Huntley B, Baxter R: Predicting spatial and temporal patterns of bud-burst and spring frost risk in north-west Europe: the implications of local adaptation to climate. Global Change Biol. 2010, 16: 1503-1514. 10.1111/j.1365-2486.2009.02095.x.View Article
  6. Franjic J, Sever K, Bogdan S, Skvorc Z, Krstonosic D, Aleskovic I: Phenological asynchronization as a restrictive factor of efficient pollination in clonal seed orchards of Pedunculate Oak (Quercus robur L.). Croat J For Eng. 2011, 31: 141-156.
  7. Kremer A, Le Corre V, Petit R, Ducousso A: Historical and contemporary dynamics of adaptive differentiation in European oaks. Molecular Approaches in Natural Resource Conservation and Management. 2010, Cambridge: Cambridge University Press
  8. Ekberg I, Eriksson G, Namkoong G, Nilsson C, Norell L: Genetic correlations for growth rhythm and growth capacity at ages 3-8 years in provenance hybrids of picea-abies. Scand J Forest Res. 1994, 9: 25-33. 10.1080/02827589409382809.View Article
  9. Derory J, Scotti-Saintagne C, Bertocchi E, Dantec LL, Graignic N, Jauffres A, Casasoli M, Chancerel E, Bodénès C, Alberto F, Kremer A: Contrasting relationships between the diversity of candidate genes and variation of bud burst in natural and segregating populations of European oaks. Heredity. 2010, 104: 438-448. 10.1038/hdy.2009.134.PubMedView Article
  10. Worrall J: Temperature - bud-burst relationships in amabilis and subalpine for provenance tests replicated at different elevations. Silvae Genetica. 1983, 32: 203-209.
  11. Vitasse Y, Delzon S, Dufrene E, Pontailler JY, Louvet JM, Kremer A, Michalet R: Leaf phenology sensitivity to temperature in European trees: Do within-species populations exhibit similar responses?. Agric Forest Meteorology. 2009, 149: 735-744. 10.1016/j.agrformet.2008.10.019.View Article
  12. Alberto F, Bouffier L, Louvet JM, Lamy JB, Delzon S, Kremer A: Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient. J Evol Biol. 2011, 24 (7): 1442-1454. 10.1111/j.1420-9101.2011.02277.x.PubMedView Article
  13. Chuine I, Cour P: Climatic determinants of budburst seasonality in four temperate-zone tree species. New Phytol. 1999, 143: 339-349. 10.1046/j.1469-8137.1999.00445.x.View Article
  14. Conover DO, Schultz ET: Phenotypic similarity and the evolutionary significance of countergradient variation. Tree. 1995, 10: 248-252.PubMed
  15. Conover DO, Duffy TA, Hice LA: The covariance between genetic and environmental influences across ecological gradients: reassessing the evolutionary significance of countergradient and cogradient variation. Ann N Y Acad Sci. 2009, 1168: 100-129. 10.1111/j.1749-6632.2009.04575.x.PubMedView Article
  16. Jensen JS, Hansen J: Geographical variation in phenology of Quercus petraea (Matt.) Liebl and Quercus robur L. oak grown in a greenhouse. Scand J Forest Res. 2008, 23: 179-188. 10.1080/02827580801995331.View Article
  17. Vitasse Y, Delzon S, Bresson C, Michalet R, Kremer A: Altitudinal differentiation in growth and phenology among populations of temperate tree species in a common garden. Can J Forest Res. 2009, 39: 1259-1269. 10.1139/X09-054.View Article
  18. Von Wuehlisch G, Krusche D, Muhs H: Variation in temperature sum requirement for flushing of beech provenances. Sylvae Genetica. 1995, 44: 343-346.
  19. Chmura D, Rozkowski R: Variability of beech provenances in spring and autumn phenology. Silvae Genetica. 2002, 51: 123-127.
  20. Howe G, Aitken S, Neale D, Jermstad K, Wheeler N, Chen T: From genotype to phenotype: unravelling the complexities of cold adaptation in forest trees. Can J Bot. 2003, 81: 1247-1266. 10.1139/b03-141.View Article
  21. Asch MV, Tienderen P, Holleman L, Visser M: Predicting adaptation of phenology in response to climate change, an insect herbivore example. Global Change Biol. 2007, 13: 1596-1604. 10.1111/j.1365-2486.2007.01400.x.View Article
  22. Ghelardini L, Santini A: Avoidance by early flushing: a new perspective on Dutch elm disease research. iForest. 2009, 2: 143-153. 10.3832/ifor0508-002.View Article
  23. Kirkpatrick M: Reinforcement and divergence under assortative mating. Proc Biol Sci. 2000, 267 (1453): 1649-1655. 10.1098/rspb.2000.1191.PubMedPubMed CentralView Article
  24. Fox GA: Assortative mating and plant phenology: evolutionnary and pratical consequences. Evolutionnary Ecol Res. 2003, 5: 1-18.
  25. Fisher RA: The correlation between relatives on the supposition of Mendelian inheritance. Trans Roz Soc Edinb. 1918, 52: 339-433.
  26. Wright S: Systems of mating. III. Assortative mating based on somatic resemblance. Genetics. 1921, 6: 144-161.PubMedPubMed Central
  27. Lande R: The influence of the mating system on the maintenance of genetic variability in polygenic characters. Genetics. 1977, 86: 485-498.PubMedPubMed Central
  28. Rosvall O, Mullin TJ: Positive assortative mating with selection restrictions on group coancestry enhances gain while conserving genetic diversity in long-term forest tree breeding. Theor Appl Genet. 2003, 107 (4): 629-642. 10.1007/s00122-003-1318-9.PubMedView Article
  29. Weis AE, Kossler TM: Genetic variation in flowering time induces phenological assortative mating: quantitative genetic methods applied to Brassica rapa. Am J Bot. 2004, 91: 825-836. 10.3732/ajb.91.6.825.PubMedView Article
  30. Jorjani H, Engström G, Strandberg E, Liljedahl LE: Genetic studies of assortative mating - a simulation study. III. Assortative mating in selected populations. Acta Agric Scand. 1997, 47: 129-137.
  31. Caisse M, Antonovics J: Evolution in closely adjacent plant populations. IX. Evolution of reproductive isolation in clinal populations. Heredity. 1978, 40: 371-384. 10.1038/hdy.1978.44.View Article
  32. Devaux C, Lande R: Incipient allochronic speciation due to non-selective assortative mating by flowering time, mutation and genetic drift. Proc Biol Sci. 2008, 275 (1652): 2723-2732. 10.1098/rspb.2008.0882.PubMedPubMed CentralView Article
  33. Gianola D: Assortative mating and the genetic correlation. Theor Appl Genet. 1982, 62: 225-231.PubMed
  34. Hayashi T: Genetic variance under assortative mating in the infinitesimal model. Genes Genet Syst. 1998, 73: 397-405. 10.1266/ggs.73.397.View Article
  35. Weis AE: Direct and indirect assortative mating: a multivariate approach to plant flowering schedules. J Evol Biol. 2005, 18 (3): 536-546. 10.1111/j.1420-9101.2005.00891.x.PubMedView Article
  36. Saintagne C, Bodenes C, Barreneche T, Pot D, Plomion C, Kremer A: Distribution of genomic regions differentiating oak species assessed by QTL detection. Heredity. 2004, 92: 20-30. 10.1038/sj.hdy.6800358.PubMedView Article
  37. Le Corre, Kremer A: Genetic Variability at Neutral Markers, Quantitative Trait Loci and Trait in a Subdivided Population Under Selection. Genetics. 2003, 164: 1205-1219.
  38. Kremer A, Le Corre V: Decoupling of differentiation between traits and their underlying genes in response to divergent selection. Heredity. 2011, 10.1038/hdy.2011.81.
  39. Le Corre V, Machon N, Petit RJ, Kremer A: Colonization with ong-distance seed dispersal and genetic structure of maternally inherited genes in forest trees: a simulation study. Genet Res. 1997, 69: 117-125. 10.1017/S0016672397002668.View Article
  40. Le Corre V, Kremer A: The genetic differentiation at quantitative trait loci under local adaptation. Mol Ecol. 10.1111/j.1365-294X.2012.05479.x.
  41. Machon N, Bardin P, Mazer S, Moret J, Godelle B, Austerlitz F: Relationship between genetic structure and seed and pollen dispersal in the endangered orchid Spiranthes spiralis. New Phytologist. 2003, 157: 677-687. 10.1046/j.1469-8137.2003.00694.x.View Article
  42. Turelli M: Heritable genetic variation via mutation-selection balance: Lerch’s zeta meets the abdominal bristle. Theor Popul Biol. 1984, 25: 138-193. 10.1016/0040-5809(84)90017-0.PubMedView Article
  43. Cornelius J: Heritabilities and additive genetic coefficients of variation in forest trees. Can J Forest Res. 1994, 24: 372-379. 10.1139/x94-050.View Article
  44. Jorjani H, Engström G, Strandberg E, Liljedahl LE: Genetic studies of assortative mating - a simulation study. Assortative mating in unselected populations. Acta Agric Scand. 1997, 47: 74-81.
  45. Schueler S, Schlünzen K, Scholz F: Viability and sensitivity of oak pollen and its implications for pollen-mediated gene flow. Trends Ecol Evol. 2005, 19: 154-161.
  46. Wright JW: Introduction to Forest Genetics. 1976, New-York: Academic Press, INC.
  47. Morgenstern EK: Geographic variation in Forest Trees: Genetic Basis and Application of Knowledge in Silviculture. 1996, Vancouver, Canada: UBC press


© Soularue and Kremer; licensee BioMed Central Ltd. 2012

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://​creativecommons.​org/​licenses/​by/​2.​0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.