Jean-Paul Soularue¹

Cyril Firmat²

Thomas Caignard¹

Armel Thöni³

Léo Arnoux¹

Sylvain Delzon¹

Ophélie Ronce⁴

Antoine Kremer¹

Antagonistic effects of assortative mating on the evolution of phenotypic plasticity along environmental gradients

Keywords: phenotypic plasticity, assortative mating, sexual selection, phenology, individual-based model, stochastic simulations, random-regression mixed effect model, quantitative genetics

Previous theory has shown that assortative mating for plastic traits can maintain genetic divergence across environmental gradients despite high gene flow. Yet, these models did not examine how assortative mating affects the evolution of plasticity. We here describe patterns of genetic variation across elevation for plasticity in a trait under assortative mating, using multiple-years observations of bud burst date in a common garden of sessile oaks. Despite high gene flow, we found significant spatial genetic divergence for the intercept, but not for the slope, of reaction norms to temperature.

We then used individual-based simulations, where both the slope and intercept of the reaction norm evolve, to examine how assortative mating affects the evolution of plasticity, varying the intensity and distance of gene flow. Our model predicts the evolution of, either suboptimal plasticity (reaction norms with a slope shallower than optimal), or hyperplasticity (slopes steeper than optimal) in the presence of assortative mating, when optimal plasticity would evolve under random mating. A co-gradient pattern of genetic divergence for the intercept of the reaction norm (where plastic and genetic effects are in the same direction) furthermore always evolves in simulations with assortative mating, consistently with our observations in the studied oak populations.

Introduction

Adaptation to variable environments requires phenotypic change, which can be achieved through genetic divergence or phenotypic plasticity (Hoffman and Sgro 2011). Phenotypic plasticity is the ability of an organism to express distinct phenotypes in different environmental conditions [START_REF] Pigliucci | Phenotypic plasticity: beyond nature and nurture[END_REF]. Phenotypic plasticity could play a key role in facilitating persistence under environmental change [START_REF] Hendry | Human influences on rates of phenotypic change in wild animal populations[END_REF][START_REF] Chevin | Phenotypic plasticity and evolutionary demographic responses to climate change: taking theory out to the field[END_REF][START_REF] Tufto | Genetic evolution, plasticity, and bet-hedging as adaptive responses to temporally autocorrelated fluctuating selection: A quantitative genetic model[END_REF][START_REF] Scheiner | The genetics of phenotypic plasticity. XV. Genetic assimilation, the Baldwin effect, and evolutionary rescue[END_REF]. This role will however depend on the adaptive nature of phenotypic plasticity and the degree of mismatch between phenotypes expressed in different environments and their optimal values [START_REF] Duputié | Phenological plasticity will not help all species adapt to climate change[END_REF][START_REF] Phillimore | Passerines may be sufficiently plastic to track temperature-mediated shifts in optimum lay date[END_REF][START_REF] Tansey | Estimating the ability of plants to plastically track temperature-mediated shifts in the spring phenological optimum[END_REF][START_REF] Gauzere | Where is the optimum? Predicting the variation of selection along climatic gradients and the adaptive value of plasticity. A case study on tree phenology[END_REF]. Understanding how phenotypic plasticity has evolved in variable and heterogeneous environment is thus critical to forecasting the maximal rate of environmental change that populations could cope with in the future [START_REF] Chevin | Adaptation, plasticity, and extinction in a changing environment: towards a predictive theory[END_REF][START_REF] Hendry | Key questions on the role of phenotypic plasticity in eco-evolutionary dynamics[END_REF][START_REF] Valladares | The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change[END_REF]. Here, we examine the idea that assortative mating could affect the evolution of phenotypic plasticity along an environmental gradient. To do so, we use a general simulation model and compare its predictions to observed patterns of quantitative genetic variation for a plastic trait involved in assortative mating.

Many theoretical studies have showed that costs of plasticity [START_REF] Van Tienderen | Generalists, specialists, and the evolution of phenotypic plasticity in sympatric populations of distinct species[END_REF][START_REF] Dewitt | Costs and limits of phenotypic plasticity[END_REF][START_REF] Sultan | meta-population structure favors plasticity over local adaptation[END_REF], unpredictable selection environments [START_REF] Moran | The Evolutionary Maintenance of Alternative Phenotypes[END_REF][START_REF] Scheiner | The genetics of phenotypic plasticity. VII. Evolution in a spatially-structured environment[END_REF][START_REF] De | Unpredictable selection in a structured population leads to genetic differentiation in evolved reaction norms[END_REF][START_REF] Tufto | The Evolution of Plasticity and Nonplastic Spatial and Temporal Adaptations in the Presence of Imperfect Environmental Cues[END_REF][START_REF] Lande | Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilation[END_REF][START_REF] Scheiner | The genetics of phenotypic plasticity. XII. Temporal and spatial heterogeneity[END_REF][START_REF] Botero | Evolutionary tipping points in the capacity to adapt to environmental change[END_REF], or limited genetic variation for the shape of the reaction norm [START_REF] Via | Genotype-environment interaction and the evolution of phenotypic plasticity[END_REF][START_REF] Scheiner | The genetics of phenotypic plasticity. X. Variation versus uncertainty[END_REF] can impede the evolution of optimal phenotypic plasticity and favor instead the local adaptation of specialist genotypes. Partially adaptive plasticity is then combined with genetic divergence through space [START_REF] Schmid | The role of phenotypic plasticity on population differentiation[END_REF]). Yet, these theoretical models also predict that very high gene flow along such gradients will erode genetic divergence [START_REF] Tufto | The Evolution of Plasticity and Nonplastic Spatial and Temporal Adaptations in the Presence of Imperfect Environmental Cues[END_REF][START_REF] Scheiner | The genetics of phenotypic plasticity. XII. Temporal and spatial heterogeneity[END_REF][START_REF] Schmid | Species' range dynamics affect the evolution of spatial variation in plasticity under environmental change[END_REF]) and lead to the evolution of a single optimal reaction norm in highly predictable environments. The common observation of genetic clines for plastic traits despite high gene flow, as observed for instance for many climate adaptations in forest trees [START_REF] Savolainen | Gene flow and local adaptation in trees[END_REF], may therefore seem paradoxical, especially given the scarce empirical evidence for plasticity costs [START_REF] Van Buskirk | The fitness costs of developmental canalization and plasticity[END_REF][START_REF] Auld | Re-evaluating the costs and limits of adaptive phenotypic plasticity[END_REF][START_REF] Murren | Constraints on the evolution of phenotypic plasticity: limits and costs of phenotype and plasticity[END_REF].

A few theoretical studies have however suggested that the observation of strong genetic divergence despite high gene flow may not be so paradoxical when considering plastic traits under strong assortative mating [START_REF] Stam | The evolution of reproductive isolation in closely adjacent plant populations through differential flowering time[END_REF][START_REF] Soularue | Assortative mating and gene flow generate clinal phenological variation in trees[END_REF], 2014). Assortative mating describes a mating pattern in which there is more similarity between the phenotypes of the mating individuals than would be expected by chance in a population [START_REF] Ennos | Pollen success, functional gender and assortative mating in an experimental plant population[END_REF]. For instance, plants can mate only with individuals flowering at the same time (Weiss et al. 2014). Assortative mating is frequently reported both in animals [START_REF] Jiang | Assortative mating in animals[END_REF] and plants (Weiss et al. 2014), and concerns many plastic traits involved in climatic adaptation, such as phenology. With assortative mating, phenotypic plasticity, affecting trait variation through space, results in non-random patterns of gene flow: for instance, among distant male plants, only those with particularly early flowering can pollinate female plants in environments where all plants flower early because of phenotypic plasticity. This filters out later flowering genotypes from the pollen pool, causing, not only reduced gene flow between distant individuals, but also biased immigration of early flowering genotypes in environments where plasticity causes early flowering. Even in the absence of any divergent selection on the trait value, this filtering results in a spatial reorganization of alleles according to the environmental gradient, and the evolution of co-gradient genetic divergence, where genetic effects align with environmental effects (Stam 1983, Soularue andKremer 2012). When the trait is under selection, this spatial sorting of genotypes due to assortative mating may also enhance or reduce the degree of local adaptation depending on the adaptive or maladaptive nature of phenotypic plasticity [START_REF] Soularue | Evolutionary responses of tree phenology to the combined effects of assortative mating, gene flow and divergent selection[END_REF]. Interestingly, with plasticity and assortative mating, genetic divergence evolves despite very high gene flow. These previous theoretical studies however considered that only the intercept, but not the slope, of the reaction norm was genetically variable and could evolve along the gradient. They therefore did not explore how assortative mating affects the evolution of plasticity along environmental gradients. [START_REF] Nonaka | Assortative mating can limit the evolution of phenotypic plasticity[END_REF] predicted that, in heterogeneous environments with high mixing, intermediate levels of assortative mating could limit the evolution of phenotypic plasticity, even in the absence of plasticity cost. Assortative mating has multifarious and antagonistic consequences on the evolution of quantitative traits, which could in particular affect the evolution of plasticity. On the one hand, assortative mating can increase the genetic variance of polygenic quantitative traits, because it generates positive linkage disequilibrium and heterozygotes deficit [START_REF] Fisher | The correlation between relatives on the supposition of Mendelian inheritance[END_REF][START_REF] Wright | Assortative mating based on somatic resemblance[END_REF] but see [START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF], which in turn can amplify the response to directional selection (see predictions by [START_REF] Fox | Assortative mating and plant phenology: evolutionary and practical consequences[END_REF], Weiss et al. 2005, Godineau et al. 2021), and promote adaptive divergence in heterogeneous landscapes despite high levels of dispersal [START_REF] Doebeli | Speciation along environmental gradients[END_REF][START_REF] Devaux | Incipient allochronic speciation due to non-selective assortative mating by flowering time, mutation and genetic drift[END_REF][START_REF] Sachdeva | Divergence and evolution of assortative mating in a polygenic trait model of speciation with gene flow[END_REF]. On the other hand, assortative mating can also constrain the maintenance of adaptive genetic diversity, because it reduces effective population size [START_REF] Devaux | Incipient allochronic speciation due to non-selective assortative mating by flowering time, mutation and genetic drift[END_REF], and increase stabilizing selection on the trait [START_REF] Kirkpatrick | Sexual selection can constrain sympatric speciation[END_REF][START_REF] Sachdeva | Divergence and evolution of assortative mating in a polygenic trait model of speciation with gene flow[END_REF], Godineau et al. 2021). When assortative mating makes it more likely for individuals with an intermediate phenotype to mate, assortative mating indeed generates a form of stabilizing sexual selection, which can oppose disruptive natural selection [START_REF] Kirkpatrick | Sexual selection can constrain sympatric speciation[END_REF][START_REF] Nonaka | Assortative mating can limit the evolution of phenotypic plasticity[END_REF]. The latter effect explains why [START_REF] Nonaka | Assortative mating can limit the evolution of phenotypic plasticity[END_REF] predicted that intermediate assortative mating could impede the evolution of phenotypic plasticity: stabilizing sexual selection disfavors genotypes expressing extreme phenotypes, even if those are favored by natural selection in different environments. Beyond some critical intensity of assortative mating, intermediate phenotypes loose their mating advantage and assortative mating does not impede the evolution of phenotypic plasticity.

The model by [START_REF] Nonaka | Assortative mating can limit the evolution of phenotypic plasticity[END_REF] considered only two habitats with either complete mixing or no gene flow. They therefore did not predict how reaction norms may differentiate along spatial environmental gradient, nor did they examine the effects of varying the intensity and distance of gene flow on such differentiation. In the present study, we set to fill these gaps by combining experimental data analysis and individual-based modeling to investigate the evolution and spatial genetic differentiation of a plastic quantitative trait with assortative mating, describing genetic variation for both the slope and the intercept of the reaction norm along some environmental gradient.

We first document genetic variation within and between populations for the plasticity of a phenology related trait, the timing of bud burst, in natural sessile oak populations along elevation gradients in the French Pyrenees. The timing of bud burst varies plastically with spring temperature [START_REF] Vitasse | Assessing the effects of climate change on the phenology of European temperate trees[END_REF], is a critical determinant of fitness and tree distributions [START_REF] Chuine | Phenology is a major determinant of tree species range[END_REF], and is thought to be under divergent selection along climatic gradients [START_REF] Gauzere | Where is the optimum? Predicting the variation of selection along climatic gradients and the adaptive value of plasticity. A case study on tree phenology[END_REF].

Co-gradient variation of the mean bud burst date has been described along these elevation gradients [START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]. We here take advantage of data collected across multiple years on the same individuals, in a common garden experiment using open-pollinated progenies, to estimate both the intercept and slope of reaction norms of bud burst date to temperature and their quantitative genetic variation within and between populations. Gene flow is very high between wind-pollinated forest tree populations [START_REF] Kremer | Historical and contemporary dynamics of adaptive differentiation in European oaks[END_REF], with pollen dispersed over very long distances [START_REF] Schueler | Viability and sunlight sensitivity of oak pollen and its implications for pollen-mediated gene flow[END_REF], Kremer et al. 2012). Due to their peculiar life history, combining high dispersal, sessile adult stages, and long lifespan, forest trees are exposed to considerable temporal and spatial variation in selection pressure. According to extant theory, we would therefore expect the evolution of optimal phenotypic plasticity to be facilitated by these specific features of tree life history, and limited genetic variation for the reaction norm of bud burst date to temperature to be maintained. Timing of bud burst is however also highly correlated with flowering time in oaks (Franjic et al. 2011), and therefore should be subject to relatively strong assortative mating. Our experimental data thus offer the opportunity to describe, for the first time, patterns of genetic differentiation for both the slope and the intercept of the reaction norm along an elevation gradient for a trait involved in assortative mating, in a situation where, in the presence of random mating, we would expect very little differentiation.

We then used individual-based simulations to investigate the evolution of phenotypic plasticity for a polygenic quantitative trait under assortative mating, in a meta-population resembling that sampled along the monitored elevation gradient. We assumed that selection on the trait value differed along the spatial gradient, but, for the sake of simplicity, ignored temporal variation in selection. We explicitly modelled gene frequency changes at loci with environment-dependent, or constant, effects on the trait affecting respectively, the slope, or the intercept of a linear reaction norm. The empirical data were used to parameterize some aspects of the model. We compared model predictions, with and without assortative mating, to observed patterns of differentiation of reaction norms along the elevation gradient. Yet, our model remains highly general and was used to explore situations with different intensities and patterns of gene flow, departing from the high gene flow observed in our empirical case study.

Methods

Variation of reaction norms for the timing of bud burst along elevation gradients in Quercus petraea

Sampled populations. Following the last deglaciation, sessile oak (Quercus petraea) has recolonized elevation gradients in Northern side of the Pyrenees Mountains, France, for approximately 300 generations [START_REF] Brewer | The spread of deciduous Quercus throughout Europe since the last glacial period[END_REF]. We sampled 10 populations from two parallel valleys (Luz and Ossau, from 42°53'N, 00°25'W to 43°45'N, 00°14'W) along a gradient extending from 131 to 1630 m above sea level (for details on the experimental set up, see Vitasse et al. 2009a, Vitasse et al. 2009b[START_REF] Niort | Population differentiation of sessile oak at the altitudinal front of migration in the French Pyrenees[END_REF][START_REF] Alberto | Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient[END_REF][START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]. Genetic differentiation at neutral markers was previously found to be low within each valley ( F SV = 0.022), as well as between valleys ( F VT = 0.003, [START_REF] Niort | Population differentiation of sessile oak at the altitudinal front of migration in the French Pyrenees[END_REF], with no significant pattern of isolation by distance [START_REF] Alberto | Imprints of natural selection along environmental gradients in phenology-related genes of Quercus petraea[END_REF]. This overall suggests that the populations of the valleys are connected by high gene flow through long distance pollen dispersal. Much higher genetic differentiation was observed at the trait level for the timing of bud burst (mean Q ST of 0.2), suggesting divergent selection between populations [START_REF] Alberto | Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient[END_REF][START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]. [START_REF] Dantec | Escape of spring frost and disease through phenological variations in oak populations along elevation gradients[END_REF] showed that early flushing enables trees to escape disease, whereas late flushing protects them against late frost damage. The strength and direction of selection on timing of bud burst might therefore be expected to vary along the elevation gradient, resulting in contrasting optimal values for bud burst date (see Gaüzere et al. 2020 for theoretical predictions). In situ, populations at high elevation flush about 52 days later than populations at low elevation [START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF].

Estimation of reaction norms.

In each population, we sampled open-pollinated progenies from individual maternal parent trees (152 parent trees in total, with 23 offspring on average per parent and a total of 2963 offspring individuals), which were planted in a common garden at sea level (Toulenne Research Station 44°34'N, 00°16'W, South Western France, see [START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]. The timing of bud burst of each tree in the common garden was recorded by scoring the trees twiceweekly during spring time and observations were repeated over five successive years (from 2010 to 2014). Daily temperature was also recorded using data loggers, (HOBO Pro RH/Temp, Onset Computer Corporation, Bourne, MA) protected from rain or direct sunlight (see Vitasse et al. 2009b for more details). Bud burst date varied as a function of spring temperature, which was measured as mean daily temperature between March 1 st and May 31 st (Vitasse et al. 2009b). Spring temperature in the common garden varied between 11.2°C and 14.3 °C across the 5 years of the survey (in the comparison, the mean spring temperature measured in situ over the same years varied between 3.9°C and 11.3°C across populations). We assumed that the timing of bud burst for each tree varied linearly with spring temperature, with intercept a and slope b. The intercept a corresponded to the bud burst timing of the tree at the mean spring temperature across the five years of survey. The differences in spring temperature between years made it possible to estimate the intercept and the slope at the single tree level.

Genetic variation for reaction norms.

Sampling open pollinated progenies allowed estimating, for both the intercept and slope of reaction norms, their total phenotypic variance among individuals, and the between-population and a within-population (V B and V W ) genetic components of this variance. We fitted a mixed random-effects regression model, to estimate the variance components of the reaction norm parameters at several levels. Variance components were estimated from the following mixed-effect model:

Z kpmij =B k + b E j + ( a p +b p E j ) + ( a pm +b pm E j ) + ( a pmi +b pmi E j ) +e kpmij (1)
where Z kpmij is the date of bud burst of individual i from the maternal family m and population p, in year j, found in block k. B k is the fixed effect of block k in the common garden, E j is the deviation of spring temperature experienced in a given year j compared to the mean spring temperature in the surveyed period, and b is the average slope of the linear response to temperature in the common garden. The deviation from this relationship was quantified by the random terms in brackets at several hierarchical levels. The levels considered were the population (subscript p), open-pollinated maternal family (subscript m) and offspring individual effects (subscript i, i.e. a permanent environmental effect, Lynch and Walsh 1998), and we also considered a final deviation e for the measurement j of an offspring i at year j (for a review of this type of random regression approaches in a similar context, see [START_REF] Gienapp | Evolutionary dynamics in response to climate change[END_REF]. The temperature E was centered on the mean temperature over the years considered, to ensure that slope variation did not have a spurious influence on the measurement of reaction norm intercept. Estimates from the two valleys (figure 1) were compared using a simple preliminary ANCOVA, regressing parameters of individual reaction norms against elevation. No significant difference in either the intercept or slope of the reaction norm was found between valleys (difference in intercept between valleys: -0.00247 day/m 95% CI:

-0.00524, 0.00032, difference in the slope of the reaction norm: 0.00005 day/degree/m 95% CI: 0.00000, 0.00011). We used the MCMCglmm R package [START_REF] Hadfield | MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R Package[END_REF]) to fit the model in (1) by a Bayesian approach. Regarding the priors, we used a zero-mean normally distributions with large variances (10 4 ) for the fixed effects, and parameter-expanded priors for the variance components, and a parameter n set to 0.002 [START_REF] Hadfield | MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R Package[END_REF]). Parameter estimates were not sensitive to changes in prior specification. Bayesian variance components estimates were very close to those obtained with restricted maximum likelihood methods (results not shown). The model was run for 110 000 iterations, including a burn-in of 10 000 iterations and a thinning interval of 100 iterations.

Autocorrelation levels in the Markov chain-Monte Carlo (MCMC) resampling were consistently below 0.07. We checked that the distribution of estimated random residuals in (1) did not show any worrying signal of strong deviations from a linear reaction norm to temperature (see figure S1).

The genetic variance between populations for the intercept of the reaction norm a was estimated as V B =var ( a p ), the variance of the population term for the intercept. The additive genetic variance within populations for the intercept was estimated as V W =4 var ( a m ) with var ( a m ) the variance among open pollinated progenies, assuming that the progenies are half-sibs (consistently with the earlier parentage analysis conducted in multiple oak stands, in Gerber et al. ( 2014)). The same variance components were estimated for the slope of the reaction norm b.

We used estimates of within-and between-population genetic variance to compute Q ST , a standardized measurement of genetic differentiation (Le Corre and Kremer 2012), separately for the intercept, and slope of the reaction norm:

Q ST = V B (V B +2 V W ) (2)
The Q ST estimates were compared to estimated F ST from neutral molecular markers. Neutral population differentiation was estimated on the basis of 16 microsatellites markers genotyped in the same 10 oak populations [START_REF] Alberto | Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient[END_REF]. According to [START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF], the average F ST value was 0.026 and the conservative 95% extreme tail of the F ST distribution was 0.047, assuming a Lewontin-Krakauer distribution (see [START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF] for details). As recommended by [START_REF] Whitlock | Evolutionary inference from QST[END_REF], we concluded that trait divergence exceeded that expected under drift-migration balance if the Q ST estimate exceeded the 95% extreme tail of the F ST distribution and we mostly focused on the comparison between the Q ST of the traits under study (i.e. the intercept and the slope of the reaction norms).

Simulation model

We used the individual-based model Metapop [START_REF] Soularue | Metapop: An individualbased model for simulating the evolution of tree populations in spatially and temporally heterogeneous landscapes[END_REF] to simulate the evolution of individual reaction norms across a spatially heterogeneous landscape, similar to the elevation gradient described above. To better understand the effect of assortative mating on the evolution of phenotypic plasticity, we focused on scenarios where, under random mating, we would expect little constraint on the evolution of an optimal reaction norm allowing to express the best phenotype in each environment. In particular, we assumed that selection varied in space, but not in time, and was perfectly predictable from environmental cues, with no intrinsic cost of plasticity.

Phenotypes and genotypes:

We assumed that some phenotypic trait (e.g. bud burst date) varied plastically as a function of some environmental cue (e.g. spring temperature), with a linear reaction norm, such that the phenotype of an individual i in the population p is:

Z pi =a i +b i E p +ϵ i (3)
where a i is the breeding value for the intercept of the reaction norm of individual i and b i is the breeding value of the slope of its reaction norm, and E p is the value of the environmental cue in population p, shared by all individuals in that population, and ϵ i is a random micro-environmental effect, drawn randomly in a Gaussian distribution with mean 0 and variance σ ϵ

2

. We assume that the slope b i and the intercept a i of the individual reaction norm are each determined by 10 freely recombining diploid loci (respectively environmentally sensitive and non-sensitive loci), without dominance or epistasis, with:

a i = ∑ l=1 10 α il f +α il m and b i = ∑ l=1 10 β il f + β il m
, where α il f is the allelic value at locus l that individual i has inherited from its father, (and α il m from its mother), which affect the intercept of its reaction norm. Similar notation (β il f and β il m ) is used for alleles affecting the slope of the reaction norm. The breeding value for the intercept of the reaction norm, a i , measures the expectation for the average trait value if offspring of our focal individual were measured in some reference environment with E p =0 . Equation (3) implies that plasticity does not affect phenotypic stability and the phenotypic variance among individuals with the same genotype in the same site.

We also simulated an additional neutral locus to contrast patterns of divergence at neutral markers and selected quantitative traits.

Natural selection:

We simulated stabilizing selection within each population on the phenotype, where the fecundity of an individual with phenotype Z in population p, was computed as:

W p ( Z)=exp [ -( Z -Z opt p ) ² 2 ω 2 ] (4)
where Z opt p is the local phenotypic optimum and ω 2 is inversely related to the intensity of stabilizing selection around that optimum. We assumed no intrinsic cost of plasticity, such that the fitness of an individual depends only on its phenotype and not directly on the slope of the reaction norm.

Landscape:

The meta-population simulated consisted of 55 patches homogeneously distributed across a 5 x 11 grid, mimicking five populations at each of 11 elevation levels (for an illustration see [START_REF] Soularue | Assortative mating and gene flow generate clinal phenological variation in trees[END_REF]. Note that this gradient could also be interpreted as varying across latitude in our relatively general model. The optimal phenotype Z opt values was assumed to vary linearly across the elevation gradient, with slope S Zopt , to simulate divergent selection. The environmental cue E p also varied linearly with elevation with slope S E . We considered that both phenotypic values and environmental values were centered such that the optimal phenotype and environmental cue would be 0 at the center of the elevation gradient. In our fictive landscape, the five populations located at the same elevation were exposed to the same macro-environmental and selective conditions. Under these assumptions, the reaction norm that allows the expression of the optimal phenotypes in all populations is characterized by an intercept a opt =0 and a slope b opt =S Zopt /S E .

Life cycle:

The model simulated non-overlapping generations for a hermaphroditic plant species with seed and pollen dispersal.

1. Fecundity selection. We assumed that the match between the phenotype of an individual and its local environment affected similarly the production of female (ovules) and male (pollen) gametes (see equation 4).

2. Gene flow. Seeds and pollen were dispersed between populations with Nm migrants per generation. Seeds were exclusively dispersed according to the stepping stone migration model [START_REF] Kimura | The stepping stone model of population structure and the decrease of genetic correlation with distance[END_REF], assuming each population received migrants from the two adjacent populations at the same elevation level and from one population at higher, and one at lower elevation. Pollen was dispersed according to, either the Wright's island migration model connecting all populations [START_REF] Wright | Evolution in Mendelian populations[END_REF], or the stepping stone migration model. The seed and pollen dispersal rates m s and m p were initialized from the Nm value, assuming m=m s +1 /2 ( m p -m p * m s ) (Lopez et al. 2008). There was no border effect when the island migration model was simulated, each population could receive pollen from other populations with equal probability. However, under the stepping stone migration model, populations at the edges of the landscapes were only connected to two or three populations, while each population from the core was connected to four populations.

We assumed a constant rate of incoming gene flow in all populations: only the number of sources for migrant pollen and seeds varied at the edge, but not the total number of migrants.

3. Recruitment. To mimic population growth during colonization, we used the following phenomenological growth model to compute the number of new recruits in each population from the number of individuals in the metapopulation in the previous generation:

N p , t =N p ,t -1 + g × I p , t × ( 1-N p , t -1 K ) (5)
where With probability τ , the offspring was generated through selfing. In case of outcrossing, the process of drawing the father of the individual depended on whether we assumed random mating or assortative mating. In the simplest case of random mating, the population of origin of the male parent was first randomly drawn from all populations contributing to the pollen cloud reaching population p, as described above for seed dispersal. Then the identity of the father was drawn with probability proportional to their fitness W k ( Z M ) in their population of origin k. In the case of assortative mating, mating could occur only if the phenotypic distance between the phenotype of the selected mother Z F and father Z M did not exceed δ, the maximal phenotypic distance between parents. The population of origin of the male parent was first randomly drawn from all populations contributing to the pollen cloud reaching the mother population, but weighting each population by the number of compatible mates for the focal mother that each contained rather than by the total number of individuals. The male parent was then drawn according to its fitness value in its patch of origin from all the individuals phenotypically similar to the female parent in that population. In the case that no male individual in the local pollen pool had a phenotype sufficiently similar to that of the female parent, another female parent in the patch was randomly drawn. The value of the maximal phenotypic distance between parents was allowed to vary between generations and, for each generation, this maximal distance was drawn in truncated Gaussian distribution centered on 0 with standard deviation σ δ . The average intensity of assortative mating was a posteriori quantified by measuring the correlation ρ between phenotypic values of male and female parents of a given offspring, within each population, and averaged across populations.

5.

Recombination. We assumed that the 21 loci (20 under selection and one neutral) were located on separate chromosomes and recombined freely.

6. Mutation. For each offspring produced, the alleles at each locus mutated with a probability μ per generation according to the k-allele mutation model [START_REF] Peng | Forward-time population genetics simulations: methods, implementation, and applications[END_REF]. Under this model, each mutation replaces the mutated allele by another of the k= 256 possible alleles at the given locus.

Initialization of simulations.

We simulated the colonization of the elevation gradient from low elevation. At the first generation, 500 founding individuals were set in each population at the lower elevation level. For each (environmentally non sensitive) a locus and each (environmentally sensitive) b locus, the 256 allelic effects were independently drawn from the Gaussian distributions

N ( a opt /(2 L a ) , σ a 2 / (2 L a ) ) or N ( b opt /(2 L b ) ,σ b 2 /(2 L b ) )
, respectively, where L a and L b are the number of loci affecting the intercept a and slope b of the reaction norm. Eighteen alleles, randomly chosen among the 256 possible alleles, were present at each locus in the starting meta-population. For each locus , we drew q the vector of initial allelic frequencies at the metapopulation level for each the 18 alleles in a Dirichlet distribution with parameter 1, thus assuming equifrequent alleles. As a result, the expected initial additive variance for respectively the intercept and slope of the reaction norm at the scale of the metapopulation were σ a 2 and σ b 2 , with cov(a, b)=0, and the mean individual reaction norm was initially optimal with intercept a opt and slope b opt . We considered the initial metapopulation to be in drift-migration equilibrium under an island model, without selection. Initial local allelic frequencies in each patch were randomly generated for all the loci as in [START_REF] Goudet | The effects of dominance, regular inbreeding and sampling design on QST, an estimator of population differentiation for quantitative traits[END_REF]. At each locus, the allelic frequencies followed a multinomial Dirichlet distribution with parameters (55,4 Nmq ) where q is the allelic frequency vector generated for the overall meta- population and 55 is the number of populations in the metapopulation [START_REF] Gaggiotti | Quantifying population structure using the F-model[END_REF].

The initial individual genotypes were then randomly assigned from the initial allelic frequencies generated, assuming no initial linkage disequilibrium and Hardy-Weinberg equilibrium within each population.

The initial values for σ a 2 and σ b 2 , the additive genetic variance for respectively the intercept and slope of the reaction norm, were set by considering the initial proportion of phenotypic variance due to the variance at loci a and b ( ϕ a and ϕ bE as described by [START_REF] Lande | Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilation[END_REF]) in a site with a given reference environmental value E 0 with:

ϕ a = σ a 2 σ a 2 +σ b 2 E 0 2 + σ ϵ 2 (6) ϕ bE = σ b 2 E 0 2 σ a 2 +σ b 2 E 0 2 + σ ϵ 2 (7)
We fixed the values of ϕ a and ϕ bE and computed the initial additive variance σ a 2 as:

σ a 2 = ϕ a σ ϵ 2 2× n A × ( 1 -ϕ a -ϕ bE ) (8)
Similarly, the initial additive variance σ b 2 was computed as:

σ b 2 = ϕ bE σ ϵ 2 2× n B × E 0 2 ( 1 -ϕ a -ϕ bE ) (9)
The reference environment cannot be the center of the gradient, at which, following our convention

E 0 =0
, and ϕ bE =0 whatever the variance for the slope of the reaction norm. We therefore set an arbitrary non null value for E 0 in the reference environment and verified that this choice had no major impact on our main conclusions.

Simulations were run for 10 000 generations. A visual inspection of simulation results suggests that the metapopulation has then reached a stationary state for most scenarios. We also kept track of the transient dynamics in the scenarios simulated. Fifty independent replicates were simulated for each parameter set.

Monitoring the evolution of phenotypic plasticity in simulations

We examined how local reaction norms differentiated along the elevation gradient and differed from the optimal reaction norm: we recorded the difference between the local mean breeding value for the intercept ā (or slope b) in each population p and the optimal intercept a opt (or optimal slope b opt ) and plotted this difference as a function of elevation. We also examined how the mean phenotype expressed in each population deviated from the local optimal phenotype and how genetic divergence for the intercept of the reaction norm contributed to total phenotypic variation along the environmental gradient. We in particular compared the slopes of the genetic and phenotypic clines in the simulations, as in the oak populations experimentally studied. As in the empirical study, we calculated Q ST values for the intercept ( Q STa ) and slope ( Q STb ) of the reaction norm separately.

Variance components were not estimated as in the empirical study, but were instead measured directly in the simulations, on the basis of the individual additive values available for each locus.

F ST values were estimated in the simulations with Genepop [START_REF] Rousset | genepop'007: a complete re-implementation of the genepop software for Windows and Linux[END_REF]) according to the method of [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF], as in the oak populations studied. Lastly, we monitored the composition and the intensity of effective gene flow among elevation levels, with an aim of identifying the drivers of plasticity evolution under assortative and random mating to help the interpretation of our results.

Choice of parameter values:

We parameterized the model with available data about the timing of bud burst in sessile oak, but also explored the parameter space for variables associated with the known evolutionary drivers of phenotypic plasticity evolution, with a specific focus on the pattern and level of gene flow (table 1). By default, we assumed moderate strength of selection, of the same magnitude than the average strength of stabilizing selection on bud set date estimated by [START_REF] Gauzere | Where is the optimum? Predicting the variation of selection along climatic gradients and the adaptive value of plasticity. A case study on tree phenology[END_REF] in oak populations in French Pyrenees, but we also assessed the effects of a stronger selection pressure. We set our default value for Nm to 10, a value typical of broad-leaved species [START_REF] Austerlitz | Effects of Colonization Processes on Genetic Diversity: Differences Between Annual Plants and Tree Species[END_REF], consistently with pairwise F ST in the studied oak populations [START_REF] Niort | Population differentiation of sessile oak at the altitudinal front of migration in the French Pyrenees[END_REF][START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]), but also considered much lower (Nm = 1) and higher (Nm = 100) migration. We set a much higher level of gene flow through pollen than seeds with m p =100 * m s in agreement with experimental data for many broad-leaved and tree species [START_REF] Ennos | Estimating the relative rates of pollen and seed migration among plant populations[END_REF]). The island model of migration for pollen dispersal fits better than the stepping-stone model the biology of wind pollinated tree species, for which considerable gene flow occurs over large distances (Bushbom 2011, Kremer et al. 2012[START_REF] Gerber | High rates of gene flow by pollen and seed in oak populations across Europe[END_REF], and is consistent with the limited isolation by distance observed in the studied populations [START_REF] Alberto | Imprints of natural selection along environmental gradients in phenology-related genes of Quercus petraea[END_REF]. By contrast the stepping stone migration model gave us the opportunity to assess the influence of assortative mating on the evolution of plasticity when gene flow though pollen is limited in distance. A low rate of selfing (τ = 0.02) was considered, as observed in oak species [START_REF] Bacilieri | Mating system and asymmetric hybridization in a mixed stand of European oaks[END_REF]. Different strengths of assortative mating were simulated by varying σ δ , resulting, for our default parameter set, in strong ( ρ = 0.62) or moderate (ρ = 0.27) correlation between parents' phenotypes (see figure S2). These two values are in the range of phenotypic correlations estimated by Weiss (2014) in different plant species (ρ = 0.046 to 0.624).

We set the initial values for ϕ a and ϕ bE so that genetic variation in the slope of the reaction norm was large and plasticity contributed largely to phenotypic variance away from the central population. We varied the ϕ a and ϕ bE values (and thus the initial relative evolutionary potential for the intercept and slope of the reaction norm) to examine whether our conclusions depended on our assumptions on the genetic architecture of the trait. The strength of divergent selection in space was also varied by changing the slope S Zopt .

Results

Estimation of reaction norms for the timing of bud burst in natural populations of sessile oaks

The estimated mean intercept of reaction norms for the timing of bud burst for Q. petraea was 116.324 days (95% credible interval (CI): 113.809, 118.841 days), corresponding to the Julian date of bud burst at the mean spring temperature. The estimated regression slope against annual spring temperature was -2.398 days/°C (CI: -2.546, -2.230 days/°C). Thus, a spring temperature lower by one degree postpones the timing of bud burst by more than two days. We found a significant variation between populations in the intercept of the reaction norm, following a co-gradient pattern:

trees from higher elevations opened their buds in average later than trees from lower elevation in the common garden (figure 1), as they did in situ. The differentiation for the intercept of the reaction norm (

Co-gradient variation evolves under assortative mating in the simulations

With random mating, the intercept of the reaction norm a remained in most cases very close to the optimal value a opt in all populations (see the black lines in the left panels of figures 2-3). Cogradient genetic variation across the environmental gradient emerged in our simulations with random mating only when pollen flow was strongly limited both in intensity and distance (e.g. with stepping stone dispersal in figure 4, see table 3 for a summary of simulation results). Even then, genetic clines were very shallow, faded with time (see figures S3, S4 and S5) and were largely blurred by stochastic variation caused by strong drift (see the large variation between replicates in figure 5). Differentiation between populations for the intercept of the reaction norm was not greater than that expected for neutral markers with random mating ( Q STa mostly lower than neutral F ST in figures 5 and 6), contrasting with our empirical observations in the oak populations (table 1).

By contrast, in all scenarios with assortative mating, co-gradient variation for the intercept of reaction norms evolved at equilibrium across the environmental gradient (see the blue and red lines in the left panels of figures 2,3,4; table 3), including cases where gene flow through pollen was strong and not distance-limited, as in the studied oak meta-population. Genetic differentiation for the intercept of the reaction norm exceeded that expected under neutrality in most cases with assortative mating (see figures 5 and 6), especially when gene flow was strong, consistently with the empirical patterns in oaks. Interestingly, with assortative mating, genetic divergence for the intercept of the reaction norm was transiently stronger after the colonization of the gradient (see figures S3-S8, which show that the gradient was fully colonized after approximately 200 generations). The predicted Q ST values for the intercept were quantitatively comparable to that estimated in the oak populations when approximately the same number of generations had passed since the colonization of the elevation gradient by oaks (i.e. 500 generations see figure 3 and table 2). Genetic clines that evolved under assortative mating were steeper with higher strength of assortative mating (compare blue and red lines in left panels of figures 2 and 4) and higher intensity of divergent selection along the environmental gradient (compare figures 2 to S9), but did not depend on the initial genetic variance in plasticity (figure S10). Phenotypic plasticity however explained most of the variation of the trait along the environmental gradient in all scenarios examined, with genetic clines contributing only marginally (see figures 3C, 3D, S11-S12).

Under assortative mating, the slope of the evolved reaction norm in the simulations is either steeper or shallower than optimal

The slope of reaction norm b remained close to optimal value b opt at equilibrium when the individuals mated randomly, with deviations consistent with the effect of drift when gene flow was very limited (figures 2 and 4, black lines in right panels, table 3). Optimal plasticity evolved regardless of the strength of stabilizing selection, the slope of phenotypic optima, the slope of the environmental gradient, the reference environment, or the level of contribution of environmentally sensitive loci to genetic differentiation (e.g. black lines in figure S9).

In contrast, with assortative mating, the slope of reaction norm deviated in general from the optimal slope (blue and red lines in figures 2, 3, 4, S9; see summary in table 3). The extent and direction of this deviation however depended on time since colonization of the gradient, the exact pattern and intensity of gene flow, the position within the range (core vs margins) and the strength of assortative mating, among other things (see also figures S3 to S8). We found both cases where plasticity was suboptimal (slope too shallow) and cases of hyperplasticity (slope too steep). With the island migration model, we often found that assortative mating decreased the level of plasticity near equilibrium compared to random mating (figure 2, colored lines in the right panels), except sometimes at the margins of the landscape (but see figure S9). With the stepping stone migration model, the slope of reaction norms sometimes evolved towards higher levels than optimal (figure 4, colored lines in the right panels). In most scenarios, deviations of the slope of the reaction norm from its optimal value were much stronger shortly after the colonization of the environmental gradient, with typically shallower slopes (figure 3 right panels) than at equilibrium (figures S3 to S8). In some cases, the slope was initially too shallow and progressively evolve to be steeper than optimal in the long term (e.g. figure S7).

In scenarios with assortative mating, the slope of the reaction norm tended to be steeper at the margins of the environmental gradient than at the core, while patterns with random mating were less clear and varied with the intensity and distance of gene flow (see figure 2 and4). As in the oak meta-population studied experimentally (figure 1,table 2), the level of differentiation among populations for the slope of reaction norms in our simulations, as estimated by Q ST , was systematically lower than the level of genetic differentiation for the intercept in the presence of assortative mating (figure 5 and6). In our reference parameter set (island model of dispersal with Nm=10), which best matches our empirical case study, the differentiation for the slope was not significantly greater than expected under neutrality, as in the oak data (figure 5, table 2). Yet, with higher intensity of gene flow (Nm=100, figure 5), or shorter dispersal distance (stepping stone, figure 6), Q ST for the slope of the reaction norm often exceeded the neutral F ST .

As a consequence of the combined evolution of non-optimal slope of reaction norms and cogradient variation for the intercept under assortative mating, the mean phenotypes expressed in peripheral populations slightly deviated from the local optima, especially when the strength of assortative mating was intermediate, gene flow high, or few generations had elapsed after colonization of the gradient, in contrast with scenarios of random mating where the optimal phenotype was always expressed in all environments (figures 3C, 3D, S11 and S12).

Discussion

Genetic variation for reaction norms along an elevation gradient in Q. petraea

Previous studies of quantitative genetic variation for bud burst timing in the same oak populations in the French Pyrenees [START_REF] Alberto | Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient[END_REF][START_REF] Alberto | Imprints of natural selection along environmental gradients in phenology-related genes of Quercus petraea[END_REF][START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]) have shown that trees from lower elevation open their buds earlier in the spring than trees from higher elevation, both in situ and when experiencing the same temperature in a common garden, a pattern of genetic divergence described as co-gradient variation. Similar patterns have also been observed at larger geographic scales, and in other species [START_REF] Alberto | Imprints of natural selection along environmental gradients in phenology-related genes of Quercus petraea[END_REF][START_REF] Robson | Flushing phenology and fitness of European beech (Fagus sylvatica L.) provenances from a trial in La Rioja, Spain, segregate according to their climate of origin[END_REF]Muller et al. 2015[START_REF] Sampaio | Does the geography of cork oak origin influence budburst and leaf pest damage?[END_REF]). Our study builds on this previous work by investigating separately the genetic variation for the intercept and slope of the reaction norm of bud burst timing to temperature. Since the intercept of these reaction norms measures the expected bud burst date at the average temperature in the studied period, our results logically confirm past findings about the variation of average bud burst dates (in particular those in [START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]): we found significant co-gradient genetic divergence among populations across elevation for the intercept, with standardized measures of differentiation as measured by Q ST exceeding expectations under neutrality ( F ST ). In contrast, our analysis shows little or no genetic variation for the slope of these reaction norms, both within and between populations, and Q ST measures for the slope not exceeding the low differentiation at neutral markers as measured by F ST . To better understand these contrasting patterns of genetic divergence for the slope and intercept of the reaction norm of plastic traits along environmental gradients, we have used a theoretical model, in which both traits are free to evolve.

Simulations show that assortative mating affects the evolution of plasticity in heterogeneous environments

The observed pattern of co-gradient variation for the intercepts of reaction norms is not consistent with the theoretical expectation that generalist genotypes with optimal reaction norms should prevail in heterogeneous habitats in the presence of strong gene flow over large distances [START_REF] Tufto | The Evolution of Plasticity and Nonplastic Spatial and Temporal Adaptations in the Presence of Imperfect Environmental Cues[END_REF][START_REF] Scheiner | The genetics of phenotypic plasticity. XII. Temporal and spatial heterogeneity[END_REF][START_REF] Schmid | Species' range dynamics affect the evolution of spatial variation in plasticity under environmental change[END_REF], as occurs in this oak meta-population. As mating is more likely to occur between trees with similar timing of bud burst, we investigated the possible role of assortative mating in the paradoxical maintenance of genetic differentiation for reaction norms despite abundant gene flow over large distances. We acknowledge that other factors, present in the studied oak populations and ignored in our model, such as temporal variation in selection (e.g. King and Hadfield 2019), overlapping generations and non-linear responses of phenology to temperature [START_REF] Wolkovich | A simple explanation for declining temperature sensitivity with warming[END_REF], are likely to affect the evolution of phenotypic plasticity in those

populations. Yet, the present simulations shed light on the rarely considered effect of assortative mating on the evolution of plasticity in spatially heterogeneous environments. The general expectation is that plasticity then evolves under the dual influence of sexual selection due to assortative mating, favoring resemblance to most abundant and fit females, and natural selection, favoring expression of phenotypes close to the local optimal phenotype [START_REF] Nonaka | Assortative mating can limit the evolution of phenotypic plasticity[END_REF].

Without gene flow, there is no conflict: the most abundant and fit females and males have a phenotype matching the local optimum. But, with the opportunity to mate with a female outside the site where the male developed, there may be a conflict between phenotypes increasing overall mating success and phenotypes maximizing survival or gamete production. [START_REF] Stam | The evolution of reproductive isolation in closely adjacent plant populations through differential flowering time[END_REF], and Soularue andKremer (2012, 2014) have already shown that assortative mating can shape cogradient differentiation patterns along environmental gradients, and oppose the effect of natural selection on genetic divergence. In their simulations, the slope of the reaction norm however did not vary among individuals, and could not evolve. The present simulations show that, in the presence of assortative mating and strong pollen flow, a co-gradient pattern of genetic divergence for the intercept of the reaction norm still evolves when both the intercept and the slope of the reaction norms are genetically variable. We further found that, with assortative mating, the average slope of the reaction norm deviates from its optimal value under random mating, being either too shallow or too steep depending on the strength assortative mating and patterns of gene flow (see table 3 for a summary of our main results). In the following, we discuss each of these findings in greater detail and compare our theoretical and empirical results.

Assortative mating leads to the evolution of co-gradient variation through a process of genetic spatial sorting when the reaction norm is free to evolve

As in [START_REF] Stam | The evolution of reproductive isolation in closely adjacent plant populations through differential flowering time[END_REF] and Soularue andKremer (2012, 2014), the evolution of co-gradient variation for the intercept of the reaction norm with assortative mating can be explained by a process of spatial genetic sorting, mediated by sexual selection. Assortative mating leads to biased pollen immigration of early flowering genotypes in environments where plasticity causes early flowering, resulting in a spatial reorganization, across the environmental gradient, of the alleles determining the intercept of the reaction norms. In particular, with assortative mating, co-gradient variation for the intercept evolves even when the slope of the reaction norm is optimal or close to optimal, and not only when plasticity is partially adaptive as found under random mating [START_REF] Schmid | The role of phenotypic plasticity on population differentiation[END_REF]. The process of genetic sorting predicts stronger co-gradient genetic divergence for the intercept of the reaction norm when assortative mating is stronger [START_REF] Stam | The evolution of reproductive isolation in closely adjacent plant populations through differential flowering time[END_REF]), which agrees with our simulation findings. We therefore expect species with strong assortative mating to display more often cogradient variation than species where assortative mating is weaker; we also expect in the same species, traits involved in assortative mating to display co-gradient variation more often than traits for which mating is random. Our model however does not explain the evolution of counter-gradient variation, as observed in other species (e.g. in beech, [START_REF] Chmura | Variability of beech provenances in spring and autumn phenology[END_REF]Rozkowski 2002, Gömöry and[START_REF] Gömöry | Trade-off between height growth and spring flushing in common beech (Fagus sylvatica L.)[END_REF], as this pattern never emerged in our simulations, neither with random nor assortative mating, when the reaction norm was free to evolve.

Despite the frequent evolution of genetic clines for plastic traits under assortative mating, our model still predicted that plasticity explained most of phenotypic variation along environmental gradients when the evolution of the reaction norm was not constrained by plasticity cost. This prediction is consistent with our empirical findings (e.g. there is a factor 10 between the slope of the phenotypic and genetic cline in bud burst date along elevation in the studied oak populations, Firmat et al.

2017

) and with general empirical observations indicating that plasticity strongly contributes to phenotypic variation in trees [START_REF] Vitasse | Quantifying phenological plasticity to temperature in two temperate tree species[END_REF], Kremer et al. 2014) and other taxa [START_REF] Merilä | Climate change, adaptation, and phenotypic plasticity: the problem and the evidence[END_REF].

Assortative mating leads to the evolution of non optimal plasticity

Our model predicts the evolution of, either suboptimal plasticity (reaction norms with a slope shallower than optimal), or hyperplasticity (slopes steeper than optimal) in the presence of assortative mating, in circumstances where optimal plasticity allowing to express optimal phenotypes in all environments would evolve under random mating. The former result is congruent with the findings of [START_REF] Nonaka | Assortative mating can limit the evolution of phenotypic plasticity[END_REF] who predict, with a two-patch model and strong gene flow, a suboptimal level of plasticity with intermediate assortative mating. The latter result is, to our knowledge, more original. Few models predict the evolution of hyperplasticity [START_REF] Scheiner | The genetics of phenotypic plasticity. XII. Temporal and spatial heterogeneity[END_REF][START_REF] King | The evolution of phenotypic plasticity when environments fluctuate in time and space[END_REF], and none of them involve assortative mating. Instead, explanations for the evolution of hyperplasticty in these models involve either different patterns of temporal and spatial variation [START_REF] King | The evolution of phenotypic plasticity when environments fluctuate in time and space[END_REF], or the benefits of increased variance in populations very far from the optimal phenotype [START_REF] Scheiner | The genetics of phenotypic plasticity. XII. Temporal and spatial heterogeneity[END_REF], see discussion in Tufto 2015), none of which applies in the present model. Interestingly, our model predicts the joint evolution of both hyperplasticity and co-gradient genetic divergence in the presence of assortative mating, while [START_REF] Schmid | The role of phenotypic plasticity on population differentiation[END_REF] predict that, under random mating, counter-gradient genetic divergence should evolve with hyperplasticity. Counter-gradient divergence is even often used to infer the existence of hyperplastic reaction norms [START_REF] Grether | Environmental change, phenotypic plasticity, and genetic compensation[END_REF], Conover and Schulz 1995[START_REF] Tansey | Estimating the ability of plants to plastically track temperature-mediated shifts in the spring phenological optimum[END_REF]). In the presence of assortative mating (as is commonly the case with phenological traits), our model suggests caution with such inference as co-gradient (and not counter-gradient) variation can also jointly evolve with hyperplasticity. Interestingly, in a study predicting the optimal bud burst date along elevation in the same oak populations as studied here, [START_REF] Gauzere | Where is the optimum? Predicting the variation of selection along climatic gradients and the adaptive value of plasticity. A case study on tree phenology[END_REF] found that the optimal reaction norm should be much shallower than the observed reaction norm in oaks: our model allows formulating the hypothesis that assortative mating could be involved in the joint evolution of hyperplasticity and co-gradient genetic divergence in these populations.

Previous models have similarly predicted local variation in plasticity depending on the position within the range [START_REF] Chevin | Adaptation to marginal habitats by evolution of increased phenotypic plasticity[END_REF][START_REF] Schmid | Species' range dynamics affect the evolution of spatial variation in plasticity under environmental change[END_REF]): these models assumed random mating and invoked (i) asymmetric gene flow at the margins, (ii) the common assumption that phenotypes are canalized at the core and plasticity there makes little contribution to phenotypic variation, and

(iii) the history of environmental change in dynamic ranges, to explain the spatial distribution of variation in plasticity. These explanations may apply too in the present model, but the fact that the contrast between core and margins is more pronounced with assortative mating suggests that nonrandom gene flow and sexual selection may also be involved in shaping these patterns. Despite differences in plasticity between core and margins, the level of differentiation among populations for the slope of reaction norms in our simulations was systematically lower than the level of genetic differentiation for the intercept, as observed in our empirical case study.

With assortative mating, interactions between sexual and natural selection shape the evolution of the slope of reaction norm

We discuss four, non-exclusive, mechanisms through which assortative mating can make the slope of reaction norms diverge from the optimal value that allows expressing the optimal phenotype in all environments with random mating. All those mechanisms potentially play out together and interact in all scenarios. While it remains at this stage difficult to disentangle their effects in the simulations, each of them can either favor increased or decreased plasticity and their relative contribution varies depending on the extent and distance of gene flow, and the relative strength of natural and sexual selection.

(i) Assortative mating affects the intensity and distance of effective gene flow. Reduced gene flow or shorter dispersal distance result in organisms experiencing more similar environments across generations, diminishing the advantage of plasticity, which was shown to prevent the evolution of optimal plasticity in several models (e.g. Bradshaw 1965[START_REF] Via | Genotype-environment interaction and the evolution of phenotypic plasticity[END_REF][START_REF] De | The influence of life-history differences on the evolution of reaction norms[END_REF][START_REF] Sultan | meta-population structure favors plasticity over local adaptation[END_REF][START_REF] Scheiner | The genetics of phenotypic plasticity. X. Variation versus uncertainty[END_REF]. We found that assortative mating decreases the distance of effective gene flow in the island migration model, in comparison with random mating. Indeed, as phenotypic divergence increases with distance, the probability of mating between individuals with similar phenotypes decreases with distance (figure S13). With very weak intensity of gene flow, we however found some situations where F ST at neutral markers was lower with assortative mating than with random mating, suggesting higher effective gene flow in the presence of assortative mating (Nm=1 in figures 5-6). A possible explanation for this surprising result may be that, with strong assortative mating, recent immigrant individuals with a phenotype deviating from the local optimum may have more difficulty finding a matching mate in their new local population than in neighboring populations, favoring more immigrant pollen (figure S14).

(ii) Co-gradient variation of the intercept of reaction norms generated by assortative mating constrains the slope of reaction norms to suboptimal levels. Given the existence of co-gradient clinal variation for the intercept of the reaction norm, natural selection should favor genotypes expressing a lower plasticity level to reach the optimal phenotype. As assortative mating generated clines of reaction norms intercept in all scenarios, through the spatial sorting effect described above, this could contribute to explain the lower plasticity level that evolves in some scenarios. Yet, this mechanism fails to explain the evolution of hyperplasticity in some other scenarios despite co-gradient variation. The evolved phenotypic clines often deviate slightly from the optimal value with assortative mating, suggesting that the evolution of reaction norms is not shaped only by natural selection.

(iii) Sexual selection can select against highly plastic genotypes. Overall, under assortative mating, locally maladapted individuals are more likely to mate with distant females (figure S13), and conversely locally adapted males are unlikely to sire offspring in sites very distant from theirs.

Moreover, if males with a steep reaction norm manage to mate with distant females because their own phenotype in their original environment matches that of the distant female, it is unlikely that their offspring, born in a different environment from their father, will express that same phenotype.

Such progeny may therefore be counter-selected. By comparison, zygotes generated from mating involving distant male parents with low plasticity and constant phenotypes across environments are more likely to display the optimal trait in recipient populations. Interestingly, as Nonaka et al.

(2014), we found that intermediate levels of assortative mating resulted in the lowest level of plasticity. As in their model, this may be because strong assortative mating drastically decreased the amount of distant gene flow, making mating between individuals in different environments very unlikely.

(iv) Spatial sorting, generated by assortative mating, also modifies the distribution of alleles with effects on the slope of the reaction norm. With assortative mating, the same process of non-random gene flow as we described to explain the evolution of co-gradient variation for the intercept of the reaction norm also acts on the slope, but with a different pattern. Genotypes with high plasticity tend to accumulate near the margins of the range, while genotypes with low plasticity preferentially accumulate near the core. This is due to biased immigration of most plastic genotypes towards the margins of the gradient, and biased immigration of least plastic genotypes towards the core. Indeed, males with high plasticity expressing a more extreme phenotype than what is locally optimal have a siring advantage when sending their pollen in adjacent populations closer to the margins and with more extreme environments. Conversely, males with a plasticity lower than optimal have a siring advantage when dispersing their pollen to adjacent populations closer to the core with less extreme environment. This phenomenon of spatial sorting mediated by assortative mating and the sexual selection it generates plays out differently with the stepping stone and island migration model of pollen: with stepping-stone migration, preferential immigration of more plastic genotypes from neighbouring populations closer to the core is compensated by immigration of less plastic genotypes from neighbouring populations closer to the margins, except at the range margins, where hyperplastic genotypes accumulate. With the island model, even at the margins, immigrant pollen contains a mixture of highly plastic genotypes from neighbouring populations and less plastic genotypes from further away.

Consequences of deviations from the optimal reaction norms

We found that assortative mating favored the evolution of non-optimal reaction norm, with phenotypes deviating from the optimal phenotype favored by natural selection, as a consequence of sexual selection on the trait in heterogeneous environment. These deviations were the highest during the transient regimes, which may be very relevant for populations of long-lived organisms such as trees that may not yet be at equilibrium. These deviations due to sexual selection were also the most visible at the margins of the landscape, adding to the migration and drift load expected in small marginal populations [START_REF] Alleaume-Benharira | Geographical patterns of adaptation within a species' range: interactions between drift and gene flow[END_REF]. Too shallow or too steep reaction norms may not allow the populations to respond adequately to a fast change in local environment, as occurs for instance with climate change [START_REF] Chevin | Phenotypic plasticity and evolutionary demographic responses to climate change: taking theory out to the field[END_REF][START_REF] Scheiner | The genetics of phenotypic plasticity. XV. Genetic assimilation, the Baldwin effect, and evolutionary rescue[END_REF]. Overall, the slope of the reaction norm however remained close to its optimal value even with assortative mating in our simulations, which suggests that plasticity would still help in a changing environment. Furthermore, the maintenance of genetic divergence between populations thanks to the process of spatial sorting mediated by assortative mating may be a source of adaptive variation if the reaction norm must evolve fast in a changing environment (see also [START_REF] Godineau | Assortative mating can help adaptation of flowering time to a changing climate: Insights from a polygenic model[END_REF].

Conclusions

While our model was strongly inspired by the case of bud burst date variation along climatic gradient in plants with high pollen dispersal, our conclusions about the effect of assortative mating on the evolution of phenotypic plasticity potentially apply to other traits and organisms.

Phenological traits in particular are often both plastic and involved in assortative mating, including for instance migration dates in many animals [START_REF] Charmantier | Climate change and timing of avian breeding and migration: evolutionary versus plastic changes[END_REF]. Two critical ingredients of our model are however that (i) dispersal is male-biased, (ii) the phenotypes of mates have been determined in different environments (e.g. by environmental conditions experienced before gathering in the same mating site). Under such assumptions, we expect the evolution of cogradient variation for the intercept of the reaction norm and weaker divergence for the slope, consistently with our observations in oak populations along an elevation gradient. We further expect assortative mating to prevent the evolution of optimal plasticity for adaptive traits across spatial environmental gradients as a consequence of conflicts between natural and sexual selection. Our

  model provides novel explanations for the evolution of hyperplasticity, a common empirical finding with little theoretical support. Further investigations are required to disentangle the multiple effects of assortative mating on gene flow and selection and their interactions on shaping the genetic variability of the intercept and the slope of reaction norms, and to assess how the present conclusions generalize for different life cycles, sex-specific dispersal patterns and types of plasticity.

Figure 1 .

 1 Figure1. Estimations of mean intercept (ā) and slope ( b) of reaction norms for the timing of bud burst in response to temperature variation in sessile oak trees (Q. petraea) as a function of the elevation of origin of the populations, measured in a common garden. Circles and squares are the mean values estimated in provenances from Gave valley and Ossau valley, respectively, each associated bar area is the 95% confidence interval of the estimation. Black lines result from subsequent linear regressions predicting the relationship between these estimations and the elevation of origin of the populations sampled. The regression slope for the intercept of reaction norms was 3.202 × 10 -3 day/m (95% CI: 1.054 × 10 -3 , 5.349 × 10 -3 ). The regression slope for the mean slope of reaction norms -2.182 × 10 -5 day/degree/m (95% CI: -5.456 × 10 -5 , 1.092 × 10 -5 ).

Figure 2 .

 2 Figure 2. Mean additive values, at generation 10 000, for the intercept (ā) and the slope ( b) of reaction norms in the presence of varying level of gene flow (Nm = 1, Nm = 10, Nm = 100) over long distance (island migration model). Red and blue lines correspond to the mean values obtained under strong (ρ = 0.62) and moderate assortative mating (ρ = 0.27). The intensity of stabilizing selection was moderate (ω 2 = 50). The slopes of the environmental gradients ( S E ) and the optimal phenotypes ( S Zopt ) were 0.5 and 0.7, respectively. Black lines are the mean values resulting from random mating (ρ = 0). Each point is the mean mean of 50 independent simulation replicates, and the 95% confidence interval of the values averaged.

Figure 3 .

 3 Figure 3. Panels A and B: mean additive values for the intercept (ā) and the slope ( b) of reaction norms at generation 500 in the presence of strong gene flow (Nm = 10) over long distance (island migration model). Red and blue lines correspond to the values obtained under strong (ρ = 0.62) and moderate assortative mating (ρ = 0.27). Black lines are the mean values resulting from random mating (ρ = 0). The intensity of stabilizing selection was moderate ( ω 2 = 50). The slopes of the environmental gradients ( S E ) and the optimal phenotypes ( S Zopt ) were 0.5 and 0.7, respectively. Each line is the mean of 50 independent simulation replicates. Error bars are the 95% confidence interval of the values averaged. Panels C and D: contribution of the differentiation at the intercept on the phenotypic differentiation across the environmental gradient (change in mean phenotypic value Z along elevation) under random (ρ = 0, panel C) and strong assortative mating (ρ = 0.62, panel D). Stars correspond to the mean phenotypic value expressed at each elevation level. Blue dashed lines represents the variation in the optimal phenotypic values. Dots are the mean intercept values at each elevation, the dotted lines being the linear regressions through these values. Dark shaded areas represent the proportion of phenotypic differentiation along the gradient (light shaded area) due to genetic differentiation for the intercept of reaction norm (a loci).The ratio between the slope of the intercept predicted against the elevation gradient and the slope of phenotypes predicted against the elevation gradient was 1.61e-4 under random mating (C) and 2e-1 under assortative mating (D).

Figure 4 .

 4 Figure 4. Mean additive values, at generation 10 000, for the intercept (ā) and the slope ( b) of reaction norms in the presence of varying level of gene flow (Nm = 1, Nm = 10, Nm = 100) between close populations only (stepping stone migration model). Red and blue lines correspond to the mean values obtained under strong (ρ = 0.62) and moderate assortative mating (ρ = 0.27). Black lines are the mean values resulting from random mating (ρ = 0). The intensity of stabilizing selection was moderate (ω 2 = 50). The slopes of the environmental gradients ( S E ) and the optimal phenotypes ( S Zopt ) were 0.5 and 0.7, respectively. Each point is the mean of 50 independent simulation replicates. Error bars are the 95% confidence interval of the values averaged.

Figure 5 .

 5 Figure 5. Genetic differentiation, under the island migration model at generation 10 000, for the intercept (a) and the slope (b) of the reaction norms ( Q ST ), and at neutral locus ( F ST ). Dark and light bars correspond to the values obtained under strong assortative mating ( ρ=0.62) and random mating (ρ=0), respectively. The intensity of stabilizing selection was moderate ( ω 2 = 50). The slopes of the environmental gradients ( S E ) and the optimal phenotypes ( S Zopt ) were 0.5 and 0.7, respectively. Each bar/line is the mean of 50 independent replicated simulations. Error bars and envelopes are the 95% confidence interval of the values averaged.

Figure 6 .

 6 Figure 6. Genetic differentiation, under the stepping stone migration model at generation 10 000, for the intercept (a) and the slope (b) of the reaction norms ( Q ST ), and at neutral locus ( F ST ). Dark and Light bars correspond to the values obtained under strong (ρ=0.62) assortative mating and random mating (ρ=0), respectively. The intensity of stabilizing selection was moderate (ω 2 = 50). The slopes of the environmental gradients ( S E ) and the optimal phenotypes ( S Zopt ) were 0.5 and 0.7, respectively. Each bar/line is the mean of 50 independent replicated simulations. Error bars and envelopes are the 95% confidence interval of the values averaged

  N p , t is the number of recruited offspring at generation t in the focal population p, N p , t -1 is the focal population size at generation t-1, I p ,t , the number of parents in other populations sending seeds reaching the patch at generation t, K , the carrying capacity of the patch, and g scales the speed of population growth. This equation implies that every established individual was replaced, but establishment of additional individuals depended on seed dispersal and declined with the N i , t -1 were m s is the migration rate of seeds and N i ,t -1 the population size in the patch i in the parental generation. The actual number of individuals in a patch was rounded to, either the integer part of (5), or to the next integer, by adding to (5) a For each new recruit in a patch, the population of origin of the mother was drawn with probability 1 -m s as the local population, and with probability m s N p ,t -1 / ∑

	n
	i=1 m s * random variable drawn in a uniform distribution U(0, 1) and discarding the decimal part of the
	resulting sum. Note that equation (5) does not allow populations to decline.
	n 4. Reproduction. i=1	N i ,t -1 as the

number of already established individuals, as did Le

[START_REF] Le Corre | Colonization with long-distance seed dispersal and genetic structure of maternally inherited genes in forest trees: a simulation study[END_REF]

. Each generation, for a given patch, the number of parents of immigrant seeds I p ,t originating from n other connected populations was calculated from I p ,t = ∑ connected population p. For each reproduction event within a population p, the female parent was first randomly drawn from the individuals of population p, weighted by their fitness values W p ( Z F ).

  Q STa = 0.173, CI: 0.045, 0.397, see table2) was above the conservative 95% threshold expected under drift ( F ST = 0.047;[START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]) and only 0.07% of the MCMC posterior distribution was below this threshold. In contrast, no significant differentiation among populations was found for the slope of the reaction norm to temperature, with no clear trend as a function of elevation (figure1, table 2). The variance of the reaction norm slope was very low both within and between populations (see table2). The low Q ST found for this parameter ( Q STb = 0.089,

CI: 0, 0.907) was also poorly estimated and clearly encompassed the conservative drift threshold ( F ST = 0.047;

[START_REF] Firmat | Evolutionary dynamics of the leaf phenological cycle in an oak meta-population along an elevation gradient[END_REF]

. Contrasting with Q STa , a large fraction (38.8%) of the MCMC posterior distribution for Q STb lies below this threshold (i.e. in the narrow interval [0, 0.047]).

strong assortative mating). The intensity of stabilizing selection was moderate (ω 2 = 50). The slopes of the environmental gradients ( S E ) and the optimal phenotypes ( S Zopt ) were 0.5 and 0.7, respectively. The plastic trait both undergoes divergent selection and affects the likelihood of assortative mating.