Denis Roze
email: roze@sb-roscoff.fr

Alexandre Blanckaert

AND THE MUTATION LOAD IN SEXUAL AND ASEXUAL POPULATIONS

Keywords: EPISTASIS, PLEIOTROPY, adaptive landscape, epistasis, evolutionary quantitative genetics, multilocus models, mutation load, stabilizing selection

whether they are published or not. The documents may come

INTRODUCTION

Although mutation represents the ultimate fuel for adaptation, it is also the source of a fitness cost for populations due to the production of sub-optimal genotypes. This "mutation load" may in turn affect many important evolutionary processes, such as the evolution of sex and recombination [START_REF] Kondrashov | Deleterious mutations and the evolution of sexual reproduction[END_REF]Barton, 1995;[START_REF] Keightley | Interference among deleterious mutations favours sex and recombination in finite populations[END_REF][START_REF] Otto | The evolutionary enigma of sex[END_REF], inbreeding depression and the evolution of mating systems [START_REF] Lande | The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models[END_REF][START_REF] Charlesworth | The genetic basis of inbreeding depression[END_REF][START_REF] Charlesworth | Evolution of plant breeding systems[END_REF], mate choice [START_REF] Rowe | The lek paradox and the capture of genetic variance by condition dependent traits[END_REF] or ploidy levels [START_REF] Otto | Recombination and the evolution of diploidy[END_REF][START_REF] Otto | Mating systems and the evolutionary transition between haploidy and diploidy[END_REF]. Different types of models have explored the effect of recurrent mutation on the average fitness of populations. The simplest model assumes that each mutation decreases fitness by a fixed factor, independently of the genetic background (multiplicative model). In that case, the average fitness of a population at mutation-selection balance (relative to the maximal possible fitness) is approximately e -U , where U is the genomic rate of deleterious mutation (e.g., [START_REF] Crow | Genetic loads and the cost of natural selection[END_REF]. This results holds for both sexual (randomly mating) and asexual populations, as long as stochastic effects can be neglected. However, it seems unlikely that most mutations have independent effects: for example, direct measures of fitness effects of mutations (alone and in combination) in microorganisms usually show a wide distribution of epistatic interactions among pairs of mutations (e.g., [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF]. [START_REF] Kimura | The mutational load with epistatic gene interactions in fitness[END_REF] explored the effects of epistasis among deleterious alleles on mean fitness, assuming that epistasis is the same for all pairs of mutations.

They showed that when interactions among mutations tend to reinforce their deleterious effect (negative epistasis) the mean fitness of sexual populations increases, while it decreases when deleterious alleles tend to compensate each other (positive epistasis).

By contrast, the mean fitness of asexual populations is not affected by epistasis and remains approximately e -U .

In a different class of models, mutations affect a given number of phenotypic traits which in turn influence fitness. [START_REF] Fisher | The genetical theory of natural selection[END_REF], [START_REF] Haldane | The causes of evolution[END_REF] and [START_REF] Wright | Evolution in populations in approximate equilibrium[END_REF] considered different models of stabilizing selection acting on a single quantitative trait influenced by many genes. In particular, Haldane and Wright derived expressions for equilibrium allele frequencies, from which mean fitness can be deduced; further developments (to which we will return below) were done in the sixties and seventies, in particular by [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF], [START_REF] Bulmer | The effect of selection on genetic variability[END_REF]Bulmer ( , 1972) ) and [START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF]Lande ( , 1977)). Interestingly, compensatory effects between mutations emerge naturally from this type of model (a mutation that displaces from the optimum can be compensated by another mutation bringing closer to the optimum), together with distributions of epistatic interactions among mutations. In his Genetical Theory of Natural Selection, Fisher also proposed a model involving multiple phenotypic traits in support of his idea that adaptation is mainly due to mutations of small effect: this geometrical model represents an n-dimensional phenotypic space where each dimension corresponds to the value of a quantitative trait, and where mutations correspond to random vectors displacing individuals in phenotypic space. Fitness is assumed to decrease monotonously as the Euclidean distance from a point corresponding to the optimal phenotype increases.

Although the initial goal of Fisher's geometrical model was to describe the dynamics of adaptation [START_REF] Fisher | The genetical theory of natural selection[END_REF][START_REF] Orr | The population genetics of adaptation: the distribution of factors fixed during adaptive evolution[END_REF]Orr, , 2000;;[START_REF] Welch | Modularity and the cost of complexity[END_REF], it has also been used to explore the effects of various parameters (population size, number of phenotypic traits under selection, shape of the fitness peak...) on the mean fitness of populations in equilibrium situations. Several authors considered a weak-mutation limit (N U 1, where N is population size and U the total rate of mutations affecting the traits under selection) so that at most one mutation segregates in the population at a given time [START_REF] Hartl | Towards a theory of evolutionary adaptation[END_REF][START_REF] Poon | Compensating for our load of mutations: freezing the meltdown of small populations[END_REF][START_REF] Sella | The application of statistical physics to evolutionary biology[END_REF][START_REF] Tenaillon | Quantifying organismal complexity using a population genetic approach[END_REF][START_REF] Sella | An exact steady-state solution of Fisher's geometric model and other models[END_REF]: the population evolves away from the optimum by fixing weakly deleterious alleles by random drift, and occasionally returns near the optimum by fixing a compensatory mutation. For example, using the fitness function W = exp -d Q , where d is the Euclidean distance from the optimum in phenotypic space and Q a parameter affecting the shape of the fitness peak, [START_REF] Tenaillon | Quantifying organismal complexity using a population genetic approach[END_REF] showed that the equilibrium mean fitness of a haploid population is

approximately (1 -1/ (2N )) n Q
, where n is the number of traits affecting fitness (often called "complexity"). Note that this result does not depend on the reproductive mode of organisms (sexual or asexual) since at most one locus is polymorphic in this low mutation limit.

It seems likely that in many organisms, however, N U 1, in which case many of the loci affecting the traits under selection may be polymorphic at the same time.

In this regime, and assuming that the position of the optimum remains constant over time, the population can be represented by a collection of points centered around the optimum, while mean fitness can usually be expressed in terms of the variance of the different trait values in the population. Since the original works of [START_REF] Haldane | The causes of evolution[END_REF] and [START_REF] Wright | Evolution in populations in approximate equilibrium[END_REF], numerous models have been used to calculate the genetic variance of a quantitative trait (or a set of traits) maintained at mutation-selection balance, and different regimes have been described (e.g., Bulmer, 1989;[START_REF] Bürger | Mutation-selection balance and the evolutionary advantage of sex and recombination[END_REF][START_REF] Johnson | Theoretical models of selection and mutation on quantitative traits[END_REF]Zhang and Hill, 2005). When mutations tend to have strong fitness effects, the distribution of phenotypic values is often highly leptokurtic (with a sharp peak at the optimum) and may present a singularity [START_REF] Waxman | Pleiotropy and the preservation of perfection[END_REF]. In the limit of very strong selection, the expected contribution to future generations of individuals located away from the optimum becomes negligible; in this case, one expects that mean fitness, relative to the fitness of individuals at the optimum (W /W 0 )

should be approximately e -U when the number of new mutations per individual is

Poisson-distributed, independently of the exact shape of the fitness peak. Indeed, a recursion for the frequency of individuals at the optimum (x 0 ) is given by x 0 = e -U x 0 W 0 /W (e -U being the probability that no mutation occurs), which yields the desired result when x 0 = x 0 (e.g., [START_REF] Kimura | The mutational load with epistatic gene interactions in fitness[END_REF].

When mutations tend to have weak fitness effects, different results can be obtained by assuming that equilibrium distributions of phenotypic values are Gaussian.

Most models exploring this regime considered a single trait under selection and a Gaussian (or quadratic) fitness function, although multivariate models have also been proposed (Lande, 1980;Turelli, 1985). Different assumptions regarding the genetic architecture of traits (in particular the number of alleles per locus) have been explored.

Continuum-of-alleles models assume that an infinite number of alleles (having different phenotypic effects) segregate at each locus; expressions for the genetic variance and mean fitness under this scenario have been obtained by assuming that the equilibrium distribution of allelic effects in the population is Gaussian at each locus (e.g., [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF][START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF]. Lande (1980) generalized these models to the case where selection acts on multiple traits, assuming a multivariate Gaussian fitness function. Although the general expression for equilibrium genetic variances (and mean fitness) is cumbersome, it takes a simpler form when selection acts independently on each trait, and in the absence of correlation of mutational effects across traits (p.294 in [START_REF] Bürger | Mutation-selection balance and the evolutionary advantage of sex and recombination[END_REF].

Neglecting linkage disequilibria among loci, mean fitness can then be written as a simple function that depends on the number of selected traits, number of loci, strength of selection and mutational variance (see Results section below). Other models considered stabilizing selection (Gaussian fitness function) acting on a single trait coded by multiple biallelic loci (e.g., Bulmer, 1972Bulmer, , 1985;;[START_REF] Barton | The maintenance of polygenic variation through a balance between mutation and stabilizing selection[END_REF]Barton, , 1989)). As shown by [START_REF] Barton | The maintenance of polygenic variation through a balance between mutation and stabilizing selection[END_REF], multiple possible equilibria exist for allele frequencies at the different loci. Although the population mean phenotype may deviate from the optimum at some of these equilibria, Barton (1989) argues that perturbations generated for example by random drift should keep the population mean near the optimum. Neglecting linkage disequilibria, and as long as the effects of mutation and drift at each locus are weaker than selection, one obtains from the expression of the equilibrium genetic variance that mean fitness should be close to e -U . More recently, [START_REF] Zhang | Multivariate stabilizing selection and pleiotropy in the maintenance of quantitative genetic variation[END_REF] investigated the maintenance of genetic variation in multiple-trait models (including correlations among traits), assuming that the fitness function is multivariate Gaussian, and that selection at each locus is sufficiently strong so that mutations remain rare in the population (only one mutant allele segregates at each locus). Results for the case of a changing optimum have been derived by [START_REF] Zhang | Fisher's geometrical model of fitness landscape and variance in fitness within a changing environment[END_REF]. Finally, effects of linkage disequilibria have been explored in both continuum-of-alleles and biallelic models, for the case of a single trait under Gaussian stabilizing selection (Bulmer, 1974;[START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF][START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF], and were shown to be minor as long as the total rate of mutation on loci affecting the trait remains small, and linkage is not too tight.

In this article, we derive several simple approximations for the mutation load in both sexual and asexual populations at high mutation rate, using a general model that allows us to play with the shape of the fitness peak (which in turn affects the average sign and magnitude of epistasis among mutations), the number of traits under selection and the degree of pleiotropy of mutations. These results are tested by individualbased simulations representing stabilizing selection on traits coded by large numbers of loci. As we will see, two different regimes are observed in asexual populations: when mutations tend to have large fitness effects, mean fitness only depends on the mutation rate and is close to e -U , while under weaker mutational effects mean fitness increases as mutational effects decreases, and can be expressed in terms of the mutation rate, the dimensionality of the landscape, the shape of the fitness peak and the average deleterious effect of mutations. By contrast, three regimes are observed in the case of sexual populations: under strong mutational effects, mean fitness is again close to e -U , while it switches to a different value that only depends on the mutation rate and the shape of the fitness peak as mutational effects decrease. Interestingly, mean fitness is the same under these two regimes when the fitness function is Gaussian (W ≈ e -U ), while W > e -U in the second regime when the fitness peak is flatter than a Gaussian around the optimum, and W < e -U when the fitness peak is sharper.

Finally, as mutational effects continue decreasing a third regime is entered, where allele frequency dynamics are dominated by mutation and drift. We also obtain expressions for linkage disequilibria and find that they can have a substantial effect on fitness, in particular when the mutation rate is high and the number of traits under selection is low. given by the vector z = (z 0 , z 1 , . . . , z n ). Following others [START_REF] Wilke | Interaction between directional epistasis and average mutational effects[END_REF][START_REF] Tenaillon | Quantifying organismal complexity using a population genetic approach[END_REF][START_REF] Gros | The evolution of epistasis and its links with genetic robustness, complexity and drift in a phenotypic model of adaptation[END_REF], we assume for simplicity that the fitness function is spherically symmetric, and is given by:

MODEL

W = exp - d Q 2V S (1)
where d = n i=1 z 2 i is the Euclidean distance in phenotypic space from the optimal phenotype (located at z i = 0 for each trait i), V S measures the strength of selection, while the parameter Q determines the shape of the fitness peak: Q = 2 corresponds to a (multivariate) Gaussian fitness function, while Q < 2 (resp. Q > 2) leads to a sharper (resp. flatter) fitness peak. Traits are coded by L biallelic loci, each locus affecting a subset m ≤ n of phenotypic traits (sampled randomly and independently for each locus among the set of n traits). As in [START_REF] Lourenço | Complexity, pleiotropy and the fitness effect of mutations[END_REF], the parameter m thus measures the number of traits affected by a single mutation ("mutation pleiotropy"),

while n is the total number of traits under selection ("complexity"): m = n corresponds to maximal pleiotropy, where each mutation affects all selected traits (as in Fisher's model). The genomic mutation rate is denoted U , so that each locus mutates at rate u = U/L. For simplicity, we assume additive effects of the different loci on phenotypic traits, and no environmental variance. Denoting 0 and 1 the two alleles at each locus, and X j an indicator variable that equals 0 if an individual carries allele 0 at locus j, and 1 otherwise, the value of trait i in a given individual is given by:

z i = L j=1 r ij X j (2)
where r ij measures the effect of allele 1 at locus j on trait i: an individual carrying alleles 0 at all loci is thus at the phenotypic optimum. For each locus j, the effect of allele 1 on each of the m traits affected by the locus is sampled in a Gaussian distribution with mean 0 and variance a 2 (therefore, r ij can be positive or negative).

Finally, we denote s the mean deleterious effect of mutations (measured on log-fitness)

in a population at the optimum, given by [START_REF] Gros | The evolution of epistasis and its links with genetic robustness, complexity and drift in a phenotypic model of adaptation[END_REF]:

s = (2a 2 ) Q 2 2V S Γ m+Q 2 Γ m 2 (3)
(where Γ is Euler's gamma function), which simplifies to s = m 2 (a 2 /V S ) in the case of a Gaussian fitness function (Q = 2, in agreement with [START_REF] Martin | A general multivariate extension of Fisher's geometrical model and the distribution of mutation fitness effects across species[END_REF].

As shown by [START_REF] Gros | The evolution of epistasis and its links with genetic robustness, complexity and drift in a phenotypic model of adaptation[END_REF], the average epistasis among pairs of mutations (at the optimum, and measured on log-fitness) is given by e = 2 -2 Q/2 s, which equals zero when Q = 2 (see also ref. [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF] but it positive for Q < 2, and negative for Q > 2.

Finally, we can note that the strength of selection V S (equation 1) can be considered as a scaling parameter: using the scaled variables

z i,s = z i / (2V S ) 1 Q to measure
phenotypic traits, fitness in terms of z i,s variables becomes independent of V S , while the variance of the distribution of mutational effects on each z i,s equals a 2 s = a 2 / (2V S )

2 Q :
therefore, a 2 and V S affect the results only through the compound parameter a 2 s . In the simulations (described next), we used V S = 1/2 so that a 2 s = a 2 . As shown by equation 3, large (resp. small) values of a 2 s imply that mutations have strong (resp.

weak) fitness effects.

Simulations. Analytical predictions are tested using individual-based simulations.

Our occur at the same rate). In order to gain execution speed, the program is parallelized (using the MPI library) to run on several processors, each processor dealing with a given segment of the genome, for all individuals; execution speed can be considerably increased when the number of processors is sufficiently large so that the probability that an event (mutation or cross-over) occurs in a given segment per generation is low (in which case most genome segments stay unchanged from one generation to the next).

At the start of the simulation each genome contains only "0" alleles, which cor-

responds

RESULTS

Asexual population. Figure 1 shows the equilibrium mean fitness of an asexual population as a function of the variance of mutational effects (on scaled traits) a 2 s (left) or of the mean deleterious effect of mutations s (right), and for different values of the shape parameter Q. As a 2 s becomes large, mean fitness converges to e -U (which is close to 0.6 for U = 0.5, dashed line). For smaller values of the scaled mutational variance, however, the population reaches a different regime where W > e -U . In this case, an approximation for mean fitness can be obtained by assuming that distributions of phenotypic traits z i are Gaussian at equilibrium. As explained in Supplementary File S1, Euclidean distances from the optimum (d) then follow a χ-distribution, which yields:

ln W = - (2V G ) Q 2 2V S Γ n+Q 2 Γ n 2 (4)
where V G is the genetic variance at equilibrium (the variance of z i , which by symmetry should be the same for all traits i); for Q = 2 (Gaussian fitness function), this

simplifies to ln W = -(n/2) (V G /V S
). An expression for V G at equilibrium is derived in Supplementary File S1. Assuming that ln W ≈ ln W , one obtains that when m = n (that is, when each mutation affects all phenotypic traits):

W ≈ exp -s n U s Q Q 2+Q (5)
When Q = 2, equation 5 simplifies to exp -n 2 U s , in agreement with the result obtained from Lande's analysis (1980) when selection acts independently on the different traits, and in the absence of mutational covariances among traits (that is, when the variance-covariance matrices representing the effects of selection and mutation on the traits are diagonal). Indeed, in an asexual population the whole genome can be considered as a single locus with many alleles (provided that the number of loci is sufficiently large), in which case Lande's analysis of the continuum-of-alleles model yields V G = √ U a 2 V S when mutation and selection covariance matrices are diagonal (the equilibrium genetic variance for each trait is not affected by the other traits under selection, see also p. 294 in ref. [START_REF] Bürger | Mutation-selection balance and the evolutionary advantage of sex and recombination[END_REF]. When m < n, equation 5 still holds when the fitness function is Gaussian (Q = 2), while W is given by a slightly more complicated expression when Q = 2 (Supplementary File S1). However, this expression only depends weakly on m (for fixed s) and remains close to equation 5 in most cases (unless Q is high and m is small). As shown by Figure 1, equation 5 provides accurate predictions for W as long as a 2 s is not too high. Supplementary Figure S1 shows the genetic variance V G at equilibrium, for the same parameter values as in Figure 1.

Sexual population. Figure 2 shows the equilibrium mean fitness of a sexual population as a function of a 2 s and s. When mutations have strong effects on fitness, one obtains again that W ≈ e -U (right-most points on Figure 2). As a 2 s decreases, a second regime is entered where W reaches a new value which is independent of a 2 s , but depends on the shape parameter Q (roughly for -4 < log 10 (a 2 s ) < -2). Finally, as a 2 s continues decreasing a third regime is entered where W increases up to W ≈ 1 as a 2 s decreases. In Supplementary File S2, we show that approximations for mean fitness under these last two regimes (a 2 s small) can be obtained by assuming that distributions of phenotypic values in the population are Gaussian. As in previous models (e.g., Bulmer, 1972;[START_REF] Barton | The maintenance of polygenic variation through a balance between mutation and stabilizing selection[END_REF] one obtains that the equilibrium frequency of allele 1 at locus j (p j ) is either (i) at p j = 1/2 (when selection at locus j is weak relative to mutation) or (ii) at one of the two symmetric equilibria where p j (1 -p j ) = X < 1/4 (where X is a function of the model's parameters), the locus being closer to fixation (for allele 0 or 1) as selection increases. Assuming that all loci are at equilibrium (ii) and neglecting linkage disequilibria among loci, one obtains the following approximation for mean fitness:

W ≈ exp - 2U Q (6)
Figure 2 shows that this expression does indeed provide a correct prediction of W for intermediate values of a 2 s (dashed lines). A better approximation can be obtained by taking into account effects of the linkage disequilibria, which is done in Supplementary File S2 using the methods developed by [START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF]. As described before [START_REF] Bulmer | The effect of selection on genetic variability[END_REF](Bulmer, , 1974;;[START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF][START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF], stabilizing selection tends to generate negative covariances among loci, which affect W through two effects: a reduction of the genetic variance V G which directly increases W , and an effect on equilibrium allele frequencies and thus on the genic variance V g (see Appendix A)

which has the opposite effect on W . As long as linkage disequilibria remain small the second effect predominates, and one obtains the following approximation for W (assuming that the number of loci L is large):

W ≈ exp - 2U Q 1 + 2U n 1 r H -1 (7)
where r H is the harmonic mean recombination rate between all pairs of loci (derived in Appendix B under our simulated genetic architecture). Figure 2 shows that equation 7 does indeed provide a slightly better approximation than equation 6 for intermediate values of a 2 s (horizontal solid lines). Interestingly, Supplementary Figure S1 shows that Q has only a weak effect on the equilibrium genetic variance V G in this intermediate regime, which remains close to the value obtained for Q = 2, neglecting linkage disequilibria: V G = 2U V S /n (which is also obtained from house-of-cards models assuming a Gaussian fitness function, e.g. Turelli, 1985;[START_REF] Bürger | Mutation-selection balance and the evolutionary advantage of sex and recombination[END_REF]. Finally, when a 2 s is very small (so that mutations have very weak fitness effects) one predicts that most loci should be at the equilibrium where

p j = 1/2, in which case W ≈ exp -s (L/4) Q 2
where L is the total number of loci. However, in this regime we expect that the relative effect of genetic drift on allele frequency dynamics could be important. In Supplementary File S4, we show that when selection at each locus becomes weaker than mutation and drift, an approximation for mean fitness is given by:

W ≈ exp -s N U 1 + 4N u Q 2 (8)
where N is population size. Figure 2 shows that equation 8 provides accurate predictions when a 2 s is very small (left part of the figures). Furthermore, combining equations 7 and 8 (solid curves) gives good predictions for all values of a 2 s as long as they remain small (roughly, < 10 -2 for the parameter values of Figure 2) so that the Gaussian approximation holds. For higher values of a 2 s , simulations indicate that distributions of phenotypic values depart from Gaussian distributions: in particular, the fourth and sixth cumulants K 4 and K 6 are positive and increase as a 2 increases (while the third and fifth cumulants stay close to zero), for all values of Q (results not shown). In Supplementary File S3, we show that using a Gram-Charlier expansion to approximate the distribution of phenotypic effects yields an expression for W in terms of the genetic variance V G and of higher cumulants (K 4 , K 6 ). Furthermore, these higher cumulants can in turn be expressed in terms of the genetic variance, assuming that p j q j is small at each locus (rare-alleles approximation, e.g. [START_REF] Barton | Adaptive landscapes, genetic distance and the evolution of quantitative characters[END_REF][START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF]: one obtains in particular

K 4 ≈ a 2 V G [3m/ (2 + m)], K 6 ≈ a 4 V G [15m/ (4 + m)]
, which fits well with the simulation results (see Figure 1 in Supplementary File S3). Finally, the methods developed by Turelli and Barton (1994) can be used to calculate V G at equilibrium in this non-Gaussian regime. As shown in Supplementary File S3, the result (which has to be obtained numerically) fits well with the simulation results as long as a 2 s is not too large, and Q ≥ 2; however, we could not find any simple analytical expression for W in this non-Gaussian regime. From the expressions of K 4 and K 6 above, one predicts that the Gaussian regime should be left more rapidly (as a 2 increases) for larger m, since K 4 and K 6 are increasing functions of m; Figure 3 (left) confirms that this is indeed the case. However, Figure 3 right indicates that for a fixed s, W is relatively insensitive to m.

Figure 4 explores the effect of linkage on mean fitness: for U = 0.5 and n = 50, the effect of linkage disequilibria on W remains small as long as the harmonic mean recombination rate r H is not too small (with L = 10000 loci, r H ≈ 0.053 for a map length of 1 Morgan, while r H ≈ 0.3 for a map length of 10 Morgans -see Appendix B). In this regime, increasing linkage tends to reduce mean fitness (through an increase in genic variance V g ), and equation 7 matches well with the simulation results

(except for Q = 1 where the fit is less good). For tighter linkage, however (roughly, R < 1) simulations results depart from the prediction from equation 7 (which assumes that the contribution of linkage disequilibria to the genetic variance remains small) and show that W increases as recombination decreases: overall, linkage thus has a non-monotonic effect on mean fitness. Interestingly, equation 7 indicates that the contribution of linkage disequilibria becomes more important when the number of traits under selection (n) is small. Figure 5 confirms this result, showing that lower values of n lead to a stronger decrease in W for both s = 0.01 and s = 0.1, while the effect is less marked (and may even be reversed) when s = 0.001. This effect of n can be understood as follows: mean fitness decreases when the sum of genetic variances for the different traits n i=1 V G,i increases. Due to the symmetry of our model, the two components of the genetic variance (the genic variance V g,i , and the contribution of linkage disequilibria d i , see Appendix A) should be the same for all traits. As shown in Supplementary File S2 (see also Turelli, 1985;[START_REF] Bürger | Mutation-selection balance and the evolutionary advantage of sex and recombination[END_REF], when the effects of linkage disequilibria are neglected the genetic variance for each trait at equilibrium is proportional to 1/n (roughly, p j q j at each locus scales with 1/m, since the strength of selection at each locus increases with m; however, a proportion m/n only of loci contribute to trait i: therefore, the effect of m cancels and V g,i scales with 1/n). As a consequence W (which depends on the product nV G ) is independent of n. The contribution of associations between loci depends on all pairwise linkage disequilibria and is proportional to V g 2 (see Supplementary File S2) and therefore to 1/n 2 . As a consequence, the overall contribution of linkage disequilibria to mean fitness is proportional to 1/n. However, pairwise linkage disequilibria may not be the only cause Effects of population size and number of loci. Many of the results shown above assume a large population size and large number of loci. Figure 6 explores the effects of varying population size: in the case of a sexual population (Figure 6A), smaller population sizes tend to decrease W when the average fitness effect of mutations is sufficiently strong (right part of the figure). However, for lower values of s one observes the opposite effect, with higher values of mean fitness for lower population sizes. With a Gaussian fitness function (Q = 2, as in Figure 6), a diffusion model can be used to predict the equilibrium genetic variance under selection, mutation and drift, assuming that mean phenotypes stay at the optimum (Bulmer, 1972). In Supplementary File S4, we obtain an approximation for mean fitness (in terms of hypergeometric functions) valid for N u 1, which fits well with simulation results for N = 500 (dotted curve in Figure 6A). In contrast to the sexual case, the effects of population size are less marked in the case of asexual organisms (Figure 6B), and W always decreases as N decreases.

These effects of finite population size can be understood as follows. Genetic drift affects mean fitness through two different effects: an effect on the genetic variance maintained at equilibrium, which is usually (but not always) lower in smaller populations, and an effect on mean phenotypes, which may be displaced from the optimum due to drift.

In the case of a Gaussian fitness function (Q = 2), these two effects can easily be separated: indeed, log-fitness equals

-n i=1 (z i -z * i ) 2 / (2V S )
, where z i is the value of phenotype i and z * i the optimal value for this phenotype (fixed to zero in our model).

Therefore, the average log-fitness is given by

-n i=1 V G,i + (z i -z * i ) 2 / (2V S ), where
V G,i is the genetic variance for trait i and z i the average value of trait i (see also Charlesworth, 2013). Assuming that ln W ≈ ln W , mean fitness can thus be written as a product of two terms:

W ≈ exp - n i=1 V G,i 2V S exp - n i=1 (z i -z * i ) 2 2V S . (9)
The first term (denoted hereafter W V G ) shows that increasing the variance of phenotypic traits tends to decrease mean fitness, while the second term (denoted hereafter W z ) shows that displacing mean phenotypes from the optimum also decreases W .

Simulations show that injecting equilibrium values of genetic variances V G,i and mean phenotypes z i (measured in the simulations) into equation 9 accurately predict W for all parameters tried with Q = 2, as long as log 10 (s) ≤ -1 (results not shown). Note that all the mathematical derivations performed in the different Supplementary Files assume that mean phenotypes are at the optimum, and therefore that W z = 1. Supplementary Figure S2 shows measures of W V G and W z for the different values of N and s explored in Figure 6. These results can be summarized as follows: in a sexual population (Figure 6A), W z stays close to 1 in most cases, but may be slightly less than 1 (around 0.97 -0.98) for N = 500: mean phenotypes thus generally stay close to their optimal values, and the mutation load is mainly generated by the term W V G .

When mutations tend to have strong fitness effects (log 10 (s) ≥ -2), drift tends to increase genetic variances V G,i (by increasing p j q j at each locus j above its value at mutation-selection balance), and thus decreases W V G -note that this effect is taken into account by the diffusion model (equation 6 in Supplementary File S4). When mutations tend to have weak fitness effects (log 10 (s) ≤ -3) however, drift tends to reduce V G,i by bringing the different loci closer to fixation, thereby increasing W V G (and hence W ). By contrast, in an asexual population drift always reduce the genetic variance (W V G increases as N decreases), but also has a substantial effect on mean phenotypes (W z decreases as N decreases). Because the second effect is stronger than the first, the overall effect of drift is to reduce the mean fitness of asexuals.

Our model is somehow artificial, however, as we assumed that the fitness optimum is located at the origin (z * i = 0 for all i). Due to the symmetry of the model, one expects that mean phenotypes should be at the origin (and thus at the optimum) in the regime where dynamics are mainly driven by mutation and drift (low s), when the number of loci is large: for example if a large, random proportion of loci are fixed for allele 1 and the other loci fixed for allele 0, mean phenotypes should be close to zero since the sum of all mutational effects (corresponding to the effects of fixed alleles 1) converges to the average effect of mutations, which is zero. It seems artificial, however, to assume that optimal trait values are precisely the values towards which the population should converge at mutation-drift equilibrium, and it seems more likely that in many cases, mean phenotypes at mutation-drift equilibrium should be far from the optimal values (decreasing the value of W z in the mutation-drift regime). In order to explore that, we modified our simulation program so that the fitness optimum is not located at the origin of the phenotypic space (z i = 0 for all i). To do this, we select a genotype sufficiently distant from the origin (in phenotypic space) so that when setting optimal phenotypic values z * i to the values coded by this genotype, fitness at the origin is less than 0.1. At the start of the simulation, the population is fixed for the genotype corresponding to the fitness optimum. At mutation-drift balance, the population is still expected to converge to the origin (due to the symmetry of the distribution of mutational effects), where fitness is lower than 0.1. Results are shown on Figures 6C and 6D in the case of a sexual and asexual population, respectively. As can be seen on these Figures, the position of the fitness optimum has little effect as long as selection is sufficiently strong (log 10 (s) ≥ -3), since in this case selection maintains the population near the optimum. For lower s, however, mutation and drift tend to bring populations away from the optimum (and closer to the origin), which decreases W z and thus decreases mean fitness (see also Supplementary Figure S2), this effect being more pronounced in asexual than in sexual populations (in sexuals, this effect becomes important only for s = 10 -5 ). Note that we could not obtain results for s = 10 -6 , because with such low mutational effects (and with 10000 loci) it was not possible to find a genotype sufficiently distant from the origin so that fitness at the origin is lower than 0.1 when setting the optimum to the phenotypic values coded by this genotype.

Interestingly, in the sexual case (Figure 6C) one still observes an intermediate range of s (s = 10 -4 , 10 -3 ) where decreasing population size increases mean fitness: in this regime, drift tends to bring each locus closer to fixation (thus reducing V G,i ), but selection still maintains mean phenotypes z i close to the optimum. Decreasing s (or decreasing N further, as shown in Supplementary Figure S2 for N = 100) generates departures of mean phenotypes from the optimum, decreasing W z (and thus W ).

Finally, Supplementary Figure S3 shows that changing the number of loci L (for a fixed U ) has generally stronger effects in the case of sexual populations, where in the low s regime, decreasing the number of loci decreases the genetic variance (and thus increases W ).

DISCUSSION

As we have seen in introduction, different forms of models have been proposed

to predict the overall effect of recurrent mutation on the mean fitness of populations.

Typically, models representing deleterious mutations (without including a phenotypic dimension) do not consider possible compensatory effects among mutations, and assume that epistasis is the same for all pairs of mutations [START_REF] Kimura | The mutational load with epistatic gene interactions in fitness[END_REF][START_REF] Bürger | Mutation-selection balance and the evolutionary advantage of sex and recombination[END_REF]. Using this type of model, [START_REF] Kimura | The mutational load with epistatic gene interactions in fitness[END_REF] predicted that epistasis should have no effect on the mean fitness of asexuals (which should remain close to e -U ), while it should increase the fitness of sexuals if it is negative, and decrease W if it is positive. On the other hand, models representing stabilizing selection on a set of quantitative phenotypic traits (such as Fisher's geometric model) provide a natural way of incorporating compensatory effects among mutations and distributions of selection coefficients and epistatic effects. In this paper, we used a general model that allows us to play with the number of traits under selection (n), the number of traits affected by a given mutation (m) and the shape of the fitness peak (Q). Interestingly, this last parameter was already shown by others to affect the average sign of epistasis among pairs of mutations: epistasis is positive on average when the fitness peak is sharper than a Gaussian function (Q < 2), while it is negative when the fitness peak is flatter than a Gaussian (Q > 2), and equal to zero when fitness is exactly Gaussian [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF][START_REF] Gros | The evolution of epistasis and its links with genetic robustness, complexity and drift in a phenotypic model of adaptation[END_REF]. Furthermore, the average deleterious effect of mutations at the optimum s can be expressed as a simple function of the variance of scaled phenotypic effects of mutations a 2 s and the parameters m and

Q (equation 3).
In the case of an asexual population, we found that mean fitness stays indeed close to e -U when mutations tend to have strong fitness effects; however, W > e -U under weaker fitness effects of mutations (Gaussian regime), due to the possibility of compensatory effects among mutations at different loci. For a fixed s, mean fitness is higher under positive epistasis (Q < 2) than under negative epistasis (Q > 2); this is probably due to the fact that selection for compensatory mutations is stronger when the fitness peak is sharp, keeping the population closer to the fitness optimum. In contrast, under negative epistasis the fitness landscape presents a plateau around the optimum, and at equilibrium many genotypes are located near the edge of this plateau, where the effect of deleterious mutations tends to be strong. Finally, we showed that under this regime mean fitness can be expressed as a simple function of U , s, n and Q (equation 5), which indicates that W decreases as the number of traits under selection n increases.

In a sexual population, mean fitness is also at e -U when mutations have very strong fitness effects: this is expected, as mutations stay at very low frequency and epistatic interactions among mutations have thus little effect. For weaker mutational effects, [START_REF] Kimura | The mutational load with epistatic gene interactions in fitness[END_REF]'s prediction is verified: positive epistasis decreases W , while negative epistasis increases it; furthermore, we showed that W converges to a very simple expression (e -2U Q ) when the phenotypic effect of mutations is sufficiently weak, so that distributions of phenotypic values in the population are nearly Gaussian. Although this expression assumes that linkage disequilibria among loci can be neglected, effects of genetic associations can be computed, showing that linkage disequilibria tend to reduce W as long as they remain weak, the importance of this effect increasing as the number of selected traits n decreases. Contrarily to results obtained previously assuming N U 1 [START_REF] Hartl | Towards a theory of evolutionary adaptation[END_REF][START_REF] Poon | Compensating for our load of mutations: freezing the meltdown of small populations[END_REF][START_REF] Tenaillon | Quantifying organismal complexity using a population genetic approach[END_REF][START_REF] Lourenço | Complexity, pleiotropy and the fitness effect of mutations[END_REF], our model predicts a higher mean fitness when the number of traits is large, through this effect on linkage disequilibrium. Figure 5 indeed shows important reductions in fitness when the number of traits decreases (for moderate to strong fitness effects of mutations, of the order 0.01 to 0.1), which may also involve higher-order genetic associations (since distributions of phenotypic values depart from Gaussian distributions for very low values of n). Rather than the "cost of complexity" often described at low mutation rate, we thus predict a "cost of simplicity" when the overall mutation rate is large (as shown by equation 7, this effect is proportional to U 2 and should thus become negligible when the overall mutation rate is small). Note that this cost does not depend on the degree of pleiotropy (m) of individual mutations, but only on the total number of traits under selection. Interestingly, these results do not depend on the precise shape of the distribution of mutational effects (which was Gaussian in the simulations), since the derivation of equations 6 and uncorrelated mutational effects on the different traits, however, and it would be interesting to extend our analysis by incorporating mutational and/or selective correlations among traits -as shown by [START_REF] Martin | A general multivariate extension of Fisher's geometrical model and the distribution of mutation fitness effects across species[END_REF], the overall effect of such correlations is to reduce the effective number of phenotypic traits (see also [START_REF] Zhang | Multivariate stabilizing selection and pleiotropy in the maintenance of quantitative genetic variation[END_REF] for an analysis of the effects of correlations among traits). Furthermore, different ways of incorporating restricted pleiotropy could be considered: in particular, [START_REF] Welch | Modularity and the cost of complexity[END_REF] proposed a model in which different sets of traits are affected by different subsets of loci (modular pleiotropy), and it would be interesting to explore how this form of modularity would affect our results.

Regarding genetic architecture, an important assumption of our model is the fact that loci are biallelic. By contrast, Lande (1980)'s analysis considered an infinite number of alleles per locus (continuum-of-alleles model); in that case, assuming diagonal mutation and selection matrices (no mutational or selective covariance among traits) and neglecting linkage disequilibria, one obtains that the mean log-fitness should be -n 2 LU s (where L is the number of loci) in the case of a Gaussian fitness function (this stems from the fact that the contribution of each locus to the genetic variance is √ u a 2 V S , where u is the per-locus mutation rate). This contrasts with the result that mean log-fitness should be ≈ -U when Q = 2 in biallelic models, and indicates that assumptions about the number of alleles per locus may have important effects.

However, a key assumption of Lande's analysis is the fact that distributions of allelic effects in the population are Gaussian at each locus, which is much stronger than the assumption that distributions of overall phenotypes are Gaussian, and may not hold for realistic values of per-locus mutation rates and mutational variances (e.g., [START_REF] Turelli | Heritable genetic variation via mutation-selection balance: Lerch's zeta meets the abdominal bristle[END_REF]. It would thus be interesting to explore how sensitive are our results to assumptions regarding the number of alleles per locus. In particular, our assumption of biallelic loci should not hold when N u ∼ 1 (in which case many alleles may segregate at the same locus).

Finally, we found that the degree of pleiotropy of mutations (m) affects mean fitness mostly through its effect on s, in regimes where W is affected by s: in particular the Gaussian regime in asexual populations (equation 5), and the regime where selection at each locus is weaker than mutation and drift in sexual populations (equation 8). This contrasts with the result obtained by [START_REF] Lourenço | Complexity, pleiotropy and the fitness effect of mutations[END_REF] in the low mutation limit (N U 1), showing that the average fitness at equilibrium depends on the total number of selected traits but not on the degree of mutational pleiotropy.

Another difference concerns the effect of population size: although the load always increases as population size decreases when N U 1 [START_REF] Hartl | Towards a theory of evolutionary adaptation[END_REF][START_REF] Poon | Compensating for our load of mutations: freezing the meltdown of small populations[END_REF][START_REF] Tenaillon | Quantifying organismal complexity using a population genetic approach[END_REF][START_REF] Lourenço | Complexity, pleiotropy and the fitness effect of mutations[END_REF], we found that in sexuals, the mutation load may actually be lower at smaller population sizes in the weak selection regime (Figure 6). As we have seen, this occurs whenever the fitness effect of mutations is small (s = 10 -4 , 10 -3 in the case of Figure 6), so that diversity at each locus (p j q j ) at mutation-selection balance is substantial, in which case drift tends to reduce this diversity and therefore reduces the genetic variance V G -the same effect has also recently been discussed by Charlesworth (2013) in a 1-trait model. For very low fitness effects of mutations, however, mutation and drift may displace mean phenotypes away from the optimum and substantially reduce mean fitness (Figures 6C andD). This last effect may only occur for a restricted range of parameters, since when s becomes very low, very large numbers of loci are required for fitness to decrease significantly (see also Charlesworth, 2013).

In general, we found that asexual populations have a higher mean fitness than sexual populations, in agreement with previous models of stabilizing selection with a constant optimum (e.g., Charlesworth, 1993): this is due to the fact that genotypes coding for optimal trait values are preserved by asexual reproduction, but broken by recombination. Nevertheless, sexuals may have a higher mean fitness than asexuals when epistasis among mutations is negative on average (Q > 2), as long as s is not too small (so that W stays close to e -U in asexuals, but reach higher values in sexuals);

this occurs for example for Q = 6, when s ≥ 10 -4 (for the parameter values used in Figures 1 and2). Sexuals may also have a higher mean fitness when the effects of drift are important: in particular, Figure 6C and 6D show that deviations of mean phenotypes from their optimal values caused by drift are generally more important in asexual than in sexual populations (for low s). Although we only observed this last effect for restricted ranges of parameters, it may become more important in the case of populations subdivided into many small demes. Previous models explored the evolution of recombination modifiers under stabilizing selection (Charlesworth, 1993;Barton, 1995), but ignored the effects of genetic drift: it would thus be interesting to extend these models to the case of finite (or spatially structured) populations. Exploring the effects of diploidy and of the mating system would be other possible extensions of the present work. Other parameters are as in Figure 1 and 2, with Q = 2. A, B: the optimum is located at the origin (in phenotypic space), that is, in z i = 0 for all traits i; C, D: the optimum is located away from the origin, so that the fitness of an individual at the origin is less than 0.1 (see text for more explanations). Solid curves on the left figures (sexual populations) correspond to the predictions from equation 8 for N = 50000, 5000 and 500 from left to right, while the dotted curves corresponds to equation 6 in Supplementary File S4 (for N = 500).

Curves on the right figures (asexual populations) correspond to equation 5. Finally, the dashed horizontal lines correspond to e -U .

  program (written in C++, and available upon request) represents a haploid population of N individuals with discrete generations. The genome of each individual is represented by a sequence of L bits (0 or 1) corresponding to the different loci. Phenotypic effects of mutation at each locus are sampled at the start of the simulation: the m traits affected by a given locus are sampled randomly (and independently for each locus) among the n selected traits, and the effect of mutation on each of these m traits is sampled from a Gaussian distribution with mean 0 and variance a 2 . At the start of a generation, phenotypic values are computed for each individual; from this, one obtains the fitness of the individual based on equation 1. For each individual of the next generation, two parents are sampled (the probability that an individual is sampled being proportional to its fitness); selfing is allowed if the same individual is sampled twice. To generate a recombinant chromosome, the number of cross-overs is sampled from a Poisson distribution with parameter R (genome map length), and the position of each cross-over along the chromosome is sampled from a uniform distribution (R = 0 corresponds to the case of an asexual population). Finally, the number of new mutations occurring in each genome is sampled from a Poisson distribution with parameter U , and the positions of mutant loci are sampled randomly; alleles are switched at mutant loci, from 0 to 1 or from 1 to 0 (mutation and back mutation thus

  to the fitness optimum. Simulations run for 50000 generations, equilibrium being reached during the first 30000 generations for most parameter values tested. Every 1000 generations, the program records the mean fitness of the population and the first 6 moments of the distribution of each phenotypic trait in the population. The results shown in the different figures correspond to averages over the last 20000 generations. Error bars were calculated using Hastings' (1970) batching method, dividing these 20000 generations into 4 batches of 5000 generations and calculating the standard error over these 4 batches; however, error bars were generally small relative to the size of symbols used in the figures and are thus not shown. All our simulation results (with the Mathematica commands used to generate the figures) can be found in the Supplementary Material.

  of the strong reduction of W shown on Figure 5 at low values of n (for s = 0.01 and 0.1): indeed, simulations indicate that for these low values of n, distributions of phenotypic traits in the population depart from Gaussian distributions (fourth and sixth cumulants become significantly different from zero), in which case the results derived in Supplementary File S2 (which assume Gaussian distributions of phenotypic traits) do not hold.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .

 12345 Figure 1. Mean fitness of an asexual population as a function of a 2 s and s. X-axes correspond to the mutational variance on scaled traits a 2 s = a 2 / (2V S ) 2 Q (left) and the

Figure 6 .

 6 Figure 6. Mean fitness of a sexual (A, C) and asexual (B, D) populations: effects

Table 1 :

 1 Parameters and variables of the model.

	n		number of phenotypic traits under selection
	m		number of traits affected by a single mutation
	Q		shape of the fitness peak (see equation 1)
	V S		strength of selection (see equation 1)
	a 2		variance of phenotypic effects of mutations
	a 2 s = a 2 / (2V S )	2 Q	variance of mutational effects on scaled traits
	L		number of loci affecting selected traits
	U		total rate of mutation per generation on selected loci
	u = U/L		mutation rate per locus
	s		mean deleterious effect of mutations (on log-fitness) at the
			fitness optimum (given by equation 3)
	R		genome map length (in Morgans)
	r H		harmonic mean recombination rate among pairs of loci
	N		population size
	z i		value of trait i in a given individual
	r ij		effect of mutation at locus j on trait i
	p j		frequency of allele 1 at locus j in the population
	D jk		linkage disequilibrium between loci j and k
	V G		genetic variance (variance of z i among individuals, which at
			equilibrium should be the same for all traits i)
	V

g genic variance (see Appendix A)

Acknowledgements. We thank Nick Barton, Harold de Vladar, Pierre-Alexis Gros, Guillaume Martin, Sylvain Glémin and two anonymous reviewers for helpful discussions and comments, and the bioinformatics and computing service of Roscoff's Bio-

Nationale de la Recherche (ANR-11-BSV7-013-03).

APPENDIX A: VARIANCE COMPONENTS

We call V G,i the variance of the value of trait i in the population: V G,i = Var [z i ]. Because we assume additive effects of the different loci on phenotypic traits (z i = j r ij X j , where the sum is over all loci j), V G,i can be decomposed into two terms (e.g., [START_REF] Lynch | Genetics and Analysis of Quantitative Traits[END_REF]

is the genic variance (p j being the frequency of allele 1 at locus j in the population), and d i = j =k r ij r ik D jk the effect of linkage disequilibria on the genetic variance

being the linkage disequilibrium between loci j and k). In the case of an asexual population, the genetic variance at equilibrium can be calculated by considering the whole genome as a single locus with a large number of alleles (Supplementary File S1). In the sexual case, V g,i and d i at equilibrium can be computed using the methods developed by Turelli and [START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF], extended to deal with multiple traits under selection (Supplementary File S2). Due to the symmetry of the model we expect that V G,i and V g,i should be the same for all traits at equilibrium, and are thus simply denoted V G and V g in the text.

APPENDIX B: HARMONIC MEAN RECOMBINATION RATE

Under the genetic setting described above, the harmonic mean recombination rate among pairs of loci r H (that appears in equation 7) can be computed as follows.

The genetic distance between adjacent loci (in Morgans) is R/ (L -1), and therefore the distance between two loci separated by i between-locus intervals is iR/ (L -1).

Furthermore, the number of different pairs of loci separated by i intervals is L -i.

Finally, the rate of recombination between two loci at genetic distance x (probability that an odd number of cross-overs occurs within the interval) is given by Haldane's mapping function: r (x) = 1 2 (1 -e -2x ), yielding:

which has to be evaluated numerically. For L = 10000 and R = 10, one obtains r H ≈ 0.3.


Camille Gervais

Diala Abu Awad

Denis Roze

Vincent Castric

Sylvain Billiard
email: sylvain.billiard@univ-lille1.fr

Genetic architecture of inbreeding depression and the maintenance of gametophytic self-incompatibility

Keywords: deleterious mutation, inbreeding depression, self-incompatibility, purging

Gametophytic self-incompatibility is a widespread genetic system, which enables hermaphroditic plants to avoid self-fertilization and mating with close relatives.

Inbreeding depression is thought to be the major force maintaining self-incompatibility; however, inbreeding depression is a dynamical variable that depends in particular on the mating system. In this paper we use multilocus, individual based simulations to examine the co-evolution of self-incompatibility and inbreeding depression within finite populations. We focus on the conditions for the maintenance of self-incompatibility when self-compatible mutants are introduced in the population by recurrent mutation, and compare simulation results with predictions from an analytical model treating inbreeding depression as a fixed parameter (thereby neglecting effects of purging within the self-compatible sub-population). In agreement with previous models, we observe that the maintenance of self-incompatibility is associated with high inbreeding depression and is facilitated by high rates of self-pollination. Purging of deleterious mutations by self-compatible mutants has little effect on the spread of those mutants as long as most deleterious alleles have weak fitness effects: in this case, the genetic architecture of inbreeding depression has little effect on the maintenance of self-incompatibility.

By contrast, purging may greatly enhance the spread of self-compatible mutants when deleterious alleles have strong fitness effects.

INTRODUCTION

Self-incompatibility (SI) is a genetic system that enables hermaphroditic plants to avoid self-fertilization and limit mating with close relatives by recognition and rejection by the pistil of pollen expressing cognate specificities. In many species, SI specificity is controlled by a single multiallelic locus, the S-locus. SI is widespread, found in more than 100 families of angiosperms [START_REF] Igic | Loss of self-incompatibility and its evolutionary consequences[END_REF] despite the fact that it entails a transmission disadvantage. Indeed, a self-compatible (SC) mutant occurring in a SI population should benefit from a higher number of potential mates, and from a transmission advantage through self-fertilization. This last advantage is equal to 2/3, and is thus higher than the 50% advantage of an allele coding for selfing in an outcrossing population [START_REF] Fisher | Average excess and average effect of a gene substitution[END_REF], because only the SC pollen contributes to the selfed offspring in heterozygous individuals for SC. Recently, it has been argued that SI may be maintained by selection acting at the species level, due to the fact that SI species diversify at higher rates [START_REF] Goldberg | Species selection maintains self-incompatibility[END_REF]. This form of selection only acts in the long term, however, and implies that SC mutations should occur very rarely.

The main advantage of SI that could possibly explain its maintenance in the short term is the avoidance of inbreeding depression (e.g., Charlesworth and Charlesworth, 1987).

However, [START_REF] Charlesworth | Evolution and breakdown of S-alleles systems[END_REF] showed that for SI to be maintained in the presence of self-compatible (SC) mutants, inbreeding depression has to be high; its minimal value depends on the number of SI alleles segregating in the population and on the selfing rate of SC mutants, but is often close to 2/3 when the number of SI alleles is large -this minimal value is higher and may even reach 1 when the number of SI alleles is not large and the selfing rate of self-compatible mutants is small to moderate, as SC mutants also benefit from a higher number of potential mates. Furthermore, the minimal value of inbreeding depression needed to maintain SI may be underestimated by [START_REF] Charlesworth | Evolution and breakdown of S-alleles systems[END_REF]'s model because inbreeding depression was treated a fixed variable, therefore neglecting the effect of purging of deleterious alleles by the self-fertilizing SC individuals. Purging is expected to facilitate the spread of SC mutants, as these mutants tend to better eliminate partially recessive deleterious alleles, thereby reducing the magnitude of inbreeding depression experienced by their selfed offspring [START_REF] Uyenoyama | Coevolution of self-fertilization and inbreeding depression. I. Mutation-selection balance at one and two loci[END_REF][START_REF] Glémin | How are deleterious mutations purged? Drift versus nonrandom mating[END_REF].

This was confirmed by a deterministic simulation model by [START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF], assuming that a given proportion of self-pollen lands on the stigma ("mass-action" pollination model, [START_REF] Holsinger | Mass-action models of plant mating systems: the evolutionary stability of mixed mating systems[END_REF] and that inbreeding depression results from recessive lethal mutations segregating at a very large (effectively infinite) number of loci [START_REF] Kondrashov | Deleterious mutations as an evolutionary factor. II. Facultative apomixis and selfing[END_REF]. In this model, invasion of a SI population by a SC mutant (i.e. the breakdown of SI) is easier than expected from [START_REF] Charlesworth | Evolution and breakdown of S-alleles systems[END_REF]'s results, in particular when the selfing rate of SC mutants is moderate to high (so that purging can occur). [START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF] also observed that in a small region of parameter space (namely, for high values of inbreeding depression and low selfing rates), the population may reach a stable, polymorphic equilibrium where both SI and SC individuals are present. This observation is important in the context of how new self-incompatibility specificities may arise through self-compatible intermediates, assuming that compensatory mutations may secondarily restore a novel SI functionality [START_REF] Uyenoyama | On the origin of self-incompatibility haplotypes: transition through self-compatible intermediates[END_REF][START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF].

Although [START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF]'s model showed that considering the joint dynamics of SI and inbreeding depression may strongly affect predictions concerning the maintenance of SI, the fact that they considered lethal mutations may overestimate the importance of purging in situations where inbreeding depression is (at least partly) generated by mildly recessive deleterious alleles. Indeed, analytical models have shown that the advantage of a selfing modifier through purging is expected to increase with the strength of selection against deleterious alleles [START_REF] Uyenoyama | Coevolution of self-fertilization and inbreeding depression. I. Mutation-selection balance at one and two loci[END_REF][START_REF] Epinat | The evolution of assortative mating and selfing with in-and outbreeding depression[END_REF], which was also observed in deterministic simulations representing the spread of a mutation affecting the selfing rate [START_REF] Charlesworth | Inbreeding depression, genetic load, and the evolution of outcrossing rates in a multilocus system with no linkage[END_REF]. More recently, Porcher and Lande (2013) showed that the effect of purging on the spread of a mutation affecting selfing is much lower when the strength of selection against deleterious alleles is mild (s = 0.05) than when it is strong (s = 1). However, Porcher and Lande (2013) only considered weak-effect selfing modifiers (changing the selfing rate by 10 -6 ), and the effect of purging may be stronger in the case of a mutation having a large effect on the selfing rate [START_REF] Charlesworth | Inbreeding depression, genetic load, and the evolution of outcrossing rates in a multilocus system with no linkage[END_REF]) such as a mutation disrupting SI. Because the results of [START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF] showed that purging may strongly limit conditions for the maintenance of SI in the presence of lethal recessive mutations, it is important to assess the generality of this result (in particular to situations where a substantial proportion of mutations are mildly deleterious) in order to better understand how self-incompatibility can be maintained in natural populations (note that [START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF] also considered situations where inbreeding depression is partly due to mildly deleterious mutations, but this component of inbreeding depression is not dynamic in their model).

All previous simulation models explicitly considering the genetic basis of inbreeding depression [START_REF] Charlesworth | Inbreeding depression, genetic load, and the evolution of outcrossing rates in a multilocus system with no linkage[END_REF][START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF]Lande, 2005, 2013) were based on [START_REF] Kondrashov | Deleterious mutations as an evolutionary factor. II. Facultative apomixis and selfing[END_REF]'s model, representing the dynamics of deleterious alleles at an infinite number of unlinked loci, in an infinite population. In principle, genetic linkage should increase the effect of purging, by increasing the association between self-compatible alleles and purged genetic backgrounds; however, whether the effect of linkage is substantial for realistic values of genomic recombination rates is unclear.

Finite population size is expected to affect the number of SI alleles maintained in the population -which depends both on population size and on the mutation rate towards new SI alleles [START_REF] Yokoyama | The expected number of self-incompatibility alleles in finite plant-populations[END_REF]. Since conditions for the maintenance of SI depend on the number of SI alleles segregating [START_REF] Charlesworth | Evolution and breakdown of S-alleles systems[END_REF][START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF][START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF], they should thus be affected by population size. Furthermore, the number of SI alleles may also change as self-compatible mutants increase in frequency, which may in turn affect the conditions needed to maintain SI and SC individuals at a polymorphic equilibrium [START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF].

In this paper, we use a multilocus individual-based simulation program to explore the conditions for the maintenance of SI when inbreeding depression is generated by deleterious alleles segregating at a large number of partially linked loci in a finite population. Our model also differs from the previous models cited above by assuming that inbreeding depression affects both seed and pollen production, which in turn affects the selfing rate of SC individuals (as it depends on the quantity of self-pollen relative to the quantity of pollen received from other individuals). We show that in many cases, conditions for the maintenance of SI are similar to those obtained from an analytical model assuming fixed inbreeding depression (that is, the effect of purging on the spread of SC mutants remains small), unless a substantial proportion of deleterious mutations have strong fitness effects. Linkage has only little effect for the parameter values tested, as long as the mean number of cross-overs per genome (i.e., genetic map length) is not too small. Finally, we almost never observe a polymorphic equilibrium, which is probably due to the fact that the number of SI alleles decreases as SC individuals increase in frequency, further enhancing the advantage of self-compatibility.

METHODS

We consider a population with a gametophytic self-incompatibility (GSI) system, i.e. fertilization is possible only if the specificity expressed by the pollen is different from the two specificities codominantly expressed in the style. Specificities are subject to negative frequency-dependent selection, since pollen bearing rare specificities can fertilize more individuals than pollen bearing more common specificities [START_REF] Wright | The distribution of self-sterility alleles in populations[END_REF]. We assume that self-incompatibility is coded by a single S-locus with many alleles (denoted S i ): a plant with genotype S i S j (i = j) is self-incompatible and can be fertilized by pollen of genotype S l , with l = i and l = j. We also consider a mutant, self-compatible allele S C segregating at the same locus: pollen carrying S C can fertilize all plants, while a plant with genotype S i S C is partially self-compatible (through its S C pollen) and can be fertilized by any pollen whose genotype is different from S i .

Finally, S C S C individuals are fully self-compatible and can be fertilized by any pollen.

A parameter α measures the proportion of pollen produced by a plant that stays on the same plant ("self-pollen"), leading to self-fertilization if it carries the S C allele. We assume that selfed offspring suffer from inbreeding depression, generated by partially recessive deleterious alleles segregating at a number of different loci. Several selective forces may affect the frequency of the S C allele: (i) automatic transmission advantage through selfing, when α > 0, (ii) transmission advantage through outcrossing, as S C pollen can fertilize all plants, (iii) negative consequences of increased homozygosity due to selfing on the mean fitness of offspring (inbreeding depression) and (iv) indirect benefits stemming from a better elimination of deleterious alleles as a consequence of increased homozygosity (purging). This last effect occurs more rapidly when deleterious alleles have stronger effects (e.g., [START_REF] Charlesworth | The genetics of inbreeding depression[END_REF]. S C is thus favoured by effects (i), (ii) and (iv), and disfavoured by effect (iii); but note that purging also tends to reduce inbreeding depression. In the following, we use a multilocus, individual-based simulation model to represent explicitly the genetic causes of effects (iii) and (iv) (inbreeding depression and purging). Results on the spread of S C are compared to the predictions of a simple analytical model where effect (iv) is ignored by treating inbreeding depression as a fixed parameter.

Analytical model. Our analytical model represents an infinite population with discrete, non-overlapping generations. We assume that n + 1 alleles segregate at the S-locus: n self-incompatible alleles and the self-compatible S C allele. We denote x 1 and x 2 the frequencies of S C S C individuals produced by selfing (x 1 ) and by outcrossing (x 2 ), and x 3 and x 4 the frequencies of S i S C individuals produced by selfing (x 3 ) and by outcrossing (x 4 ), where S i can be any SI allele. We suppose that selfed individuals produce fewer gametes (inbreeding depression), the number of male and female gametes produced by outcrossed individuals being proportional to W 2 = W 4 = 1, while the number of gametes produced by selfed individual is proportional to

W 1 = W 3 = 1 -δ;
mean fecundity W is thus given by 1 -δ(x 1 + x 3 ). The selfing rate a i (proportion of selfed seeds) of an individual of type i (1, 2, 3 or 4) is given by:

a i = α θ i W i αW i + (1 -α)W (1)
where θ i = 1 for i = 1, 2 and θ i = 1/2 for i = 3, 4. The numerator of equation 1 represents the quantity of compatible self-pollen (bearing allele S C ), while the denominator is the total quantity of pollen received by the individual. The frequency of allele S C among gametes is given by q = i θ i W i x i /W , while, by symmetry, each SI allele is present in frequency p = (1 -q)/n. From this, frequencies at the next generation are given by:

W x 1 = W 1 a 1 x 1 + W 2 a 2 x 2 + 1 2 W 3 a 3 x 3 + 1 2 W 4 a 4 x 4 W x 2 = W 1 (1 -a 1 ) q x 1 + W 2 (1 -a 2 ) q x 2 + [W 3 (1 -a 3 ) x 3 + W 4 (1 -a 4 ) x 4 ] q 2 (1 -p) W x 3 = 1 2 W 3 a 3 x 3 + 1 2 W 4 a 4 x 4 W x 4 = W 1 (1 -a 1 ) (1 -q) x 1 + W 2 (1 -a 2 ) (1 -q) x 2 + 1 2 W 3 (1 -a 3 ) x 3 + 1 2 W 4 (1 -a 4 ) x 4 + (1 -x 1 -x 2 -x 3 -x 4 ) q 2 (1 -p) . (2)
In the Supplementary Material, we use a local stability analysis to determine the values of α, n and δ for which S C increases in frequency when rare (which involves solving a fourth-order equation in δ numerically).

Multilocus simulations. Individual-based, multilocus simulations were used to explore the conditions for the maintenance of SI when deleterious alleles segregate at a large number of partly linked loci, and when the number of SI alleles evolves freely by mutation, selection and drift. The simulation program (written in C++, and available from Dryad) corresponds to a modified version of the program described in Roze and [START_REF] Michod | Deleterious mutations and selection for sex in finite, diploid populations[END_REF], representing a population of N diploid individuals (the parameters used in the simulation model are summarized in Table 1). Each individual possesses two copies of a linear chromosome, and a S-locus located at the mid-point of the chromosome. We assume that a maximum of k + 1 alleles can segregate at the S-locus:

k self-incompatible alleles S i (1 ≤ i ≤ k)
W i = (1 -hs) N he (1 -s) N ho ,
where N he and N ho are the number of mutations in the heterozygous and homozygous state within its genome. If the maternal plant carries at least one self-compatible allele S C , its selfing rate a i is calculated as:

a i = α γ ii W i α γ ii W i + 1-α N -1 j =i γ ij W j (3)
where γ ij is 0, 1 or 2, and represents the number of individual j's S-alleles which are compatible with individual i. In the case of selfing, the offspring's genome is formed from two of i's recombinant gametes; otherwise an individual j is sampled randomly with a probability proportional to W j , and contributes as a father only if one of its S-alleles is compatible with those of the mother (if not, another individual is sampled until a compatible partner is found). To form a recombinant gamete, the number of cross-overs occurring along the chromosome is drawn from a Poisson distribution with parameter L (genome map length, in Morgans) and the position of each cross-over is sampled from a uniform distribution.

At the beginning of each simulation individuals are free of deleterious mutations and heterozygous for randomly sampled SI alleles. During the first 7,500 generations, the number of SI alleles segregating in the population is allowed to reach equilibrium by considering only mutation between SI alleles (U SC = 0, U = 0). Deleterious mutations are then introduced and allowed to reach mutation-selection balance over the next 7,500 generations. Finally, during the next 500,000 generations, self-compatible mutants are also introduced. Because loss of SI is irreversible, simulations are stopped after 50,000 generations if S C is fixed in the population. Every 100 generations, different variables are measured from the population: the frequency of allele S C , the effective number of SI alleles present and the level of inbreeding depression. The effective number of SI alleles n e (measured before allele S C is introduced) corresponds to the number of alleles that would yield the same genetic diversity at the S-locus if all alleles were present in frequency 1/n e ; it is calculated as n e = 1/ k i=1 p 2 i , where p i is the frequency of allele S i . Inbreeding depression is measured as δ = 1 -W s /W o , where W s and W o are the mean fecundities of selfed and outcrossed offspring, respectively (estimated by creating 100 selfed and 100 outcrossed offspring from randomly sampled parents, without taking into account the compatibility between their S-alleles). For each simulation run, δ and n e are averaged over the last 50 samples before introduction of self-compatible mutants, while the frequency of S C is averaged over the last 300 samples of the simulation (last 30,000 generations). The minimal value of δ necessary to maintain SI in the population is determined by running simulations with increasing values of U for each set of parameters values (U SI , U SC , s, h, N , L and α). When the frequency of S C stays lower than 0.05 throughout the simulation, SI is considered maintained. to SI), one expects that allele S C should necessarily become fixed after a sufficiently long time. However, changing the number of generations with U SC > 0 to 10 5 or to 10 6 did not lead to significant differences in the threshold values of δ (results not shown), suggesting that as one enters the area above the points in Figure 1, the expected fixation time of S C quickly reaches extremely high values (i.e. SI is stably maintained). In the same vein, increasing the mutation rate towards S C (U SC ) from 10 -5 to 10 -4 has very little effect on the results (not shown).

RESULTS

Mildly deleterious mutations.

Curves on Figure 1 correspond to the predictions derived from the analytical model (neglecting the effect of genetic associations between the S-locus and selected loci), where the number of SI alleles n is set to the average effective number of alleles n e measured in simulations corresponding to the critical δ. Overall, simulation results are qualitatively and quantitatively consistent with these analytical expectations: maintenance of SI is always observed when inbreeding depression is high (δ > 2/3), the critical δ being lower when the self-pollination rate α is higher, in agreement with previous results neglecting effects of purging [START_REF] Charlesworth | Evolution and breakdown of S-alleles systems[END_REF][START_REF] Uyenoyama | On the origin of self-incompatibility haplotypes: transition through self-compatible intermediates[END_REF]. The effect of population size N on the critical δ is due to the fact that the effective number of SI alleles is reduced in smaller populations through the loss of low-frequency alleles by drift (e.g., [START_REF] Yokoyama | The expected number of self-incompatibility alleles in finite plant-populations[END_REF]: on average n e at the critical δ equals 11.9, 23.0 and 36.4 for N = 500, 2000 and 5000, respectively. Lower effective numbers of SI alleles favour the spread of S C by increasing the transmission advantage of S C through outcrossing, as S C pollen never encounters incompatible pistils and can fertilize every potential mate in the population [START_REF] Charlesworth | Evolution and breakdown of S-alleles systems[END_REF][START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF][START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF].

Simulations for different values of U SI (rate of mutation towards new SI alleles) show that increasing U SI (with N U SI = 0, 0.02 and 0.2) has similar effects as increasing N (Figure S1).

Overall, the match between the simulation results and predictions from the analytical model indicates that purging has little effect on the spread of S C for these parameter values; however, the analytical model systematically underestimates the critical δ for SI to be maintained, which is probably a consequence of purging (the discrepancy remains slight, but becomes more apparent as α increases, and thus as the selfing rate of self-compatible individuals becomes more important). The effect of purging on the spread of S C should be more important when selection against deleterious alleles is stronger (e.g., [START_REF] Charlesworth | Inbreeding depression, genetic load, and the evolution of outcrossing rates in a multilocus system with no linkage[END_REF], as deleterious alleles are eliminated more rapidly when present in homozygotes, and as the benefit of being associated with chromosomes carrying fewer deleterious alleles is stronger. In agreement with this prediction, Figure 2 shows that discrepancies between analytical and simulation results becomes more important for higher values of s and h. One can also see from Figure 2 that, at the critical δ, s has little effect on the effective number of SI alleles n e maintained in the population (the curves on Figure 2A are almost superposed), while higher values of h lead to lower values of n e . This is likely due to the fact that n e is affected by background selection (reduction in diversity due to selection at linked sites, e.g., [START_REF] Charlesworth | The effect of deleterious mutations on neutral molecular variation[END_REF]. Indeed, background selection is stronger when the deleterious mutation rate U is higher (e.g., [START_REF] Hudson | Deleterious background selection with recombination[END_REF], and higher values of U are needed to reach the critical δ when h is increased, leading to stronger background selection effects at the critical δ. By contrast, δ is only weakly dependent on s (as long as population size is sufficiently large, e.g., [START_REF] Bataillon | Inbreeding depression due to mildly deleterious mutations in finite populations: size does matter[END_REF], and s has thus little effect on the strength of background selection at the critical δ. Finally, we find that varying map length L has relatively little effect on the results as long as it is sufficiently large (roughly, L > 5 using our default parameter values shown in Table 1). Figure S2 shows that reducing L decreases the effective number of SI alleles (background selection) and increases selection for S C through purging (higher discrepancy between analytical and simulation results when L is lower), both effects reducing the parameter range where SI is maintained.

DISCUSSION

Self-incompatibility (and especially gametophytic SI) is widespread among angiosperms, despite the fact that self-compatible mutants should benefit from a direct transmission advantage. Different models have shown that, provided that it is sufficiently high, inbreeding depression can allow the stability of SI despite the recurrent occurrence of SC mutants [START_REF] Charlesworth | Evolution and breakdown of S-alleles systems[END_REF][START_REF] Uyenoyama | On the origin of self-incompatibility haplotypes: transition through self-compatible intermediates[END_REF][START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF][START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF]. However, these models differ in the way inbreeding depression is introduced: some models treat inbreeding depression as a fixed parameter, thereby neglecting the effects of purging within SC backgrounds [START_REF] Charlesworth | Evolution and breakdown of S-alleles systems[END_REF][START_REF] Uyenoyama | On the origin of self-incompatibility haplotypes: transition through self-compatible intermediates[END_REF][START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF], while the model by [START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF] explicitly represents the genetic architecture of inbreeding depression (infinite number of unlinked loci subject to recessive lethal mutations) and shows that purging can dramatically decrease the parameter range where SI is maintained, in particular when the selfing rate of SC mutants is high. The importance of this effect should, in principle, depend on the genetic basis of inbreeding depression, since purging may occur over just a few generations when deleterious alleles are highly deleterious, but much more slowly when mutations tend to have weak fitness effects. In this paper, we compare predictions from an analytical model assuming fixed inbreeding depression to the results of multilocus simulations where deleterious alleles occur along a linear genetic map, in order to assess the importance of purging on the spread of SC mutants. When inbreeding depression is mainly due to weak-effect mutations, purging has limited effects. In most cases, maintenance of SI mainly depends on the number of SI alleles segregating in the population, the rate of self-pollination and inbreeding depression, independently of population size and the genetic architecture of inbreeding depression -a similar result was obtained recently by Porcher and Lande (2013) on the evolution of self-fertilization through weak-effect modifiers. Note that our model does not incorporate pollen limitation, which would tend to favour the loss of SI [START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF]. However, our results relative to the effect of purging should remain valid in the presence of pollen limitation: pollen limitation adds direct selection for self-compatibility, but should not modify the indirect effect of deleterious alleles on SC mutants.

Whilst our current knowledge on the genetic basis of inbreeding depression remains fragmentary, several lines of evidence suggest an important role of mutations of small effects [START_REF] Carr | Recent approaches into the genetic basis of inbreeding depression in plants[END_REF][START_REF] Charlesworth | The genetics of inbreeding depression[END_REF], although a study by [START_REF] Fox | Experimental evolution of the genetic load and its implications for the genetic basis of inbreeding depression[END_REF] showed a rapid reduction of inbreeding depression after several generations of inbreeding, indicating a potentially important effect of strongly deleterious mutations. More experimental work is thus necessary to assess whether purging is susceptible to significantly affect the spread of self-compatible mutants within self-incompatible populations.

Finally, we almost never observe stable polymorphic equilibria involving both SI and SC alleles (except for restricted cases involving very high inbreeding depression and low rates of self-pollination). This stands in contrast to previous models involving infinite populations [START_REF] Uyenoyama | On the origin of self-incompatibility haplotypes: transition through self-compatible intermediates[END_REF][START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF][START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF], where wider regions of parameter space allowing polymorphism are observed. A possible explanation for this discrepancy is that, in infinite populations, the number of SI alleles stays constant and is not affected by the frequency of SC, while in our model the number of SI alleles decreases as SC increases in frequency (since the size of the SI sub-population decreases). The decrease in number of SI alleles tends to favour SC, whose frequency can further increase until reaching fixation. This result is consistent with the fact that SC alleles are rarely found in natural SI populations [START_REF] Stone | Molecular mechanisms underlying the breakdown of gametophytic self-incompatibility[END_REF], although some cases have been reported where self-incompatibility appears to be quantitative rather than qualitative, with some partially self-compatible alleles [START_REF] Mena-Ali | Segregation analyses of partial selfincompatibility in self and cross progeny of Solanum carolinense reveal a leaky Sallele[END_REF][START_REF] Paape | Evolutionary genetics of an S-like polymorphism in Papaveraceae with putative function in selfincompatibility[END_REF]. The lack of polymorphic equilibria should impose restrictions on the evolution of new SI specificities. In the present model we assume that new SI alleles appear in a single mutational step, but in reality the evolution of a new specificity involves at least two mutations: one affecting the protein expressed by the pollen and the other the receptor expressed by the pistil (both genes being part of the S-locus). Most scenarios for the evolution of new specificities rely on an intermediate step involving a SC mutant present at an intermediate frequency in the population (e.g., [START_REF] Uyenoyama | On the origin of self-incompatibility haplotypes: transition through self-compatible intermediates[END_REF][START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF], which should become more difficult in the absence of polymorphic equilibrium (unless the mutation rate at the S-locus is sufficiently high, so that a compensatory mutation can appear before SC reaches fixation). Extending our simulation model to explicitly represent the pollen and pistil components of the S-locus (in order to explore conditions for the evolution of new SI specificities with dynamical inbreeding depression) would be an interesting extension of the present work.

  and a self-compatible allele S C ; k is fixed at 100 in all simulations. Every generation, each S i allele mutates to any of the k -1 other SI alleles with probability U SI and to allele S C with probability U SC (loss of SI). There is no reverse mutation from S C to S i , i.e. loss of self-incompatibility is irreversible. Deleterious mutations occur at a rate U per haploid genome, that is, the number of new deleterious mutations on each chromosome is drawn from a Poisson distribution with parameter U , while the position of each new mutation along the chromosome is sampled from a (continuous) uniform distribution -the number of sites at which deleterious alleles may segregate is thus effectively infinite. All deleterious mutations have the same selection and dominance coefficients (s and h respectively), although we also consider cases where a given proportion of mutations are lethal (as detailed below). Reproduction occurs as follows: for each of the N individuals of the next generation the maternal parent is sampled randomly among all individuals of the previous generation, the probability of sampling parent i being proportional to its fecundity

  Figure1shows the minimal value of inbreeding depression δ necessary to maintain SI in the population as a function of the rate of selfpollination α, for different values of population size N . In all figures, the highest value of inbreeding depression (in the absence of S C ) for which we observed that S C invades the population (i.e., reaches frequency 0.05) is slightly below the points, at a distance equivalent to the size of points (not shown). In almost all simulations the frequency of S C at equilibrium was either close to zero or equal to 1, i.e. polymorphic equilibria for S C were almost never observed; we will return to this point in the Discussion. Note that because the loss of SI is irreversible in our model (no back mutation from SC

  Adding a proportion of nearly recessive lethal mutations. Using a deterministic model where mutations generating inbreeding depression are lethal and nearly recessive,[START_REF] Porcher | Loss of gametophytic self-incompatibility with evolution of inbreeding depression[END_REF] found much more stringent conditions for the maintenance of SI than those shown on Figures1-2. We modified our simulation program so that a proportion λ of deleterious mutations are nearly recessive lethals (selection and dominance coefficients: s l = 1 and h l = 0.02, respectively) while the other mutations (in proportion 1-λ) are mildly deleterious (s = 0.05 and h = 0.2). As shown by Figure3, the discrepancy between analytical and simulation results becomes stronger as λ increases, and the parameter range where SI is maintained is considerably reduced. In particular, allele S C invades the population for much larger values of inbreeding depression when the rate of self-pollination α is high, so that individuals carrying S C frequently self. The proportion of lethal mutations has almost no effect on the effective number of SI alleles maintained in the absence of S C , and therefore the analytical predictions are nearly unaffected by λ (curves on Figure3are nearly superposed).

Figure 1 .Figure 3 .

 13 Figure 1. Critical inbreeding depression above which the SC mutant cannot in-

Table 1 :

 1 Parameters and default values used in simulations.Rate of mutation fromS i to any S j (with j = i) U SI 10 -5Rate of mutation from S i to S C U SC 10 -4

	Description

Acknowledgements. This work was supported by the French Agence Nationale de la Recherche (ANR-11-BSV7-013).