Olivier David
email: olivier.david@inrae.fr

Gaëlle Van Frank

Isabelle Goldringer

Pierre Rivière

Michel Turbet Delof

Bayesian inference of natural selection from spatiotemporal phenotypic data

Keywords: Adaptation, Evolution, Quantitative genetics, Statistics, Wheat

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Selection is the process by which favourable types become more frequent in successive generations (Barton et al., 2007, p. 457). Selective pressures may vary in space and in time [START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF][START_REF] Papaïx | Influence of cultivated landscape composition on variety resistance: an assessment based on wheat leaf rust epidemics[END_REF][START_REF] Troth | Selective trade-offs maintain alleles underpinning complex trait variation in plants[END_REF][START_REF] Garnault | Spatiotemporal dynamics of fungicide resistance in the wheat pathogen Zymoseptoria Tritici in France[END_REF] and these variations may influence the evolution of various features of organisms, for example local adaptation [START_REF] Kawecki | Conceptual issues in local adaptation[END_REF][START_REF] Débarre | Quantifying the effects of migration and mutation on adaptation and demography in spatially heterogeneous environments[END_REF][START_REF] David | Adaptive diversification in heterogeneous environments[END_REF], specialisation [START_REF] Poisot | A conceptual framework for the evolution of ecological specialisation[END_REF], dispersal [START_REF] Friedenberg | Experimental evolution of dispersal in spatiotemporally variable microcosms[END_REF] or phenotypic heterogeneity [START_REF] Sasaki | The evolutionarily stable phenotype distribution in a random environment[END_REF][START_REF] Ackermann | A functional perspective on phenotypic heterogeneity in microorganisms[END_REF]. Thus, these variations are involved in numerous practical issues, for example they may generate biological diversity in particular when considering divergent selection in subdivided populations [START_REF] Dobzhansky | Nothing in biology makes sense except in the light of evolution[END_REF][START_REF] Kocher | Adaptive evolution and explosive speciation: the cichlid fish model[END_REF][START_REF] Enjalbert | Dynamic management of crop diversity: From an experimental approach to on-farm conservation[END_REF], they may allow some genotypes to better develop in some environments [START_REF] Rhoné | Insight into the genetic bases of climatic adaptation in experimentally evolving wheat populations[END_REF][START_REF] Gautier | A whole genome Bayesian scan for adaptive genetic divergence in West African cattle[END_REF][START_REF] Pariaud | Aggressiveness and its role in the adaptation of plant pathogens[END_REF][START_REF] Fabre | Epidemiological and evolutionary management of plant resistance: optimizing the deployment of cultivar mixtures in time and space in agricultural landscapes[END_REF][START_REF] Ropars | Diversity and mechanisms of genomic adaptation in penicillium[END_REF] or they may be involved in the emergence of resistances to pesticides and drugs [START_REF] Papaïx | Influence of cultivated landscape composition on variety resistance: an assessment based on wheat leaf rust epidemics[END_REF]REX Consortium, 2013;[START_REF] Garnault | Spatiotemporal dynamics of fungicide resistance in the wheat pathogen Zymoseptoria Tritici in France[END_REF]. Natural selection has a strong impact on the management and breeding of crop species for instance in crop diversity conservation. In in situ conservation approaches, defined in the International Treaty on Plant Genetic Resources for Food and Agriculture (ITPGRFA 2009) as "the conservation of ecosystems and natural habitats and the maintenance and recovery of viable populations of species in their natural surroundings and, in the case of domesticated or cultivated plant species, in the surroundings where they have developed their distinctive properties", populations of wild or domesticated species are grown over time in environments submitted to all evolutionary forces including natural selection [START_REF] Rhoné | Insight into the genetic bases of climatic adaptation in experimentally evolving wheat populations[END_REF][START_REF] Enjalbert | Dynamic management of crop diversity: From an experimental approach to on-farm conservation[END_REF]. Such populations evolve and may adapt to changing climatic conditions, emerging diseases and agricultural practices depending on the type, strength and direction of selection. Based on this principle, different strategies have been proposed focusing on the conservation and renewing of crop genetic diversity (dynamic management, [START_REF] Enjalbert | Dynamic management of crop diversity: From an experimental approach to on-farm conservation[END_REF] or on the improvement of the population (evolutionary plant breeding, [START_REF] Suneson | An evolutionary plant breeding method[END_REF]. More recently, decentralised participatory breeding approaches have developed to produce populations adapted to local contexts and practices with low input requirements in traditional or industrialised agricultural contexts [START_REF] Sperling | A framework for analyzing participatory plant breeding approaches and results[END_REF][START_REF] Ceccarelli | Farmer participation in barley breeding in Syria, Morocco and Tunisia[END_REF][START_REF] Dawson | Collaborative plant breeding for organic agricultural systems in developed countries[END_REF]. Characterisation of selective pressures on these populations is important to anticipate their evolution and optimise their management.

Selection can be inferred in various ways, in particular using retrospective methods (Barton et al., 2007, Chap. 19;Lefèvre et al., 2016, p. 665). This approach consists in measuring quantitative traits on individuals sampled in a population at several time points. Phenotypic evolution is then modelled taking account of selection and parameters are estimated by fitting the model to the data. For example such methods have been used to quantify selection for Darwin's finches [START_REF] Grant | Predicting microevolutionary responses to directional selection on heritable variation[END_REF]. Unlike prospective methods, they do not require measures of fitness components [START_REF] Lande | The measurement of selection on correlated characters[END_REF]Barton et al., 2007, p. 525).

The main retrospective method that has been used is based on the breeder's equation [START_REF] Lande | Natural selection and random genetic drift in phenotypic evolution[END_REF][START_REF] Grant | Predicting microevolutionary responses to directional selection on heritable variation[END_REF]Gillespie, 1998, Chap. 5). It assumes that genetic variance is constant, but this assumption may not be satisfied in practice [START_REF] Mcguigan | Studying phenotypic evolution using multivariate quantitative genetics[END_REF]. A few quantitative genetics models allowing genetic variance to evolve have been developed [START_REF] Cavalli-Sforza | Evolution of continuous variation: Direct approach through joint distribution of genotypes and phenotypes[END_REF][START_REF] Débarre | Multidimensional (co)evolutionary stability[END_REF][START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for evolution of selfing[END_REF]. However they are not simple because multiple effects may affect traits besides selection, in particular environmental effects, mutation or drift. As a result they are not easy to fit to data. Few statistical methods have been developed to fit such evolutionary models to phenotypic time series.

This article considers populations evolving in several locations. The objective is to develop a Bayesian retrospective method for inferring how selective pressures vary between locations and between generations. For the sake of simplicity, our study is restricted to organisms with non-overlapping generations that reproduce by selfing or asexually, for example autogamous annual plants. Individuals are characterised by a single trait that evolves under the effect of natural selection and mutation. As the model we want to use is not available in the literature, we develop our own model (Section 2). Model parameters are estimated using a Bayesian framework (Section 3). This method is assessed using simulated data (Section 4) and is applied to data from a wheat participatory plant breeding programme (Section 5).

Model

Notations and assumptions

We consider populations located in several sites and model their phenotypic evolution under the effect of mutation and selection. Selective pressures may vary between locations and between generations. Locations are indexed by i ∈ {1, . . . , I}. Populations are assumed to be nested within locations and are indexed by (i, j), j ∈ {1, . . . , J i }.

Populations are assumed to be isolated and large enough that demographic stochasticity (or drift) can be neglected.

The life cycle is divided into a reproduction phase and a growth one (Fig. 1). For the sake of clarity, reproduction is divided into a birth phase and a mutation one, although this distinction is only formal. Generations are non-overlapping, are indexed by t ∈ N and start after reproduction. Thus, an offspring of an adult of generation t belongs to generation t + 1 after the mutation phase but still belongs to generation t after the birth phase.

In the following, the symbols µ, σ 2 and τ 2 are often used to denote the expectation, variance and inverse of the variance of a random variable respectively, using subscripts reminding the random variable, dropping unnecessary subscripts and avoiding double subscripts. For example the expectation of Z ijtk is denoted by µ Zijt rather than by µ Z ijtk , if it depends on i, j and t and if it does not depend on k. We let p(X = x) denote the value of the density of the random variable X at point x. We use the common notation N (µ, σ 2 ) for the normal distribution with mean µ and variance σ 2 .

Phenotype

Adults are characterised by an adaptive trait Z ∈ R affecting fitness. This trait is assumed to be the sum of genetic and environmental effects (Gillespie, 1998, p. 105;[START_REF] Barton | The infinitesimal model: Definition, derivation, and implications[END_REF]. Consider an adult individual sampled in population (i, j) at generation t. Its trait Z ijt is assumed to be equal to

Z ijt = G ijt + A it + B ijt + C ijt .
In this equation, G ijt represents a random genetic effect. A it represents an environmental effect that is common to all the individuals of location i. B ijt is an environmental effect that is common to all the individuals of population (i, j) and that is peculiar to this population. The effects A it and B ijt can be regarded as fixed parameters, for which we will assume a prior distribution in Section 3. C ijt represents a local environmental effect that is peculiar to the individual and that is assumed to be normally distributed:

C ijt ∼ N (0, σ 2 Cijt ).
Note that the variance σ 2 Cijt is assumed to depend on the population and the generation. The effects G ijt and C ijt are assumed to be independent, thus, the expectation µ Zijt and variance σ 2

Zijt of Z ijt satisfy

µ Zijt = µ Gijt + A it + B ijt , σ 2 Zijt = σ 2 Gijt + σ 2 Cijt .
(1)

Birth

During the birth phase, adults produce offsprings and die. Individuals are assumed to reproduce asexually or by selfing. Fecundity is assumed to have a normal shape [START_REF] Cavalli-Sforza | Evolution of continuous variation: Direct approach through joint distribution of genotypes and phenotypes[END_REF][START_REF] Johnson | Theoretical models of selection and mutation on quantitative traits[END_REF][START_REF] Zhang | Genetic variability under mutation selection balance[END_REF]. Thus, an adult with phenotype z ∈ R in location i at generation t gives birth to F it (z) offsprings, where F it (z) satisfies

F it (z) ∝ exp -τ 2 F i (z -µ F it ) 2 /2 , µ F it = θ i + ε it , ε it ∼ N (0, σ 2 ε ).
This fecundity function leads to a stabilising selection toward the optimal phenotype µ F it with a strength quantified by τ F i . Optimal phenotypes fluctuate randomly around θ i , which represents the mean optimal phenotype over generations for location i [START_REF] Sasaki | The evolutionarily stable phenotype distribution in a random environment[END_REF][START_REF] Scheiner | Bet-hedging as a complex interaction among developmental instability, environmental heterogeneity, dispersal, and life-history strategy[END_REF]. Note that F it is a deterministic function of the parental phenotype given ε it and that populations in a given location are assumed to undergo the same selective pressures. Selective pressures are allowed to vary spatially, since τ F i and µ F it depend on i, and temporally, since µ F it depends on t.

G ijt

Observations t

? - Birth µ F it , τ F i G ijt Mutation -G ij(t+1) σ M - Growth A i(t+1) , B ij(t+1) C ij(t+1) G ij(t+1)
Observations t + 1

Mutation and growth

After birth, offspring genetic effects are modified by mutation. Their distribution after mutation is assumed to be normal around their value before mutation with variance σ 2 M . During growth, individuals undergo environmental effects, express their phenotype and are assumed to survive with the same probability. Thus, the survival probability does not influence the evolution of populations.

Genetic evolution

The distribution of genetic effects evolves under the effect of selection and mutation (Fig. 1). For t ∈ N, i ∈ {1, . . . , I} and j ∈ {1, . . . , J i }, consider an adult sampled in population (i, j) at generation t with genetic effect G ijt and trait Z ijt . Likewise consider an offspring sampled independently in the population at the end of generation t after birth and let G ijt denote its genetic effect. Individuals are assumed to be homozygous or haploid and to transmit their genetic effects to their offsprings but not their environmental effects. Thus, for g ijt ∈ R, the probability density of genetic effects after birth satisfies

p(G ijt = g ijt ) = R F it (z ijt ) p(Z ijt = z ijt , G ijt = g ijt ) dz ijt R R F it (z ijt ) p(Z ijt = z ijt , G ijt = g ijt ) dz ijt dg ijt , = w it (g ijt ) p(G ijt = g ijt ) R w it (g ijt ) p(G ijt = g ijt ) dg ijt , (2)
where

w it (g ijt ) = R F it (z ijt ) p(Z ijt = z ijt | G ijt = g ijt ) dz ijt and p(Z ijt = z ijt | G ijt = g ijt ) denotes the density of Z ijt given {G ijt = g ijt }.
For g ij(t+1) ∈ R, the density of genetic effects after mutation writes

p(G ij(t+1) = g ij(t+1) ) = R p(G ij(t+1) = g ij(t+1) | G ijt = g ijt ) p(G ijt = g ijt ) dg ijt ,
where G ij(t+1) is the genetic effect after mutation of the sampled offspring, with distribution N (g ijt , σ 2 M ) given {G ijt = g ijt }. The distribution of genetic effects does not change during growth.

If the initial distribution of genetic effects in the adult population is normal:

G ij0 ∼ N (µ Gij0 , σ 2
Gij0 ), then the distributions of genetic effects among adults in the following generations remain normal (Supplementary material A):

∀t ∈ N, G ij(t+1) ∼ N (µ Gij(t+1) , σ 2 Gij(t+1) ), (3)
where µ Gij(t+1) and σ 2 Gij(t+1) satisfy the following recursions:

µ Gij(t+1) = τ 2 Gijt µ Gijt + τ 2 Sijt (µ F it -A it -B ijt ) τ 2 Gijt + τ 2 Sijt , (4)
τ 2 Gij(t+1) = (τ 2 Gijt + τ 2 Sijt ) τ 2 M τ 2 Gijt + τ 2 Sijt + τ 2 M , (5)
τ 2 Sijt = τ 2 Cijt τ 2 F i τ 2 Cijt + τ 2 F i
,

where τ 2 Gijt = σ -2 Gijt , τ 2 M = σ -2
M , τ Sijt is the overall selection efficiency and τ 2 Cijt = σ -2 Cijt . Overall selection efficiency increases with selection strength τ F i and decreases with environmental variance σ 2

Cijt . Genetic variance σ 2 Gij(t+1) decreases with selection efficiency τ Sijt and increases with mutation variance σ 2 M . In the recursions above, it is assumed that τ Gij0 > 0 and τ M > 0, so that the ratios are well defined. However we allow τ F i to be equal to 0, which corresponds to the case when evolution is neutral.

Bayesian inference

Observations

Our objective is to estimate model parameters from phenotypic data. For t ∈ {0, . . . , T }, i ∈ {1, . . . , I} and j ∈ {1, . . . , J i }, K ijt > 1 adults are sampled independently in population (i, j) at generation t and their trait is measured (Fig. 1). Measurement errors are not confounded with local environmental effects since they do not affect trait values and do not interact with selection. However it is difficult to distinguish them from genetic and environmental effects without replicating measurements on some individuals. Thus, measurement errors are here neglected.

For k ∈ {1, . . . , K ijt }, the observation Z ijtk for individual k follows a normal distribution:

Z ijtk ∼ N (µ Zijt , σ 2
Zijt ), where the expressions for µ Zijt and σ 2

Zijt result from (1), ( 4) and (5).

Prior distribution

A possible joint prior distribution of model parameters is composed of the following independent distributions. Vague normal prior distributions are placed on mean optimal phenotypes and initial genetic means:

θ i ∼ N (α θ , β 2 θ ), µ Gij0 ∼ N (α G , β 2 G )
, where the parameters of these distributions are known constants.

Vague prior distributions on standard deviations are usually specified as half-Student or half-normal distributions [START_REF] Gelman | Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper)[END_REF]. Thus, vague half-normal prior distributions are placed on the standard deviations of mutation effects and of the effects affecting optimal phenotypes:

σ M ∼ N + (0, γ 2 M ), σ ε ∼ N + (0, γ 2 ε )
, where γ M and γ ε are known constants and N + (µ, σ 2 ) is the normal distribution restricted to positive values with parameters µ and σ.

Hierarchical distributions are placed on between-population environmental effects (Robert, 2007, p. 460):

A it ∼ N (0, σ 2 A ), σ A ∼ N + (0, γ 2 A ), B ijt ∼ N (0, σ 2 B ), σ B ∼ N + (0, γ 2 B )
, where γ A and γ B are known constants. Such hierarchical distributions shrink estimates (Robert, 2007, p. 472).

The prior distributions of selection strengths are chosen to allow these parameters to take low values when evolution is neutral and large values when selection is present. Assuming that selection strengths follow vague log-normal distributions does not seem satisfactory in the neutral case since ln(τ F i ) → -∞ when τ F i → 0. Thus, half-normal prior distributions are placed on selection strengths:

τ F i ∼ N + (0, γ 2 F )
, which place a substantial mass on low values. Thus, a large estimate of τ F i in an analysis will be due to the data rather than the prior distribution. However, our analyses showed that posterior distributions could be sensitive to the value of γ F used. Thus, rather than considering γ F as a known constant, it is considered as an unknown parameter that is estimated from the data [START_REF] Gelman | Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper)[END_REF]. It is written as

γ F = ω/σ Z ,
where σZ is a prior value for the overall phenotypic standard deviation, which takes account of all possible sources of variation, and ω is an unknown parameter that is assumed to follow a vague half-normal prior distribution:

ω ∼ N + (0, γ 2 ω ),
where γ ω is a known constant.

As only the phenotypic variances are observed, the standard deviations of local environmental effects may be difficult to separate from the genetic standard deviations. As these parameters are numerous, a hierarchical distribution is placed on these parameters to facilitate their estimation:

ln(σ Cijt ) ∼ N (ln(ρ i ), γ 2 C ), ρ i ∼ N + (0, κ 2 ), γ C ∼ N + (0, η 2 ),
where κ and η are known constants and ρ i is the median of σ Cijt . This prior distribution shrinks the estimates of σ Cijt toward ρ i , thereby providing more robust estimates (Robert, 2007, p. 472).

Likewise as only the phenotypic variances are observed, the initial genetic standard deviations may be difficult to estimate. Thus, these parameters are not assumed to follow vague prior distributions but are assumed to be approximately proportional to the medians of the standard deviations of local environmental effects:

σ Gij0 ∼ N + (λ ρ i , γ 2 G ),
where λ is a known constant and γ G is a known constant with a moderate value. For example when λ is fixed to one, initial genetic standard deviations and environmental ones are assumed to have similar orders of magnitude. An example of set of prior distributions is given in Supplementary material E.

Posterior distribution

In the Bayesian framework, the prior distribution is updated into a posterior distribution, that writes (Robert, 2007, p. 22)

p(φ | y) = p(φ) p(y | φ) p(y) ,
where φ denotes the vector of model parameters, y the vector of observations, p(φ) the prior distribution of φ, p(y | φ) the likelihood and p(y) the marginal distribution of y. This posterior distribution cannot be calculated explicitly for our model but it can be estimated using Markov Chain and Monte Carlo (MCMC) methods. Such methods simulate values of φ according to a Markov chain that converges in distribution to p(φ | y) (Robert, 2007, p. 302). MCMC methods are here implemented using the R package rjags [START_REF] Plummer | rjags: Bayesian Graphical Models using MCMC[END_REF][START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. Details on our implementation of MCMC methods are given in Supplementary material B and a script showing how to implement the methods in R is given in Supplementary material C.

Simulation study

Data

The ability of the developed method to infer evolutionary parameters in heterogeneous environments was assessed using simulations. Data were simulated for experimental set-ups similar to that of the wheat application (Section 5). There were T consecutive generations and 10 locations, with 10 populations per location. At each generation, 20 individuals were sampled in each population and their trait was measured. Data were simulated for two quantitative traits using the model of Section 2. For the first (respectively second) trait, parameter values were similar to those of spike weight (respectively plant height) (Section 5). The initial genetic means were equal to µ Gij0 = µ G0 -δ, j = 1, . . . , 5, µ Gij0 = µ G0 + δ, j = 6, . . . , 10, for i = 1, . . . , 10.

Four experimental set-ups, defined by the combinations of two two-level factors, were considered. The first factor quantified the amount of information provided by the set-up, so that the data were either weakly informative or informative. Set-ups providing more information had more generations, smaller between-population environmental variances and a larger δ value. Specifically, parameter values were equal to T = 4, σ A = 0.4, σ B = 0.3, γ C = 0.2, δ = 0.3 for the weakly informative set-up and trait 1, to T = 6, σ A = 0.2, σ B = 0.2, γ C = 0.1, δ = 0.6 for the informative set-up and trait 1, to T = 4, σ A = 100, σ B = 50, γ C = 0.2, δ = 100 for the weakly informative set-up and trait 2, and to T = 6, σ A = 30, σ B = 30, γ C = 0.1, δ = 200 for the informative set-up and trait 2. The second factor indicated whether selection was absent or present. When selection was present, τ F i = 1 (respectively τ F i = 0.003) for trait 1 (respectively trait 2) and selective pressures varied both spatially and temporally. In particular, mean optimal phenotypes depended on the location: θ i = µ G0 -∆, i = 1, . . . , 5, θ i = µ G0 + ∆, i = 6, . . . , 10. The parameter values used in the simulations are given in Supplementary material D. For each combination of parameter values, four data sets were simulated independently.

Analysis

The data were analysed using the model of Sections 2 and 3 with γ ω = 0.2. Mutation effects were assumed to be a priori smaller than environmental effects [START_REF] Lefèvre | The rate of polygenic mutation[END_REF][START_REF] Houle | Comparing mutational variabilities[END_REF][START_REF] Johnson | Theoretical models of selection and mutation on quantitative traits[END_REF], since γ M had a smaller value than γ A , γ B and κ (Supplementary material D). Likewise, γ ε had a smaller value than γ A , γ B and κ to help to separate temporal variations of optimal phenotypes from environmental effects. Parameter λ was fixed to 1. The values of the other parameters of prior distributions are given in Supplementary material D.

The sensitivity of posterior distributions to variations of γ ω values was investigated. To reduce computing time, the posterior means for γ ω = 0.1 and γ ω = 0.3 were estimated using importance sampling and the posterior distribution for γ ω = 0.2 as proposal distribution (Carlin and Louis, 2008, p. 182).

The main evolutionary parameters were estimated: selection strengths and mean optimal phenotypes, that characterised spatial variations of selective pressures; σ ε , that characterised temporal variations of selective pressures; mutation standard deviation. Differences between estimates and true parameter values were quantified using the average absolute error (aae). For example, for selection strengths, experimental set-up s and a given value of γ ω , aae was equal to aae(s, γ ω ) = 1 4

4 r=1 1 10 10 i=1 | τ F i (s, r, γ ω ) -τ F i (s)| ,
where s = 1, . . . , 4, τ F i (s) was selection strength for location i and set-up s and τ F i (s, r, γ ω ) was the posterior mean of this parameter for set-up s, replication r and parameter value γ ω .

Results

The estimates of selection strengths and mutation standard deviation showed lower errors and were less sensitive to variations of prior distributions for informative set-ups than for weakly informative set-ups (Figs. 2 and 3; Figs. D.1 and D.3 in Supplementary material D). For weakly informative set-ups, these estimates both increased with γ ω , presumably because an increase in selection strengths, that reduced within-population diversities, was counterbalanced by an increase in mutation standard deviation, that increased within-population diversities, in order to give within-population diversities consistent with the data.

Estimation errors were rather large for mean optimal phenotypes (Figs. 2 and 3;Figs. D.2 and D.4 in Supplementary material D) and for σ ε and trait 2 (Fig. 3). The 95% Highest Posterior Density (HPD) credible intervals of σ ε were rather wide, for example this interval was equal to [0.000,0.304] (respectively [0.143,99.412]) for trait 1 (respectively trait 2), the informative set-up with selection, replicate 1 and γ ω = 0.2.

Application to wheat data

Data

A participatory plant breeding programme on bread wheat was set up in France by farmers and facilitators from the Réseau Semences Paysannes1 and researchers from INRA (GQE -le Moulon) to develop heterogeneous populations adapted to farmers' local conditions, practices and objectives [START_REF] Dawson | Collaborative plant breeding for organic agricultural systems in developed countries[END_REF]. Some of these populations were exchanged among farmers. They then evolved several years in the same farms in Table 1: Wheat data. The number of years corresponds to the mean number of consecutive years a population is sown in a farm, while a trial is defined as a farm x year combination.

field trials and sometimes underwent farmers' mass selection. Whenever populations were cultivated on-farm, some characters such as plant height, spike length and spike weight were measured on several individuals per population (10 to 50).

The data used for this analysis consisted of a subset of the data gathered since the beginning of the project in 2006 (Tab. 1). Two variables were analysed: spike weight and plant height. The populations selected (i) evolved for at least three years in a given farm, (ii) were measured for at least three consecutive years, (iii) had at least 10 observations per year and (iv) were not subjected to mass selection. We studied the selective pressures undergone by these populations grown in different environments and the influence of natural selection on their evolution. In this application, locations were farms, trials were farm x year combinations and on average populations were sown and measured for 3.5 and 4.3 consecutive years for spike weight and plant height respectively. Measurement errors, that were neglected in our model, were likely to be smaller for spike weight than for plant height, since the former was measured in the laboratory while the latter was measured in the field. More information on the data collected in each farm is given in Supplementary material E.

Analysis

As wheat is an autogamous annual plant, the data were analysed using the model of Section 2 and the prior distributions given in Supplementary material E. As in the simulation study, γ M and γ ε had a smaller value than γ A , γ B and κ, and λ was fixed to 1. Parameter γ ω took two possible values, γ ω = 0.1 or γ ω = 0.4, in order to study the influence of this parameter on posterior distributions.

Results

Selection strengths and mean optimal phenotypes. The posterior distributions of τ F i and θ i were sensitive to changes in γ ω values for both variables (Figs. 4 and5). For spike weight, the HPD credible intervals of selection strengths contained values close to 0 for all the locations when γ ω = 0.1 (Fig. 4,a), so that the presence of selection was not significant. On the contrary when γ ω = 0.4, the presence of selection was significant in 4 locations (Fig. 4,b). The significant estimate of selection strength in location 5 could be a false positive since there were only 6 populations observed for 3 years in this location (Supplementary material E). The posterior medians of mean optimal phenotypes were 0 1 2 3 a q q q q q q q q q q q q τ F12 τ

F11 τ F10 τ F9 τ F8 τ F7 τ F6 τ F5 τ F4 τ F3 τ F2 τ F1 0 1 2 3 b q q q q q q q q q q q q τ F12 τ F11 τ F10 τ F9 τ F8 τ F7 τ F6 τ F5 τ F4 τ F3 τ F2 τ F1
0.5 1.0 1.5 2.0 2.5 3.0 c q q q q q q q q q q q q θ 12 θ 11 θ 10 θ 9 θ 8 θ 7 θ 6 θ 5 θ 4 θ 3 θ 2 θ 1 0.5 1.0 1.5 2.0 2.5 3.0 d q q q q q q q q q q q q θ 12 θ 11 θ 10

θ 9 θ 8 θ 7 θ 6 θ 5 θ 4 θ 3 θ 2 θ 1
Figure 4: Analysis of wheat data: selection strengths and mean optimal phenotypes for spike weight. Estimates (•) and 95% HPD credible intervals (solid lines) of selection strengths in 12 locations for γ ω = 0.1 (a) and γ ω = 0.4 (b). Sample location means (×), estimates (•) and 95% HPD credible intervals (solid lines) of mean optimal phenotypes in the same locations for γ ω = 0.1 (c) and γ ω = 0.4 (d). Results are in g for mean optimal phenotypes and in g -1 for selection strengths.

close to 2 when γ ω = 0.1, while they ranged from 1.320 to 2.199 and tended to be smaller than sample means when γ ω = 0.4 (Fig. 4, c andd). For plant height, the presence of selection was not significant for both values of γ ω , but estimates of selection strengths and mean optimal phenotypes varied with this parameter (Fig. 5).

Environmental and mutation standard deviations. When γ ω = 0.1, the credible intervals of σ A , σ B and ρ i did not contain values close to 0 (Fig. 6). The estimates of σ A were larger than those of σ B . The estimates of ρ i appeared more variable between locations for spike weight than for plant height. The posterior median of σ M was equal to 0.033 (respectively 0.132) when γ ω = 0.1 (respectively γ ω = 0.4) for spike weight. It was equal to 21.316 (respectively 22.705) when γ ω = 0.1 (respectively γ ω = 0.4) for plant height. The credible intervals of σ ε were rather wide for both values of γ ω . The 95% 0.000 0.002 0.004 0.006 0.008

a q q q q q q q τ F7 τ F6 τ F5 τ F4 τ F3 τ F2 τ F1
0.000 0.002 0.004 0.006 0.008

b q q q q q q q τ F7 τ F6 τ F5 τ F4 τ F3 τ F2 τ F1
600 800 1000 1200 1400 1600 c q q q q q q q θ 7 θ 6 θ 5 θ 4 θ 3 θ 2 θ 1 600 800 1000 1200 1400 1600 d q q q q q q q θ 7 θ 6 Results are in mm for mean optimal phenotypes and in mm -1 for selection strengths.

0.0 0.2 0.4 0.6 a q q q q q q q q q q q q q q q q ρ 12 ρ 11

ρ 10 ρ 9 ρ 8 ρ 7 ρ 6 ρ 5 ρ 4 ρ 3 ρ 2 ρ 1 σ B σ A σ M
σ ε 0 50 100 150 200 b q q q q q q q q q q q ρ 7 ρ 6 credible interval of γ C was equal to [0.341,0.434] (respectively [0.334,0.441]) for spike weight (respectively plant height) and γ ω = 0.1, and thus did not contain values close to 0. Thus, the standard deviations of local environmental effects varied between years and populations, in agreement with observed phenotypic standard deviations (Fig. 7, b andd).

ρ 5 ρ 4 ρ 3 ρ 2 ρ 1 σ B σ A σ M σ ε
Population means and standard deviations. Estimates of population phenotypic means and standard deviations were consistent with observed values (Fig. 7). Genetic means and standard deviations tended to evolve toward equilibrium values [START_REF] Cavalli-Sforza | Evolution of continuous variation: Direct approach through joint distribution of genotypes and phenotypes[END_REF][START_REF] Johnson | Theoretical models of selection and mutation on quantitative traits[END_REF][START_REF] Zhang | Genetic variability under mutation selection balance[END_REF].

Discussion

Selection. This article develops a quantitative genetics model involving selection. Here fecundity is assumed to have a normal shape so that selection favours an interme- q q q q q q 0 1 2 3 4 700 800 900 1100 c Year Mean q q q q q 0 1 2 3 4 50 60 70 80 90 100 d Year Standard deviation q q q q q Figure 7: Analysis of wheat data: population means and standard deviations. Sample phenotypic means ( ), estimated phenotypic means ( ) and estimated genetic means (•) for a population of location 2 and spike weight (a). Estimated phenotypic (respectively genetic) means are the posterior medians of µ Zijt (respectively µ Gijt ). Sample phenotypic standard deviations ( ), estimated phenotypic standard deviations ( ) and estimated genetic standard deviations (•) for the same population and spike weight (b). Estimated phenotypic (respectively genetic) standard deviations are the posterior medians of σ Zijt (respectively σ Gijt ). Sample phenotypic means ( ), estimated phenotypic means ( ) and estimated genetic means (•) for a population of location 5 and plant height (c). Sample phenotypic standard deviations ( ), estimated phenotypic standard deviations ( ) and estimated genetic standard deviations (•) for the same population and plant height (d). In these analyses, γ ω = 0.1. Results are in g (respectively mm) for spike weight (respectively plant height). diate trait value. This type of selection appears common in nature [START_REF] Kingsolver | The strength of phenotypic selection in natural populations[END_REF] and leads to explicit formulae for the distribution of observations, which simplifies the implementation of statistical methods. Fecundity may also be modelled using a quadratic function [START_REF] Débarre | Multidimensional (co)evolutionary stability[END_REF]. This model takes account not only of stabilising selection but also of directional selection, which favours extreme trait values in one direction, and of diversifying selection, which favours extreme trait values. However it assumes that selection is weak and that genetic effects have a normal distribution at each generation. On the contrary, the model with a normal fecundity function does not assume that selection is weak and assumes that genetic effects have a normal distribution at the first generation only. This is why we here restricted ourselves to a normal fecundity function.

Selection on a trait may be real or apparent [START_REF] Johnson | Theoretical models of selection and mutation on quantitative traits[END_REF][START_REF] Zhang | Genetic variability under mutation selection balance[END_REF]Barton et al., 2007, p. 526). Selection is real when fitness differences are caused by the trait of interest and apparent when the trait of interest is neutral but evolves under the effect of the selection acting on correlated traits. Our model is univariate and is not able to distinguish between real or apparent selections. The targets of selection could be better studied by measuring many traits and by developing multivariate models allowing for the correlations between traits [START_REF] Lande | The measurement of selection on correlated characters[END_REF][START_REF] Troth | Selective trade-offs maintain alleles underpinning complex trait variation in plants[END_REF]. However this approach can hardly measure all the traits that are under real selection.

In principle, selection may affect both fecundity and survival. Here it is assumed to affect fecundity only to adapt the model to the available data. As our data are collected on adults, they provide more information on the selective pressures occurring during reproduction than during growth. In addition it is difficult to separate the selective pressures occurring during these phases using a single observation time per generation.

In the simulation study, estimates tend to be less sensitive to prior distributions when the data are more informative. Posterior distributions combine a prior information with the information brought by the data. Thus, they may be sensitive to changes in prior distributions when the data provide little information. This is the case in the wheat application, so that it is difficult to infer how selection varies between farms and between years. This is presumably because phenotypic time series are short and perturbed by large environmental effects. It is recommended to carry out a sensitivity analysis to prior distributions. If the results are sensitive to prior distributions, one may report this sensitivity and conclude that the available data do not allow to estimate parameters precisely. Alternatively, additional information may be introduced in the analysis by collecting more data (Carlin and Louis, 2008, p. 181).

For the sake of simplicity, the significance of estimated selection strengths are here assessed by examining if HPD credible intervals contain low values. Testing if these parameters vanish formally is difficult because zero is the boundary of the parameter space [START_REF] Müller | Model selection in linear mixed models[END_REF][START_REF] Baey | Asymptotic distribution of likelihood ratio test statistics for variance components in nonlinear mixed effects models[END_REF]. More work is needed to provide rigorous tests of selection strengths.

Environmental effects. In our model, environmental effects are decomposed into between-and within-population effects, with standard deviations σ A , σ B and σ Cijt . Both effects appear present in the wheat application. The standard deviation of local effects σ Cijt is allowed to vary between generations to take account of environmental stochasticity (Caswell, 2001, Chap. 14). In the wheat application, this stochasticity could be due to variations in climatic conditions, soil fertility or farming practices over years.

Model assumptions. Our model assumes that populations are isolated and that individuals reproduce by selfing or asexually and are homozygous or haploid. These assumptions are partly justified in our application since wheat plants mainly reproduce by selfing, but not fully however since new genotypes may appear in populations because plants are not fully homozygous, because plants reproduce by a mixture of self fertilisation and outcrossing and because of seed dispersal. Thus, our estimates of mutation effects could capture various diversifying effects besides mutation. The selfing assumption could be relaxed by using the model developed by [START_REF] Cavalli-Sforza | Evolution of continuous variation: Direct approach through joint distribution of genotypes and phenotypes[END_REF], which is based on outcrossing.

This article considers large populations and neglects demographic stochasticity, thus the developed method should be applied to small populations with caution. Populations are not very large in the wheat application, so that demographic stochasticity should be present. However its influence on the estimation of selection parameters should be limited since locations hold several populations. Evaluating the relative importances of evolutionary forces remains an open issue in evolutionary biology [START_REF] Wright | Evolution in mendelian populations[END_REF][START_REF] Lande | Natural selection and random genetic drift in phenotypic evolution[END_REF]Barton et al., 2007, Chap. 19). Our model allows to investigate the relative importances of selection and mutation but not those of selection and demographic stochasticity.

Figure 2 :Figure 3 :

 23 Figure 2: Analysis of simulated data: sensitivity analysis for trait 1, selection strengths (a), mean optimal phenotypes (b), σ ε (c) and σ M (d). The quality of estimates is quantified by average absolute error (aae). : Weakly informative set-up without selection; •: Informative set-up without selection; : Weakly informative set-up with selection; : Informative set-up with selection. Estimation errors for mean optimal phenotypes and σ ε are relevant in the presence of selection only.

Figure 5 :

 5 Figure5: Analysis of wheat data: selection strengths and mean optimal phenotypes for plant height. Estimates (•) and 95% HPD credible intervals (solid lines) of selection strengths in 7 locations for γ ω = 0.1 (a) and γ ω = 0.4 (b). Sample location means (×), estimates (•) and 95% HPD credible intervals (solid lines) of mean optimal phenotypes in the same locations for γ ω = 0.1 (c) and γ ω = 0.4 (d). Results are in mm for mean optimal phenotypes and in mm -1 for selection strengths.

Figure 6 :

 6 Figure 6: Analysis of wheat data: environmental and mutation standard deviations. Estimates (•) and 95% HPD credible intervals (solid lines) for spike weight (in g, a) and plant height (in mm, b). In these analyses, γ ω = 0.1.

?

  Figure 1: Life cycle diagram, observation times and genetic composition for population j of location i. The life cycle is decomposed into birth, mutation and growth phases. The population is observed when individuals are adults. Individuals are characterised by an adaptive trait that has a genetic component. The random variable G ijt is the genetic effect of an adult sampled in the population at generation t. The random variables G ijt and G ij(t+1) are the genetic effects before and after mutation respectively of an offspring sampled in the population at the end of generation t. Selection favours individuals with trait values close to µ F it with a strength of τ F i . Mutation effects with standard deviation σ

M increase genetic variance. Environmental effects A i(t+1) , B ij(t+1) and C ij(t+1) occur during growth.

www.semencespaysannes.org

Acknowledgements

This work was supported by INRA (project AgriBio 4). We are grateful to the INRA MIGALE bioinformatics platform (http://migale.jouy.inra.fr) for providing computational resources. We thank Charles Ladmiral and two anonymous reviewers for helpful comments.