Abu Diala

Awad

Denis Roze
email: roze@sb-roscoff.fr

Effects of partial selfing on the equilibrium genetic variance, mutation load and inbreeding depression under stabilizing selection

Keywords: adaptive landscape, epistasis, evolutionary quantitative genetics, multilocus population genetics, self-fertilization

The mating system of a species is expected to have important effects on its genetic diversity. In this paper, we explore the effects of partial selfing on the equilibrium genetic variance V g , mutation load L and inbreeding depression δ under stabilizing selection acting on a arbitrary number n of quantitative traits coded by biallelic loci with additive effects. Overall, our model predicts a decrease in the equilibrium genetic variance with increasing selfing rates; however, the relationship between self-fertilization and the variables of interest depends on the strength of associations between loci, and three different regimes are observed. When the U/n ratio is low (where U is the total haploid mutation rate on selected traits) and effective recombination rates are sufficiently high, genetic associations between loci are negligible and the genetic variance, mutation load and inbreeding depression are well predicted by approximations based on single-locus models. For higher values of U/n and/or lower effective recombination, moderate genetic associations generated by epistasis tend to increase V g , L and δ, this regime being well predicted by approximations including the effects of pairwise associations between loci. For yet higher values of U/n and/or lower effective recombination, a different regime is reached under which the maintenance of coadapted gene complexes reduces V g , L and δ. Simulations indicate that the values of V g , L and δ are little affected by assumptions regarding the number of possible alleles per locus.

INTRODUCTION

Genetic diversity maintained within populations plays an important role in defining their adaptive potential (for a species to evolve, there must be heritable phenotypic variation on which selection can act). The ultimate source of this diversity is mutation, with a substantial proportion of new mutations being of a slightly deleterious nature [START_REF] Eyre-Walker | The ditribution of fitness effects of new mutations[END_REF]: hence, a corollary to the maintenance of genetic diversity is the existence of a mutation load, defined as the reduction in mean fitness of a population relative to the fitness of an optimal genotype [START_REF] Haldane | The effect of variation on fitness[END_REF]. Furthermore, the fact that most deleterious alleles are partially recessive causes inbred offspring to have a lower fitness (on average) than outbred ones, as they tend to carry higher numbers of homozygous mutations (inbreeding depression, [START_REF] Charlesworth | Inbreeding depression and its evolutionary consequences[END_REF].

By affecting the average degree of homozygosity of individuals and the efficiency of recombination between loci, the reproductive system of a species is expected to have an important influence on the effect of selection against deleterious alleles, and thus on the mutation load, inbreeding depression and level of diversity maintained within populations. One mating system that has received considerable attention is self-fertilization, a reproductive strategy occurring at various rates in an important proportion of plant and animal species [START_REF] Jarne | Animals mix it up too: the distribution of selffertilization among hermaphroditic animals[END_REF][START_REF] Goodwillie | The evolutionary enigma of mixed mating systems in plants: occurence, theoretical explanations, and empirical evidence[END_REF][START_REF] Igic | The distribution of plant mating systems: study bias against obligately outcrossing species[END_REF]. Self-fertilization, and inbreeding in general, may have different effects on genetic polymorphisms depending on the strength of selection acting on them [START_REF] Glémin | Mating systems and the efficacy of selection at the molecular level[END_REF]. When directional selection against deleterious alleles is sufficiently strong relative to drift (N e s 1), the increased homozygosity caused by inbreeding is expected to improve the efficiency of selection against those alleles (purging), reducing the mutation load and inbreeding depression [START_REF] Lande | The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models[END_REF][START_REF] Charlesworth | Inbreeding depression, genetic load, and the evolution of outcrossing rates in a multilocus system with no linkage[END_REF]. At the other extreme, polymorphism at neutrally-behaving loci (N e s 1)

should also be lowered by inbreeding, as the effective population size is reduced by identity-by-descent within loci [START_REF] Pollak | On the theory of partially inbreeding finite populations. I. Partial selfing[END_REF] and by stronger interference effects between loci -background selection, hitchhiking [START_REF] Nordborg | Structured coalescent processes on different time scales[END_REF][START_REF] Glémin | Adaptation and maladaptation in selfing in outcrossing species: new mutations versus standing variation[END_REF]Roze, 2016). In intermediate regimes (N e s ∼ 1), however, the reduction in N e due to inbreeding may cause an increased frequency of deleterious alleles (because selection is less effective), which may explain the higher π N /π S ratio observed in various selfing species compared with their outcrossing relatives [START_REF] Brandvain | Genomic identification of founding haplotypes reveals the history of the selfing species Capsella rubella[END_REF][START_REF] Burgarella | Molecular evolution of freshwater snails with contrasting mating systems[END_REF], and other references listed in Table 1 of [START_REF] Hartfield | Evolutionary genetic consequences of facultative sex and outcrossing[END_REF].

Most classical results on the effects of selfing on genetic diversity, mutation load and inbreeding depression are based on single-locus models, and thus neglect the effects of linkage disequilibria and other forms of genetic associations among loci. Previous analytical and simulation models showed that intermediate selfing rates generate correlations in homozygosity between loci, termed "identity disequilibria" [START_REF] Weir | Mixed self and random mating at two loci[END_REF][START_REF] Vitalis | Two-locus identity probabilities and identity disequilibrium in a partially selfing subdivided population[END_REF], which tend to reduce the efficiency of purging when deleterious alleles are partially or fully recessive (an effect called "selective interference" by [START_REF] Lande | High inbreeding depression, selective interference among loci, and the threshold selfing rate for purging recessive lethal mutations[END_REF]. When the number of highly recessive mutations segregating within genomes is sufficiently high, these correlations in homozygosity may entirely suppress purging unless the selfing rate exceeds a given threshold [START_REF] Lande | High inbreeding depression, selective interference among loci, and the threshold selfing rate for purging recessive lethal mutations[END_REF][START_REF] Scofield | Mitosis, stature and evolution of plant mating systems: low-Φ and high-Φ plants[END_REF][START_REF] Kelly | Mutation-selection balance in mixed mating populations[END_REF][START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]. Linkage disequilibrium corresponds to another form of association between loci that may also affect the efficiency of selection: in particular, selection may be strongly limited by Hill-Robertson effects in highly selfing populations, due to the fact that selfing reduces the efficiency of recombination between loci -recombination having no effect when it occurs in homozygous individuals (Kamran-Disfani and Agrawal, 2014; [START_REF] Hartfield | Limits to adaptation in partially selfing species[END_REF]. Epistatic interactions represent another possible source of linkage disequilibrium between selected loci. [START_REF] Charlesworth | Multilocus models of inbreeding depression with synergistic selection and partial self-fertilization[END_REF] considered a model in which epistasis between deleterious alleles is fixed and synergistic (the effects of mutations alone being smaller than when combined with others), and showed that the effect of the selfing rate on the load and inbreeding depression may be non-monotonic under this form of epistasis, with an increase in both variables above a (high) self-fertilization threshold. However, although models with fixed epistasis have lead to important insights, epistatic interactions are known to vary across pairs of loci, and this variation may have important evolutionary consequences [START_REF] Phillips | Beyond the average: the evolutionary importance of gene interactions and variability of epistatic effects[END_REF][START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF]. Interestingly, several aspects of the complexity of epistatic interactions (such as possible compensatory effects between deleterious alleles, i.e., reciprocal sign epistasis) are captured by models of stabilizing selection acting on quantitative traits, such as Fisher's geometric model [START_REF] Fisher | The genetical theory of natural selection[END_REF]. Furthermore, the distributions of epistasis generated by this type of model seem compatible with our empirical knowledge on epistasis [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF].

Only a few models have explored the effect of self-fertilization on genetic variance for quantitative traits at equilibrium between mutation and stabilizing selection.

Modeling a quantitative trait coded by additive loci, [START_REF] Wright | The genetical structure of populations[END_REF] showed that, in the absence of selection, the genetic variance for the trait is increased by a factor 1 + F (where F is the inbreeding coefficient), due to the increased homozygosity of the underlying loci. Selection will oppose this increase in variance, however, by elim-inating genotypes that are too far from the optimum (purging). Stabilizing selection is also known to generate positive linkage disequilibrium between alleles at different loci having opposite effects on the trait [START_REF] Bulmer | The effect of selection on genetic variability[END_REF], the immediate consequence of which is to reduce the genetic variance. These linkage disequilibria should also affect the efficiency of selection at each locus, and thus indirectly affect the genetic variance. Lande (1977) proposed a model of stabilizing selection acting on a single trait coded by additive loci in a partially selfing population, in which a Gaussian distribution of allelic effects is assumed to be maintained at each locus. He found that, as the increase in variance due to homozygosity is exactly compensated by the effect of purging, and the decrease in variance caused by linkage disequilibria is exactly compensated by the decreased efficiency of selection acting at each locus due to these linkage disequilibria, overall the equilibrium genetic variance is not affected by the selfing rate of the population. More recently, [START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for the evolution of selfing[END_REF] extended this model to multiple selected traits, and used a method developed by [START_REF] Kelly | Mutation-selection balance in mixed mating populations[END_REF] to take into account the effects of correlations in homozygosity across loci by splitting the population into selfing age classes (corresponding to classes of individuals having the same history of inbreeding), while assuming a Gaussian distribution of allelic effects at each locus within each class.

Numerical iterations of the model showed that above a threshold selfing rate, a different regime is reached, in which strong compensatory associations between alleles at different loci reduce the genetic variance and may generate outbreeding depression (i.e., lower fitness of outcrossed offspring relative to selfed offspring).

The hypothesis made by Lande (1977) and [START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for the evolution of selfing[END_REF] of a Gaussian distribution of allelic effects maintained at each locus (either in the whole population or in each selfing age class) has been criticized on the grounds that it implicitly assumes an unrealistically high mutation rate per locus and/or very weak fitness effects of mutations [START_REF] Turelli | Heritable genetic variation via mutation-selection balance: Lerch's zeta meets the abdominal bristle[END_REF]. [START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for the evolution of selfing[END_REF] also considered an infinitesimal model (in which traits are coded by an infinite number of loci, selection having a negligible effect on allele frequencies at each locus), and showed that a similar threshold pattern emerges, although the effect of selfing on the genetic variance and inbreeding depression above the threshold differs between the two models (in particular, outbreeding depression is not observed in the infinitesimal model). However, the effect of selfing on the genetic variance of quantitative traits under more general assumptions regarding the strength of selection at the underlying loci remains unclear.

In this paper, we introduce partial self-fertilization into previous models of stabilizing selection acting on quantitative traits coded by biallelic loci [START_REF] Latter | Natural selection for an intermediate optimum[END_REF]Bulmer, 1972;[START_REF] Barton | The maintenance of polygenic variation through a balance between mutation and stabilizing selection[END_REF]Barton, , 1989;;[START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF][START_REF] Blanckaert | Epistasis, pleiotropy and the mutation load in sexual and asexual populations[END_REF]. Charlesworth and Charlesworth (1995) had extended such models to take complete selfing into account, and found that the genetic variance under stabilizing selection should be lower under complete selfing than under random mating. The present paper generalizes these results to arbitrary selfing rates and multiple selected traits. Assuming additive effects of alleles on phenotypes (no dominance or epistasis on phenotypic traits), we develop approximations incorporating the effects of pairwise associations between selected loci, and compare these approximations with results from individualbased simulations. Our results indicate that different regimes are possible depending on the effect of genetic associations on the genetic variance, this effect increasing as the overall mutation rate U and selfing rate σ increase, and decreasing as the genome map length R, mean fitness effect of mutations s and number of selected traits n increase.

When U and σ are sufficiently low and R, s and n sufficiently high, the effect of associ-ations is negligible and the mutation load, inbreeding depression and genetic variance are well predicted by classical expressions ignoring associations. As the strength of associations increases, a second regime is entered in which the overall effect of associations is to reduce purging, thereby increasing the genetic variance, mutation load and inbreeding depression; this "weak association" regime is generally well predicted by our approximations which include the effects of pairwise associations between loci. For yet higher U , σ and/or lower R, s or n, a third regime is reached in which strong associations between loci caused by compensatory effects among mutations reduce the genetic variance, load and inbreeding depression. Although our approximations break down in this "strong association" regime, the approximation proposed by Charlesworth and Charlesworth (1995) provides accurate results under complete selfing when the mutation rate is sufficiently high and the mean fitness effect of mutations sufficiently low.

MODEL

Genotype-phenotype map. The parameters and variables of our model are summarized in Table 1. We consider a diploid population of size N with discrete generations.

Offspring are produced by self-fertilization with probability σ, and by random union of gametes with probability 1 -σ. The fitness of an organism represents its overall relative fecundity (assumed very large for all individuals), and depends on the values of n quantitative phenotypic traits under stabilizing selection. In the following we use subscripts α, β, γ... to denote phenotypic traits, while subscripts i, j, k... denote loci.

The value of trait α in a given individual is denoted z α , and can be decomposed into a genetic and an environmental component:

z α = g α + e α (1)
where the environmental component e α is sampled from a Gaussian distribution with mean zero and variance V e (the same for all traits). The genetic component g α ("breeding value") is controlled by a large number of biallelic loci with additive effects. The two alleles at each locus are denoted 0 and 1, while X M i and X P i are defined as indicator variables that equal zero if the individual carries allele 0 at locus i on its maternally (X M i ) or paternally (X P i ) inherited chromosome, while they equal 1 if allele 1 is present.

We also assume that g α = 0 in an individual homozygous for allele 0 at all loci, so that:

g α = i=1 r αi X M i + X P i (2)
where is the number of loci affecting phenotypic traits, and r αi the effect on phenotype α of changing the allelic state of one gene copy at locus i from 0 to 1 (note that r αi may be negative).

Following [START_REF] Chevin | Fisher's model and the genomics of adaptation: restricted pleiotropy, heterogenous mutation, and parallel evolution[END_REF], [START_REF] Lourenço | Complexity, pleiotropy and the fitness effect of mutations[END_REF] and [START_REF] Blanckaert | Epistasis, pleiotropy and the mutation load in sexual and asexual populations[END_REF], a parameter m measures the degree of pleiotropy of mutations: each locus affects a subset of m phenotypic traits, sampled randomly (and independently for each locus) among the n traits. Therefore, m = 1 means that each locus affects a single trait, while m = n corresponds to full pleiotropy (each locus affecting all traits), as in Fisher's geometric model [START_REF] Fisher | The genetical theory of natural selection[END_REF]. We assume that the distribution of mutational effects r αi over all loci affecting trait α has average zero and variance a 2

(the same for all traits); if locus i does not affect trait α, then r αi = 0. For simplicity, we consider a fully isotropic model with no mutational covariance between traits. Finally, u denotes the mutation rate from allele 0 to allele 1 and from allele 1 to allele 0 at each locus, while U = u is the haploid mutation rate over all loci (per generation).

From the previous definitions, and assuming that population size is sufficiently large, mean trait values are given by:

z α ≈ g α = 2 i=1 r αi p i (3)
where p i is the frequency of allele 1 at locus i. As we assume no G × E interaction, the variance in trait α is given by:

V p, α = V g, α + V e (4)
where V g, α is the variance in g α (genetic variance). In the next subsection, we show how V g, α can be expressed in terms of genetic associations within and between loci.

Genetic associations and decomposition of the genetic variance. Genetic associations are defined as in [START_REF] Kirkpatrick | General models of multilocus evolution[END_REF]. In particular, the centered variables ζ i,∅ and ζ ∅,i are defined as:

ζ i,∅ = X M i -p i , ζ ∅,i = X P i -p i . (5)
Furthermore, products of ζ i,∅ , ζ ∅,i variables are denoted:

ζ U,V = i∈U ζ i,∅ j∈V ζ ∅,j (6)
where U and V represent sets of loci. For example, for U = V = {i}, we have:

ζ i,i = X M i -p i X P i -p i (7)
while for U = {i, j} and V = {i}:

ζ ij,i = X M i -p i X M j -p j X P i -p i . (8)
Finally, genetic associations D U,V are defined as averages of ζ U,V variables over all individuals:

D U,V = E [ζ U,V ] (9)
where E stands for the average over all individuals in the population. We also define DU,V as (D U,V + D V,U ) /2, and write DU,∅ as DU (for simplicity). In particular, D i,i is a measure of excess homozygosity (due, for example, to non-random mating) at locus i (D i,i = 0 at Hardy-Weinberg equilibrium). As shown in Supplementary File S1, it can be written as D i,i = F p i q i , where F is the inbreeding coefficient (probability of identity by descent between two alleles present at the same locus in the same individual). The association Dij corresponds to the linkage disequilibrium between loci i and j (association between alleles present on the same haplotype), while Di,j is the association between alleles at loci i and j present on different haplotypes of the same individual. We will see that the association D ij,ij also appears in the computations, and can be expressed as φ ij p i q i p j q j , where φ ij is the probability of joint identity be descent at loci i and j. The quantities φ ij and F enter into the definition of the identity disequilibrium between loci i and j, given by [START_REF] Weir | Mixed self and random mating at two loci[END_REF], which will appear in some of our results.

G ij = φ ij -F 2 (
From these definitions, and using equations 2 and 3, the genetic variance for trait α can be written as:

V g, α = E (g α -g α ) 2 = E i,j r αi r αj (ζ i,∅ + ζ ∅,i ) (ζ j,∅ + ζ ∅,j ) (10)
where the last sum is over all i and j (including i = j). Using the fact that Dii = p i q i , one obtains from equation 9:

V g, α = 2 i=1 r αi 2 (p i q i + D i,i ) + 2 i, j =i r αi r αj Dij + Di,j . (11)
Following previous usage (e.g., Bulmer, 1985), we will call genic variance (denoted V 0 g, α ) the quantity 2 i=1 r αi 2 p i q i , corresponding to the genetic variance in a population with the same allele frequencies, but in the absence of genetic association (within and between loci). As shown by equation 11, excess homozygosity tends to increase the genetic variance through the term in D i,i . The second term of equation 11 (the effect of between-locus associations) tends to be negative under stabilizing selection, since the allele increasing the value of trait α at locus i tends to be associated with the allele decreasing its value at locus j (e.g., [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]. However, below we show that that excess homozygosity and associations between loci also affect equilibrium allele frequencies, and thus the genic variance.

Fitness function. Most of the results derived in this paper assume an isotropic, Gaussian fitness function, the fitness of an individual being given by:

W = exp - d 2 2ω 2 (12)
where ω 2 measures the strength of selection, and d = n α=1 z 2 α is the Euclidean distance (in phenotypic space) between the individual's phenotype and the optimum, which we assume is located at z = (0, 0, . . . , 0). From equation 12, the fitness associated with a given genotype (obtained by averaging over environmental effects) is also Gaussian, and given by:

W g = W g,max exp - d g 2 2V s (13) with V s = ω 2 + V e , W g,max = (ω 2 /V s ) n/2
(the mean fitness of an optimal genotype), and

d g = n α=1
g 2 α (the Euclidean distance between the breeding value of the individual and the optimum). Under our mutational model, the mean reduction in log W g caused by a heterozygous mutation present in an optimal genotype is:

s = ma 2 2V s (14)
(e.g., Martin and Lenormand, 2006b). Under our assumption of additivity of phenotypic effects it is easy to show that the reduction in log W g caused by a homozygous deleterious allele (in an optimal genotype) is four times the reduction caused by the same allele in the heterozygous state. Provided that most mutations have weak fitness effects (so that log (1 -s) ≈ -s), the dominance coefficient of deleterious alleles is thus close to 0.25 at the fitness optimum (see [START_REF] Manna | Fitness landscapes: an alternative theory for the dominance of mutation[END_REF] for more general results on dominance in Fisher's geometric model).

The effect of the shape of the fitness peak will be explored using a generalized version of equation 13 (e.g., Martin and Lenormand, 2006a;[START_REF] Tenaillon | Quantifying organismal complexity using a population genetic approach[END_REF]:

W g = W g,max exp - d g √ 2V s Q . ( 15
)
Gaussian fitness (equation 13) thus corresponds to Q = 2, while the fitness peak is sharper around the optimum when Q < 2, and flatter when Q > 2. Importantly, Q affects the average dominance coefficient of deleterious alleles, making them more dominant for Q < 2 and more recessive for Q > 2 [START_REF] Manna | Fitness landscapes: an alternative theory for the dominance of mutation[END_REF], as well as the average epistasis (on fitness) between alleles, positive for Q < 2, and negative for Q > 2 [START_REF] Gros | The evolution of epistasis and its links with genetic robustness, complexity and drift in a phenotypic model of adaptation[END_REF]. Approximations for the mutation load and inbreeding depression can be derived for Q = 2 as long as the distribution of breeding values in the population is approximately Gaussian.

Individual-based simulations. In order to verify the analytical results obtained, individual-based simulations were run using a C++ program described in Supplementary File S5 (and available from Dryad), in which the genome of each individual consists of two copies of a linear chromosome carrying equidistant biallelic loci affecting the n traits under selection. Another version of the program was used to consider a different genetic architecture, under which an infinite number of alleles are possible at each locus (see Supplementary File S5).

RESULTS

Neglecting associations between loci. In the following section we show that genetic associations between loci may be neglected when the haploid genomic mutation rate U is sufficiently low. In this case, equation 11 simplifies to:

V g, α ≈ 2 i=1 r αi 2 (p i q i + D i,i ) . ( 16
)
Expressions for p i q i and D i,i at equilibrium, assuming weak selection (V g, α V s ) and neglecting associations among loci are derived in Supplementary File S1. To leading order, D i,i ≈ F p i q i where F = σ/ (2 -σ) is the inbreeding coefficient. Neglecting associations between loci and assuming that mean phenotypes are at the optimum (g α = 0), the effect of selection on p i q i is given by:

∆ sel p i q i ≈ -s i (1 + 3F ) (1 -2p i ) 2 p i q i ( 17
)
where

s i = n α=1 r 2 αi / (2V s )
is the heterozygous effect of a mutation at locus i on log fitness in an optimal genotype. Furthermore, because mutation changes p i to

p i (1 -u) + u (1 -p i
), the change in p i q i due to mutation is (to the first order in u):

∆ mut p i q i ≈ u (1 -2p i ) 2 . ( 18
)
In regimes where genetic drift can be neglected, ∆ sel p i q i = -∆ mut p i q i at mutationselection balance, leading to either p i = 1/2 or:

s i p i q i ≈ u 1 + 3F , (19)
in agreement with results of previous biallelic models under random mating (e.g., Bulmer, 1972;[START_REF] Barton | The maintenance of polygenic variation through a balance between mutation and stabilizing selection[END_REF]. A stability analysis indicates that the equilibrium given by equation 19 is stable when p i q i < 1/4 (that is, when s i (1 + 3F ) > 4u), otherwise

p i = 1/2 is stable.
When all loci are at the equilibrium where p i q i < 1/4, summing both sides of equation 19 over i yields, using

s i = n α=1 r 2 αi / (2V s ): n α=1 V 0 g, α ≈ 4V s U 1 + 3F ( 20
)
where again V 0 g, α = 2 n α=1 r 2 αi p i q i is the genic variance. By symmetry the equilibrium genic variance should be the same for all traits, and thus:

V 0 g, α ≈ 4V s U n (1 + 3F ) = 2V s U n 2 -σ 1 + σ . ( 21
)
From equations 16 and 21, the equilibrium genetic variance is:

V g, α ≈ V 0 g, α (1 + F ) ≈ 4V s U n (1 + σ) . (22)
When σ = 0 and n = 1, equation 22 is equivalent to the result of previous biallelic models (e.g., [START_REF] Latter | Natural selection for an intermediate optimum[END_REF]Bulmer, 1972) and to Turelli's house-of-cards approximation [START_REF] Turelli | Heritable genetic variation via mutation-selection balance: Lerch's zeta meets the abdominal bristle[END_REF].

Assuming that the variance in log-fitness is small, mean fitness is approximately W ≈ e log Wg . Defining the mutation load L as the reduction in W relative to the average fitness of an optimal genotype, one obtains from equation 13:

L = 1 - W W g,max ≈ 1 -exp - n α=1 V g, α 2V s . (23)
Equations 22 and 23 yield:

L ≈ 1 -exp - 2U 1 + σ . ( 24
)
Inbreeding depression δ measures the mean fitness of selfed offspring, relative to the mean fitness of outcrossed offspring. Under the same assumptions, it is given by:

δ = 1 - W self W out ≈ 1 -exp - n α=1 V self g, α -V out g, α 2V s (25)
where V self g, α and V out g, α are the genetic variances for trait α among selfed and outcrossed offspring, respectively (e.g., [START_REF] Lande | The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models[END_REF]. The intralocus as-

sociation D i,i among selfed offspring is D self i,i = 1 2 (p i q i + D i,i ) and therefore V self g, α = V 0 g, α 1 + 1 2 (1 + F ) , while V out g, α = V 0 g, α , yielding (using equation 21): δ ≈ 1 -exp - n α=1 V g, α 4V s ≈ 1 -exp - U 1 + σ . (26)
Equations 24 and 26 are equivalent to the classical expressions obtained for the load and inbreeding depression at mutation-selection balance when the dominance coefficient h of deleterious alleles is set to 0.25 (e.g., [START_REF] Charlesworth | Inbreeding depression and its evolutionary consequences[END_REF], in agreement with the fact that h ≈ 0.25 under Gaussian stabilizing selection when mutations have additive effects on phenotypes (see previous section).

Figure 1 shows that the mutation load is well predicted by equation 24 when andS3 show that the genetic variance and inbreeding depression follow similar patterns. Drift may have significant effects on genetic variation, however, when N e s is ≈ 1 or lower. Following Bulmer (1972), a diffusion model can be used to compute the expected value of p i q i under selection, mutation and drift, provided that the effects of associations between loci are neglected.

As explained in Supplementary File S2, the result can then be integrated over the distribution of s i across loci to obtain the equilibrium genetic variance, inbreeding depression and mutation load. Figures 1 and S1 -S3 show that drift increases V g , L and δ in regimes where p i q i tends to stay small at most loci at the deterministic equilibrium (s = 10 -2 , 10 -3 in Figure 1), and has the opposite effect in regimes where p i q i is high (s = 10 -4 in Figure 1). Simple approximations can be obtained when the effect of selection is negligible at most loci (see Supplementary File S2), which provide accurate predictions when N e s is sufficiently low, or when s u so that p i = 1/2 at most loci at the deterministic equilibrium (Figures S1 -S4 (respectively) for N = 5,000 and σ = 1 (right-most points in Figure 1B). Replacing

N by N e (1 + F ) in the diffusion model provides predictions that closely match the simulation results for s = 10 -4 and 10 -3 , suggesting that the initial discrepancy was indeed caused by background selection reducing N e (results not shown). However, for s = 10 -2 , the diffusion model still performs poorly despite the corrected N e . This implies that the discrepancy between analytical and simulation results is more likely due to interactions among loci, and possibly also to deviations of mean phenotypes from the optimum caused by genetic drift (that are not taken into account in the analysis).

In Supplementary File S3, we derive expressions for the genetic variance, mutation load and inbreeding depression (for both the mutation-selection and the mutationdrift regimes) under the generalized fitness function given by equation 15. In the mutation-selection regime (s i 1/N e , u at most loci), one obtains:

V g V s ≈ 4U Q (1 + σ) Γ n 2 Γ Q+n 2 2 Q (27) (where Γ is Euler's Gamma function), while L ≈ 1 -exp - 4U Q (1 + σ) (28) δ ≈ 1 -exp - 4U Q (1 + σ) 3 -σ 2 Q 2 - 2 -σ 2 Q 2 , ( 29
)
these equations being equivalent to equations 22, 24 and 26 when Q = 2. As shown by Figure 2, equations 27 -29 provide good predictions of the simulation results when the population size and number of loci are sufficiently large (and selfing is not too high). As Q increases, the fitness peak becomes flatter around the optimum, and the equilibrium genetic variance increases (Figure 2B). However, despite increasing the genetic variance, higher values of Q lead to lower mutation loads due to the fact that deleterious alleles are more often eliminated when present in combination within the same genome: this corresponds to the classical result that negative epistasis reduces the mutation load in sexually reproducing populations (e.g., [START_REF] Kimura | The mutational load with epistatic gene interactions in fitness[END_REF][START_REF] Kondrashov | King's formula for the mutation load with epistasis[END_REF]. Indeed, the average epistasis between deleterious alleles equals zero for Q = 2, but becomes negative when Q > 2, and positive when Q < 2 [START_REF] Gros | The evolution of epistasis and its links with genetic robustness, complexity and drift in a phenotypic model of adaptation[END_REF]. By contrast, inbreeding depression is less affected by Q, δ slightly increasing or decreasing as Q increases, depending on the selfing rate. for m = 1, due to the higher variance in fitness effects of mutations, causing a larger fraction of loci to be substantially affected by both selection and drift (this effect being captured by the diffusion model). Similarly, the effect of m is more pronounced when s = 10 -2 and σ = 1, as N e is greatly reduced by background selection when selfing is high, causing higher proportions of loci to be substantially affected by drift.

As shown by Figure 3B, the number of selected traits n has only little effect on the load in the mutation-drift regime (s = 10 -4 ), in agreement with equation B9 in Supplementary File S2. However, while the diffusion model also predicts very little effect of n in the mutation-selection regime (s = 10 -2 ), larger effects are observed in the simulations, with larger deviations from the analytical predictions (and higher load)

for lower values of n. These deviations are caused by associations between loci (which are neglected in equation 24 and in the diffusion model). In the next subsection, we

show that the relative effect of these associations is indeed stronger when n is lower, and derive an approximation including the effect of pairwise genetic associations that better matches the simulation results.

Effects of associations between loci. In Supplementary File S1, we derive approximations for the effects of associations between pairs of loci on the genetic variance at mutation-selection balance, under a Gaussian fitness function (Q = 2). For this, we assume that these associations remain weak, and neglect the effects of all associations involving more than two loci. As shown by equation 11, associations Dij , Di,j (between alleles at loci i and j, either on the same or on different haplotypes) directly affect the genetic variance. At equilibrium, these associations are approximately given by:

Dij ≈ 1 1 -F 1 ρ ij + 2F ∆ sel Dij (30) Di,j ≈ F 1 -F 1 ρ ij + 2 ∆ sel Dij , (31)
where again F = σ/ (2 -σ), ρ ij is the recombination rate between loci i and j, and ∆ sel Dij is the change in Dij and Di,j due to selection:

∆ sel Dij ≈ - n α=1 r αi r αj V s (1 + F ) 2 + G ij p i q i p j q j . ( 32
)
The term G ij in equation 32 represents the identity disequilibrium between loci i and j: the correlation in identity by descent between loci, generating a correlation in homozygosity (Weir and Cockerham, 1973, Supplementary File S1). Equation 32shows that stabilizing selection generates a positive association between alleles at loci i and j that tends to displace phenotypes in opposite directions (allele 1 with allele 1, and allele 0 with allele 0 if n α=1 r αi r αj < 0): the effect of the deleterious allele at locus i is then partially compensated by its associated allele at locus j (e.g., 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]. This effect of selection is strengthened by homozygosity (and correlations in homozygosity between loci) caused by selfing.

As may be seen from equations 30 and 31, Di,j ≈ F Dij when loci are tightly linked (ρ ij 1), as expected from separation of timescales arguments (e.g., [START_REF] Nordborg | Structured coalescent processes on different time scales[END_REF]Roze, 2016). However, our approximations diverge as recombination tends to zero (or as the selfing rate tends to 1), due to the assumption that genetic associations remain weak.

We show in Supplementary File S1 how Dij and Di,j can be summed over all pairs of loci in order to compute their overall direct effect on the genetic variance (second term of equation 11). These associations depend on recombination rates through the terms in 1/ρ ij in equations 30 and 31, and also through the identity disequilibrium G ij in equation 32. However, because G ij only weakly depends on the recombination rate, its average over all pairs of loci is often very close to the value obtained under free recombination, provided that the genome map length is not too small (see Supplementary Figure S5). In the following, we thus approximate G ij by its expression for freely recombining loci, denoted G:

G = 4σ (1 -σ) (4 -σ) (2 -σ) 2 . ( 33
)
By contrast, linkage has more effect on the average of 1/ρ ij over all pairs of loci, corresponding to the inverse of the harmonic mean recombination rate between all pairs of loci (denoted ρ H thereafter). Assuming that the number of loci is large, one obtains for the direct effect of linkage disequilibria on the genetic variance (see Supplementary File S1):

2 i, j =i r αi r αj Dij + Di,j ≈ - 2 V s 2 + σ (1 -σ) (2 -σ) (4 -σ) 1 ρ H + 2σ V 0 g, α 2 (34)
where the genic variance V 0 g, α may be replaced by its expression to leading order, given by equation 21. Equation 34shows that the immediate effect of associations between alleles with compensatory phenotypic effects is to reduce the genetic variance (since this term is negative). The fraction in equation 34 is an increasing function of σ, which implies that self-fertilization increases the strength of associations, thus decreasing V g, α . However, because the genic variance is expected to decrease with σ (equation 21), the direct effect of linkage disequilibria on V g, α may remain approximately constant (or even slightly decrease) as σ increases from zero.

Associations between loci do not only affect V g, α through equation 34, however, but also affect the equilibrium allele frequencies and the excess homozygosity D i,i at each locus. The effect on D i,i is mainly driven by identity disequilibria: indeed, neglecting associations between 3 or more loci, one obtains (see Supplementary File S1):

D i,i ≈ F 1 -2 j =i s j G ij p j q j p i q i . ( 35
)
Equation 35 is equivalent to equation 5 in [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF] (which is expressed to the first order in p j ), as can be noted by replacing s and h in Roze ( 2015) by 4s j and 1/4. It shows that identity disequilibria reduce the excess homozygosity at each locus: this is due to the fact that homozygotes at locus i are more likely to be also homozygous at locus j, while homozygotes at locus j have a lower fitness than heterozygotes when deleterious alleles are partially recessive. Identity disequilibria thus tend to reduce the genetic variance through this effect on D i,i , by an amount corresponding to the sum of the term in G ij in equation 35 over all pairs of loci. Approximating G ij by its expression for freely recombining loci (equation 33), one obtains that this effect reduces

V g, α by approximately -n F G V 0 g, α
2 / (2V s ), where again V 0 g, α may be replaced by the expression given by equation 21, to leading order (see Supplementary File S1).

Finally, associations between loci affect equilibrium allele frequencies (p i q i ) at each locus. As shown in Supplementary File S1, both the linkage disequilibria generated by epistasis and the identity disequilibria caused by partial selfing reduce the efficiency of purging, thereby increasing p i q i and thus the genic variance. Indeed, an expression for the effect of selection on p i q i that includes the effects of pairwise associations is, to leading order:

∆ sel p i q i ≈ -s i (1 -2p i ) 2 p i q i + 3D i,i -6 (1 + F ) j =i s j G ij p i q i p j q j + 2 (1 -2p i ) 2 j =i a ij (1 + 2F ) Dij + Di,j + a ij 16 (1 + σ) (2 + σ) (2 -σ) 2 (4 -σ) p i q i p j q j ( 36
)
with a ij =α r αi r αj / (2V s ), and where Dij , Di,j and D i,i are given by equations 30, 31 and 35. The first line of equation 36 is equivalent to equation 6 in [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF],

showing that identity disequilibria reduce the efficiency of purging by decreasing the excess homozygosity (D i,i ), and by two additional effects represented by the term in 6 (1 + F ) (see [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF] for interpretation of these effects). The term on the second and third lines (proportional to a 2 ij ) represents the effect of epistasis between loci: this term also reduces purging, since selection against deleterious alleles is less efficient when these alleles are partially compensated by alleles at other loci.

An expression for the genic variance at mutation-selection balance is given by equation A65 in Supplementary File S1. From this, one obtains for the genetic variance:

V g, α = 4V s U n (1 + σ) 1 + 2U σ (1 -σ) [6 + σ (2 -σ)] (2 -σ) (4 -σ) (1 + σ) 2 + 4U n (2 -σ) (2 + σ) (1 -σ) (4 -σ) (1 + σ) 1 2ρ H + 2 + σ (1 -σ) (2 -σ) 2 + σ (1 -σ) (37)
where the terms in U between the brackets correspond to the effect of between-locus associations. The first of these terms (on the first line of equation 37) represents the effect of identity disequilibria, while the term in U/n on the second line represents the effect of epistasis (compensatory effects between alleles at different loci). Both terms are positive, indicating that the overall effect of interactions between loci is to increase the genetic variance, due to the fact that correlations in homozygosity and compensatory effects between mutations both reduce the efficiency of purging (equation 36). Furthermore, while the effect of identity disequilibria scales with U , the effect of epistasis scales with U/n: indeed, it becomes less and less likely that alleles at different loci have compensatory effects on all of the traits as the dimensionality of the fitness landscape increases. Finally, the effect of epistasis is more strongly affected by linkage between loci (through the term in 1/ρ H ); the effect of linkage on the term in U representing the effect of linkage disequilibria is weaker, and has been neglected in equation 37. Under random mating (σ = 0), equation 37 simplifies to:

V g, α = 4V s U n 1 + 2U n 1 ρ H + 2 (38)
which takes a similar form as equation 4.16 in Turelli and [START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF] in the case of a single selected trait (n = 1). Supplementary Figure S6 shows how the equilibrium genetic variance and its different components vary with the selfing rate, in a regime where both identity disequilibria and epistasis have significant effects (U = 1).

Assuming that W ≈ e log Wg , an approximation for the load at mutation-selection balance is 1 -exp [-n V g / (2V s )], where V g is the genetic variance given by equation 37 (the same for all traits). A slightly better approximation can be obtained by using

W ≈ e log Wg 1 + 1 2 Var [log W g ] ,
where Var [log W g ] is the variance in log fitness in the population [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]. To leading order, it is given by (see Supplementary File S1):

Var [log W g ] ≈ 2s U 1 + 3σ 1 + σ + 4U 2 σ (1 -σ) (4 -σ) (1 + σ) 2 + 8U 2 n (2 -σ) (2 + σ) (4 -σ) (1 + σ) 2 , ( 39
)
simplifying to 2U (s + 4U/n) in the absence of selfing. The first term of equation 39 represents the sum of single-locus contributions to the variance in log fitness, while the second and third term correspond to the effects of identity disequilibria and epistasis (respectively), both increasing the variance in fitness. The mutation load is then given by:

L ≈ 1 -1 + 1 2 Var [log W g ] exp - nV g 2V s . ( 40
)
Similarly, we show in Supplementary File S1 that an expression for inbreeding depression including the effects of pairwise associations between loci is:

δ ≈ 1 -1 + 1 2 ∆Var [log W g ] exp - nV g 4V s (41)
where ∆Var [log W g ] is the difference in variance in log fitness between selfed and outcrossed offspring, given by:

∆Var [log W g ] ≈ 7sU 1 + σ + U 2 σ (1 -σ) (4 -σ) (1 + σ) 2 + 2U 2 n (10 -σ) (2 -σ) (4 -σ) (1 + σ) 2 (42)
(the terms in Var [log W g ] and ∆Var [log W g ] in equations 40 and 41 are often small, however, and may thus be neglected). After replacing V g , Var [log W g ] and ∆Var [log W g ] by the expressions given by equations 37, 39 and 42, the approximations obtained for the load and inbreeding depression include terms in U 2 representing the effect of identity disequilibria, and terms in U 2 /n representing the effect of epistasis between loci.

The terms in U 2 are identical to the terms representing the effect of identity disequilibria in a model with purely multiplicative selection against deleterious alleles (no epistasis) when setting the dominance coefficient h of deleterious alleles to 1/4 (equations 11 and 14 in [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]. The novelty here thus corresponds to the effect of epistasis (compensatory effects between deleterious alleles), that tends to increase V g , L, δ by reducing the efficiency of purging.

Figure 3B shows that equations 37, 39 and 40 capture the increase in load observed in the simulations as the number of traits n decreases (see Supplementary Figure S7 for the genetic variance and inbreeding depression). Note that the harmonic mean recombination rate ρ H between pairs of loci under our simulated genetic architecture (linear chromosome with equally spaced loci) can be obtained from: [START_REF] Blanckaert | Epistasis, pleiotropy and the mutation load in sexual and asexual populations[END_REF], yielding ρ H ≈ 0.42 for = 1,000 and R = 20. Figure 4 shows that for low or moderate selfing rates, decreasing the genome map length from R = 20 to R = 1 increases the mutation load, by increasing the strength of linkage disequilibria caused by epistasis, that in turn reduce the efficiency of purging. In this regime, equations 37, 39 and 40 provide an accurate prediction for the load (see Supplementary Figure S8 for genetic variance and inbreeding depression). At high selfing rates, however, a different regime is entered, in which the assumption of weak genetic associations breaks down. As can be seen in Figure 4, in this regime (which spans a broader parameter range under tighter linkage) the load decreases more rapidly as σ increases. Increasing linkage tends to reduce the mutation load when the selfing rate is high, although the effect of R vanishes when σ = 1. When linkage is extremely tight, the approximations given above break down for all values of σ: as shown by Figures 4 andS8, decreasing R has a non-monotonic effect on the genetic variance, load and inbreeding depression when selfing is small to moderate, the lowest values of V g , L and δ being reached when R = 0 (in which case selfing has no effect).

1 ρ H = 2 ( -1) -1 i=1 2 ( -i) 1 -exp -2i R -1 (43) (see Appendix 2 in
An approximation for the genetic variance under complete linkage can be obtained by treating the whole genome as a single locus with a very large number of possible alleles, and assuming a Gaussian distribution of allelic effects in the population (Lande, 1977;Supplementary File S4). This yields:

V g 2V s ≈ 2 sU 2n 2 2 + σ (1 -σ) , (44)
L ≈ 1 -exp - nV g 2V s , δ ≈ 1 -exp - nV g 4V s . ( 45
)
Note that equation 44 is equivalent to equation 3A in Charlesworth and Charlesworth (1995) when σ = 1 and n = 1. As shown by Figures 4 andS8, equations 44 and 45 only slightly overestimate V g , L and δ when σ = 1 and/or R = 0. As shown below, better predictions are observed for higher values of U/n and lower values of s.

The effects of identity disequilibria between loci remain negligible for the parameter values used in Figures 3 and4. As shown by Figure 5, identity disequilibria become more important for higher values of the mutation rate U . Indeed, the relative effects of identity disequilibria on the load can be deduced from the differences between the three curves in each panel of Figure 5, the red curves showing the predicted mutation load in the absence of epistasis, but taking into account identity disequilibria (obtained by removing the terms in U 2 /n from equations 37, 39 and 40, leading to an expression equivalent to equation 11 in [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]. The difference between the black and red curves thus represents the predicted effect of identity disequilibria on the load, while the difference between the red and green curves corresponds to the additional effect of epistasis. Simulations indicate that the change in regime observed above a threshold selfing rate (around σ = 0.5 for U = 1 in Figure 5) is due to epistasis, since this threshold is not observed in simulations without epistasis (red dots). Supplementary Figure S9 shows that this threshold pattern is little affected by population size N , as long as the effects of drift remain small. Similarly, the results only weakly depend on the number of loci , as long as the mutation rate per locus u = U/ is small enough so that p i q i < 1/4 at most loci (see Supplementary Figure S10 for distributions of allele frequencies in simulations with = 1,000 and = 10,000). Figure S9 also shows that the results are little affected by the degree of pleiotropy of mutations m, as long as s remains constant. However, s does affect V g , L and δ in the regime where our approximations break down. As shown by Figure 6, decreasing s lowers the threshold selfing rate above which our approximations are not valid and results in lower equilibrium mutation loads (see Supplementary Figure S11 for results on V g and δ). Figures 6 and S11 also show that, when s is sufficiently small, the single-locus, Gaussian model (equations 44 and 45, dotted curves on the figures) provides accurate predictions for V g , L and δ under complete selfing (σ = 1).

In Figure 7 we show that decreasing the number of traits under selection n decreases the threshold selfing rate above which our approximations break down (see Supplementary Figure S12 for inbreeding depression and scaled genetic variance). Below the threshold, the mutation load decreases as n increases, as predicted by our analytical results (although our approximations become less precise for low n and high U ), while n has the opposite effect above the threshold. Overall, we observe that in this second regime (in which interactions between loci have important effects), the mutation load generally increases with the number of selected traits, the fitness effects of mutations s, the mutation rate U and recombination rate (through the parameter R). However, Figure 8 shows that the effects of these parameters on inbreeding depression are more complicated. In particular, outbreeding depression (negative δ) may occur in regimes where the effects of epistasis are particularly strong (high U , low n)

and when the selfing rate is moderate to high (above 0.5 but below 1), outbreeding depression becoming stronger when s, U and R increase (the approximation derived from equation 37 fails for all values of σ for the parameter values used in Figure 8, and is not shown here). Supplementary Figures S13 andS14 show that for the same parameter values, V g and L always increase when s, U and R increase.

DISCUSSION

The response of a population to environmental change depends critically on its genetic diversity. Our results predict that the level of genetic variation maintained at equilibrium under stabilizing selection acting on quantitative traits is generally lower in more highly selfing population, due to more efficient purging (although increasing selfing may sometimes increase genetic variation, for example when mutations have weak fitness effects, as shown by Figure S3). This finding agrees with Charlesworth and Charlesworth's (1995) theoretical prediction that fully selfing populations should maintain lower genetic variance for quantitative traits under stabilizing selection than fully outcrossing ones, and with several empirical studies comparing levels of genetic variation for morphological traits in closely related pairs of plant species with contrasted mating systems (Charlesworth and Charlesworth, 1995;[START_REF] Geber | Inheritance and natural selection on functional traits[END_REF]Bartkowska and Johnston, 2009 and references therein). We also show that the lower level of variation present in more highly selfing populations is associated with lower values of the mutation load and inbreeding depression. The meta-analysis carried out by [START_REF] Winn | Analysis of inbreeding depression in mixed-mating plants provides evidence for selective interference and stable mating systems[END_REF] showed that inbreeding depression is indeed lower in highly selfing plant species compared to species with lower selfing rates, while no significant difference is observed between species with low vs. intermediate selfing rates. It has been put forth that correlations in homozygosity between selected loci may suppress purging at moderate selfing rates ("selective interference", [START_REF] Lande | High inbreeding depression, selective interference among loci, and the threshold selfing rate for purging recessive lethal mutations[END_REF][START_REF] Winn | Analysis of inbreeding depression in mixed-mating plants provides evidence for selective interference and stable mating systems[END_REF]; this, however, would imply that a large number of segregating deleterious alleles have very low dominance coefficients, generating very high inbreeding depression [START_REF] Kelly | Mutation-selection balance in mixed mating populations[END_REF][START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF], which seems unlikely. Another possible explanation for the lack of purging at intermediate selfing rates involves epistasis (compensatory effects between mutations coding for the same quantitative trait, [START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for the evolution of selfing[END_REF].

Our analysis of the effects of epistasis (under assumptions that differ from those made in Lande and Porcher's model) shows that different regimes are possible, and outlines how the parameters affect transitions between these regimes.

In our model, the effect of epistasis on the equilibrium genetic variance V g is inversely proportional to effective recombination rates between selected loci, and scales with U/n (where n is the number of selected traits and U the total mutation rate on those traits). Indeed, U/n determines the number of segregating "interacting" mutations, that is, mutations with epistatic fitness effects. As n tends to infinity, all mutations become orthogonal in phenotypic space (with independent fitness effects), and our results converge to the results from previous population genetics models without epistasis (e.g., [START_REF] Charlesworth | Inbreeding depression and its evolutionary consequences[END_REF][START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]. When U/n is small and map length R is sufficiently large, associations between loci have little effect. Under Gaussian stabilizing selection (Q = 2), the average coefficient of epistasis between mutations (on fitness) is zero (Martin and Lenormand, 2006b) while the dominance coefficient of deleterious alleles (in an optimal genotype) is close to 0.25 under the assumption of additive effects on phenotypes. In this case, we found that classical deterministic expressions based on single-locus models (hence neglecting the variance in epistatic interactions) provide accurate predictions for the mutation load L and inbreeding depression δ. Simple approximations are also obtained under the more general fitness function given by equation 15, confirming that the mutation load is an increasing function of the average coefficient of epistasis between mutations [START_REF] Kimura | The mutational load with epistatic gene interactions in fitness[END_REF][START_REF] Kondrashov | King's formula for the mutation load with epistasis[END_REF][START_REF] Phillips | Beyond the average: the evolutionary importance of gene interactions and variability of epistatic effects[END_REF][START_REF] Blanckaert | Epistasis, pleiotropy and the mutation load in sexual and asexual populations[END_REF]. Neglecting the effect of associations between loci also allowed us to explore the effects of drift using diffusion methods. As in previous studies (e.g., [START_REF] Charlesworth | Why we are not dead one hundred times over[END_REF][START_REF] Blanckaert | Epistasis, pleiotropy and the mutation load in sexual and asexual populations[END_REF], we found that drift may lower the mutation load by reducing V g . However, this result probably strongly depends on the assumption that mutations may increase or decrease phenotypic traits with the same probability (no mutational bias): indeed, previous works showed that drift may increase the load in the presence of a mutational bias by displacing mean phenotypes away from the optimum [START_REF] Zhang | The anomalous effects of biased mutation revisited: mean-optimum deviation and apparent directional selection under stabilizing selection[END_REF][START_REF] Charlesworth | Why we are not dead one hundred times over[END_REF]. Given that partial selfing reduces effective population size, it would be of interest to study the combined effects of drift and mutational bias in models with selfing.

The variance in epistasis has stronger effects as the U/n ratio increases and/or as the effective recombination rate decreases (i.e. due to selfing). Our results showed that two different regimes are possible. (1) When genetic associations (linkage disequilibria) generated by epistasis stay moderate, the overall effect of epistasis is to increase V g , L and δ by decreasing the efficiency of selection against deleterious alleles. This regime is generally well described by our model taking into account the effects of associations between pairs of loci. This result bears some similarity with the result obtained by [START_REF] Phillips | Beyond the average: the evolutionary importance of gene interactions and variability of epistatic effects[END_REF], showing that the variance in epistasis between deleterious alleles increases the mutation load. Equation 2.1 in [START_REF] Phillips | Beyond the average: the evolutionary importance of gene interactions and variability of epistatic effects[END_REF] is not fully equivalent to our expression for the load under random mating, however, possibly due to different assumptions on the relative orders of magnitude of s i , s j and e ij (where s i and s j are the strength of selection at loci i and j and e ij is epistasis between those loci) and how they covary. Nevertheless, it is interesting to note that both results become equivalent if Var [e ij / (s i s j )] in [START_REF] Phillips | Beyond the average: the evolutionary importance of gene interactions and variability of epistatic effects[END_REF] is replaced by Var [e ij /s 2 ] (where Var stands for the variance across all pairs of loci), using the fact that Var [e ij ] = 4s 2 /n in our model (with a Gaussian fitness function). ( 2) Increasing the value of U/n and/or reducing effective recombination rates or s generates a transition to a different regime in which the effect of the variance in epistasis switches, reducing V g , L and δ.

Because our analytical approach fails in this regime (presumably due to higher-order associations between loci), it is more difficult to obtain an intuitive understanding of the selective mechanisms involved. However, it is likely that selection operates on multilocus genotypes (comprising combinations of alleles with compensatory effects) that can be maintained over many generations due to high selfing rates and/or low recombination. A similar transition from genic to genotypic selection as recombination decreases was described by [START_REF] Neher | Competition between recombination and epistasis can cause a transition from allele to genotype selection[END_REF], using a haploid model in which epistasis is randomly assigned to genotypes.

Although our results show some qualitative similarities with those obtained by [START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for the evolution of selfing[END_REF] -e.g., the same transition between regimes occurs in both models as selfing increases -several differences can be observed. In particular, Lande and Porcher's model predict little or no effect of selfing on V g below the threshold selfing rate corresponding to the change in regime, and an abrupt change in V g at the threshold (except in their infinitesimal model). A step change such as this is never observed in our model, even for parameter values at which the effect of drift should be negligible at most loci. These differences between the models are not due to the different genetic architectures considered (biallelic vs. multiallelic): indeed, Supplementary Figures S15 andS16 show that assuming biallelic loci or an infinite number of possible alleles per locus in our individual-based simulations yields very similar results (for = 1,000 and = 10,000). Rather, they must be due to Lande and Porcher's assumption of a Gaussian distribution of allelic effects maintained at each locus in each selfing age class, implicitly assuming a sufficiently high mutation rate per locus u and low fitness effect of mutations s [START_REF] Turelli | Heritable genetic variation via mutation-selection balance: Lerch's zeta meets the abdominal bristle[END_REF]. In our multiallelic simulations (with u = 10 -5 to 10 -3 and s = 0.01), the number of alleles maintained at each locus is not sufficiently large to generate a Gaussian distribution of segregating allelic effects (see Figures S15 andS16). One may also note that the effect of the number of selected traits n seems different in both models (compare Lande and Porcher's Figure 5 and 6 to our Figure 7), but this is due to the fact that the overall mutation rate U is proportional to n in Lande and Porcher's model (while U is fixed in Figure 7). Increasing both n and U in order to maintain a constant U/n ratio, we indeed observed that the transition between regimes occurs at lower selfing rates when n is larger, as in Lande and Porcher's Figure the degree of pleiotropy of mutations (Lande and Porcher assume no pleiotropy). Our model allowed us to explore the effects of pleiotropy through the parameter m, showing that pleiotropy mostly affects the results through its effect on s (equation 14). The equilibrium genetic variance thus depends on m in regimes where V g is affected by s, in particular when N e s ≈ 1 or lower (Figures 1, S1 -S4), and when genetic associations are strong (Figure 6). However, pleiotropy may have stronger effects under different assumptions regarding the genetic architecture of traits, for example when different sets of traits are affected by different sets of loci (modular pleiotropy, [START_REF] Welch | Modularity and the cost of complexity[END_REF]. The effects of selective or mutational covariance among traits would also be interesting to explore: indeed, such covariances decrease the effective number of selected traits (Martin and Lenormand, 2006b), potentially increasing the importance of associations between loci.

In the regime where genetic associations generated by epistasis reduce V g (regime (2) mentioned above), outbreeding depression may occur due to the lower fitness of recombinants between selfing lineages maintaining coadapted gene complexes (Figure 8), a result shared with [START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for the evolution of selfing[END_REF] Gaussian Allele Model. In our additive model of phenotypic effects, outbreeding depression should only be expressed in F2 individuals (that is, among the offspring of an individual produced by a cross between different selfing lineages), once recombination has disrupted compensatory associations between alleles at different loci. This explains why outbreeding depression is not observed under complete (or nearly complete) selfing in Figure 8, as all outcrossed individuals are F1 hybrids between selfing lineages. Outbreeding depression between lineages collected from the same geographical location has been observed in highly selfing plants [START_REF] Parker | Outbreeding depression in a selfing annual[END_REF][START_REF] Volis | Epistasis in natural populations of a predominantly selfing plant[END_REF] and Caenorhabditis nematodes [START_REF] Dolgin | Inbreeding and outbreeding depression in Caenorhabditis nematodes[END_REF][START_REF] Gimond | Outbreeding depression with low genetic variation in selfing Caenorhabditis nematodes[END_REF]. In all cases, estimated selfing rates are higher than those leading to δ < 0 in our simulations, however, and outbreeding depression was observed in F1 offspring of crosses between inbred lines of nematodes. The occurrence of outbreeding depression at higher selfing rates may be partly explained by the fact that experimental crosses were often performed between genetically different lines; by contrast, in our simulations the parents of an outcrossed individual may share the same genotype (in particular when the number of genetically different selfing lineages is reduced due to the low effective size of highly selfing populations), reducing the magnitude of outbreeding depression. However, the occurrence of outbreeding depression in F1 individuals must involve dominance effects which are absent from our model.

Exploring the effects of dominance/recessivity of mutations on phenotypic traits would be an interesting extension of this work.

Due to the lower genetic diversity of self-fertilizing populations, it has been suggested that they should be less able to adapt to a changing environment (e.g., [START_REF] Stebbins | Self fertilization and population variability in higher plants[END_REF][START_REF] Williams | Natural Selection: Domains, Levels, and Challenges[END_REF][START_REF] Takebayashi | Is self-fertilization an evolutionary dead end? Revisiting an old hypothesis with genetic theories and a macroevolutionary approach[END_REF]. In the absence of epistasis, existing models indeed predict that selfing populations should have lower rates of adaptation than outcrossing ones [START_REF] Glémin | Adaptation and maladaptation in selfing in outcrossing species: new mutations versus standing variation[END_REF][START_REF] Hartfield | Limits to adaptation in partially selfing species[END_REF]. When compensatory effects between mutations are possible, however, a substantial amount of genetic variance may be hidden by genetic associations between loci in highly selfing populations [START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for the evolution of selfing[END_REF], the present study). After a change in environment, this variance may be liberated by rare outcrossing events, increasing the short-term evolutionary response of highly (but not fully) selfing populations. Exploring how selfing affects adaptation under directional selection, and more generally how the variability of epistatic interactions between loci may influence the evolution of mating systems represents a natural next step of this work.

� (�) � = ������ � = �� -� � = �� -� � = �� -� ��� ��� ��� ��� ��� σ ���� ���� ���� ���� � (�) � = �����
■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ��� ��� ��� ��� ��� σ ���� ���� ���� ���� � A � = � � = � � = �� ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■

N

  e s is sufficiently large (for U = 0.1 and n = 50), and generally decreases as selfing increases -results for different numbers of loci are shown in Supplementary Figure S1, while Supplementary Figures S2

  ). In this mutation-drift regime, V g , L and δ are nearly independent of σ when N e u 1 (the increase in variance caused by excess homozygosity being exactly compensated by the reduction in variance caused by the lower effective population size), or increase with σ for larger values of N e u. The discrepancies between analytical and simulation results observed in Figure 1 at high selfing rates are partly due to the reduction in effective population size N e caused by background selection, which is not accounted for in the diffusion model. An estimation of N e using the equilibrium diversity at a neutral locus (with an infinite number of possible alleles) at the mid-point of the chromosome (as in Roze, 2016) yielded an N e of approximately 740, 300 and 200 for s = 10 -4 , 10 -3 and 10 -2

Figure 3

 3 Figure3shows the effects of the parameters m and n (for Q = 2). The degree

Figure 1 .Figure 2 .

 12 Figure 1. Mutation load L as a function of the selfing rate σ. Black curve: approx-

Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 3456 Figure 3. Mutation load L as a function of the selfing rate σ for s = 10 -2 (top, filled

Figure 7 .

 7 Figure 7. Mutation load L as a function of the selfing rate σ, for different values

Figure 8 .

 8 Figure 8. Inbreeding depression δ as a function of the selfing rate σ, for n = m = 5

Table 1 :

 1 Parameters and variables of the model. Genic variance for trait α (2 i r 2 αi p i q i )

	N	Population size
	σ	Selfing rate
	n	Number of selected traits
	m	Degree of pleiotropy of mutations
	a 2	Variance of mutational effects on selected traits
	V e	Environmental variance (on selected traits)
	ω 2	Strength of stabilizing selection on phenotypic traits
	V s = ω 2 + V e	Strength of stabilizing selection on breeding values g α
	Q	Shape of the fitness peak (equation 15)
		Number of loci affecting selected traits
	u	Mutation rate per locus
		Mutation rate (per haploid genome) on loci affecting
	U = u
		selected traits
	R	Genome map length
		Harmonic mean recombination rate between pairs of loci
	ρ H
		affecting selected traits
		Average heterozygous effect of mutations on log fitness (in
	s
		an optimal genotype)
	z α	Value of phenotypic trait α (in a given individual)
	g α , e α	Genetic and environmental components of trait α

Acknowledgements. We thank Aneil Agrawal, Patrice David, Sylvain Glémin, Frédéric Guillaume, Thomas Lenormand, Emmanuelle Porcher, Ophélie Ronce, Joëlle Ronfort and an anonymous reviewer for helpful discussions and comments, and the bioinformatics and computing service of Roscoff's Biological Station (Abims platform) for computing time. This work was supported by the French Agence Nationale de la Recherche (project SEAD, ANR-13-ADAP-0011 and project SexChange, ANR-14-CE02-0001).