Eloïse Vanhoenacker

Linnéa Sandell

Denis Roze
email: roze@sb-roscoff.fr

Stabilizing selection, mutational bias and the evolution of sex

Keywords: adaptive landscape, epistasis, evolutionary quantitative genetics, multilocus population genetics, reproductive systems

Stabilizing selection around a fixed phenotypic optimum is expected to disfavor sexual reproduction, since asexually reproducing organisms can maintain a higher fitness at equilibrium, while sex disrupts combinations of compensatory mutations. This conclusion rests on the assumption that mutational effects on phenotypic traits are unbiased, that is, mutation does not tend to push phenotypes in any particular direction. In this paper, we consider a model of stabilizing selection acting on an arbitrary number of polygenic traits coded by bialellic loci, and show that mutational bias may greatly reduce the mean fitness of asexual populations compared with sexual ones in regimes where mutations have weak to moderate fitness effects. Indeed, mutation and drift tend to push the population mean phenotype away from the optimum, this effect being enhanced by the low effective population size of asexual populations. In a second part, we present results from individual-based simulations showing that positive rates of sex are favored when mutational bias is present, while the population evolves towards complete asexuality in the absence of bias. We also present analytical (QLE) approximations for the selective forces acting on sex in terms of the effect of sex on the mean and variance in fitness among offspring.

INTRODUCTION

Various possible evolutionary benefits of sexual reproduction have been proposed in order to explain the widespread occurrence of this reproductive mode among eukaryotes (e.g., Agrawal, 2006;Otto, 2009;[START_REF] Hartfield | Current hypotheses regarding the evolution of sex and recombination[END_REF]. These broadly fall into two categories: direct selective advantages of meiotic recombination, in particular in terms of DNA repair (e.g., [START_REF] Bernstein | Genetic damage, mutation, and the evolution of sex[END_REF][START_REF] Bernstein | Is meiotic recombination an adaptation for repairing DNA, producing genetic variation[END_REF], or indirect benefits stemming from the disruption of linkage disequilibria and other forms of genetic associations through recombination and segregation. Breaking genetic associations affects the mean fitness of offspring when the fitness effect of alleles depends on the genetic background (dominance, epistasis); it may also affect the variance in fitness among offspring, and thus the response to selection. In the absence of dominance or epistasis and under random mating, stochastic events occurring in finite populations tend to generate negative genetic associations -negative linkage disequilibrium between selected loci [START_REF] Hill | The effect of linkage on limits to artificial selection[END_REF][START_REF] Felsenstein | The evolutionary advantage of recombination[END_REF] and excess heterozygosity in diploids (e.g., [START_REF] Balloux | The population genetics of clonal and partially clonal diploids[END_REF]. Breaking these negative associations increases the variance in fitness among offspring and the efficiency of natural selection, favouring higher rates of sex or recombination [START_REF] Otto | The advantages of segregation and the evolution of sex[END_REF]Barton, 1997, 2001;[START_REF] Barton | Evolution of recombination due to random drift[END_REF][START_REF] Barton | The Hill-Robertson effect and the evolution of recombination[END_REF]Martin et al., 2006;[START_REF] Michod | Deleterious mutations and selection for sex in finite, diploid populations[END_REF].

Multilocus simulation programs showed that selection for recombination generated by such stochastic effects may be strong when sex is rare, but decreases rapidly as the baseline rate of sex in the population increases [START_REF] Keightley | Interference among deleterious mutations favours sex and recombination in finite populations[END_REF][START_REF] Hartfield | The role of advantageous mutations in enhancing the evolution of a recombination modifier[END_REF]Roze, 2014).

Genetic associations may also be produced by deterministic forces: in particular, dominance and epistatic interactions between alleles affecting fitness are known to be widespread (e.g., [START_REF] De Visser | The evolution of sex: empirical insights into the roles of epistasis and drift[END_REF][START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF][START_REF] Halligan | Spontaneous mutation accumulation studies in evolutionary genetics[END_REF][START_REF] Manna | The high-throughput yeast deletion fitness data and the theories of dominance[END_REF], and represent another source of linkage disequilibria or deviations from Hardy-Weinberg equilibrium. In randomly mating populations living in a constant environment, breaking associations generated by dominance or epistasis decreases the mean fitness of offspring (segregation or recombination load), generating a short-term cost for sex and recombination [START_REF] Barton | A general model for the evolution of recombination[END_REF]Charlesworth and Barton, 1996;[START_REF] Otto | The advantages of segregation and the evolution of sex[END_REF] -this short-term cost may turn into a short-term benefit when mating is non-random or when the environment changes in space or time [START_REF] Lenormand | The evolution of recombination in a heterogeneous environment[END_REF][START_REF] Otto | The advantages of segregation and the evolution of sex[END_REF][START_REF] Lenormand | Self-fertilization and the evolution of recombination[END_REF][START_REF] Gandon | The evolution of sex and recombination in response to abiotic or coevolutionary fluctuations in epistasis[END_REF]Agrawal, 2009). In a longer term, sex is generally beneficial when interactions cause a negative curvature of the fitness function (e.g., negative epistasis, partially recessive deleterious alleles), generating negative genetic associations that limit the efficiency of selection [START_REF] Barton | A general model for the evolution of recombination[END_REF][START_REF] Otto | The advantages of segregation and the evolution of sex[END_REF]. Multilocus simulations including fixed epistasis between loci have suggested that epistatic interactions may only play a secondary role in the evolution of recombination, however, stochastic (Hill-Robertson) effects being often stronger (Otto and Barton, 2001;[START_REF] Keightley | Interference among deleterious mutations favours sex and recombination in finite populations[END_REF]. Nevertheless, epistatic interactions are known to vary across pairs of loci (e.g., [START_REF] Phillips | Beyond the average: the evolutionary importance of gene interactions and variability of epistatic effects[END_REF][START_REF] De Visser | The evolution of sex: empirical insights into the roles of epistasis and drift[END_REF][START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF], and this variation (which should generally disfavor recombination, [START_REF] Otto | Deleterious mutations, variable epistatic interactions, and the evolution of recombination[END_REF] has not been considered in recent multilocus simulation studies on the evolution of sex and recombination.

Models of stabilizing selection acting on quantitative phenotypic traits represent a simple way of introducing distributions on epistatic interactions (on fitness), including possible compensatory effects between mutations (indeed, a mutation displacing a phenotypic trait away from the optimum can be compensated by another mutation having the opposite effect on the trait). Interestingly, the predicted distribution of fitness effects and epistatic interactions among mutations obtained from classical models such as Fisher's geometric model of adaptation with a Gaussian shaped fitness function have been shown to accurately describe empirical distributions of epistasis in bacteria and viruses [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF], justifying the use of such models to explore the effects of the variance in epistasis. Selection for recombination under stabilizing, directional or fluctuating selection acting on one or several polygenic traits has been explored by previous simulation models ( [START_REF] Smith | Selection for recombination in a polygenic model[END_REF], 1988;[START_REF] Kondrashov | Evolution of amphimixis and recombination under fluctuating selection in one and many traits[END_REF]. They showed that while recombination is disfavored under stabilizing selection around a fixed optimum, environmental change may favor recombination. A mathematical analysis based on the infinitesimal model was proposed by Charlesworth (1993) (see also Appendix 2 in [START_REF] Barton | A general model for the evolution of recombination[END_REF], showing that recombination increases the phenotypic variance by breaking negative genetic associations generated by epistatic interactions among loci, thereby increasing the speed of adaptation.

As in most evolutionary quantitative genetics models, the studies just mentioned assume unbiased mutational effects on phenotypic traits: mutations are always as likely to increase as to decrease the value of a given trait. Several authors explored the effect of mutational bias on quantitative traits (e.g., [START_REF] Waxman | The anomalous effects of biased mutation[END_REF][START_REF] Zhang | The anomalous effects of biased mutation revisited: mean-optimum deviation and apparent directional selection under stabilizing selection[END_REF]Charlesworth, 2013a,b), and showed that such a bias may significantly reduce the mean fitness of populations in regimes where drift has substantial effects at loci coding for the traits, by displacing mean phenotypes away from their optimal values (thereby introducing a component of directional selection). Although the effect of mutational bias has only been explored in sexual populations, it should in principle be stronger in asexual populations, due to their reduced effective population size caused by interference effects between loci. This may generate selection for sex and recombination in the absence of environmental change.

In this paper, we explore the effect of mutational bias in a simple, isotropic model of stabilizing selection acting on an arbitrary number of phenotypic traits, in a haploid, facultatively sexual population. We first assume a fixed rate of sex in the population, and show that mutational bias may strongly reduce the mean fitness of populations in which sex is rare or absent, provided that mutations affecting phenotypic traits have weak to moderate fitness effects. We then introduce genetic variation for the rate of sex, and show that the equilibrium rate of sex is an increasing function of the degree of mutational bias. Finally, we use the methods of [START_REF] Barton | A general model for the evolution of recombination[END_REF] and Charlesworth and Barton (1996) to express different components of selection for sex in terms of the effect of sex on the mean fitness and additive variance in fitness among offspring, and show that these expressions provide correct predictions when selection is sufficiently weak.

METHODS

Life cycle. The different parameters and variables of the model are summarized in Table 1. We consider a population of N haploid organisms with discrete generations.

Each individual may generate a fraction of its offspring asexually (by mitosis), the remaining fraction being produced sexually. In the last case, gametes are produced by mitosis and fuse at random in the population to form zygotes, which immediately undergo meiosis to produce haploid juveniles. We will first consider that all individuals invest equally into sexual reproduction, the parameter σ representing the rate of sex in the population (proportion of sexually produced offspring): σ = 0 corresponds to obligate asexual reproduction, and σ = 1 to obligate sex. In a second step (described below), we will introduce genetic variation for the rate of sex. We assume that individuals are hermaphroditic (generating both male and female gametes) and produce very large (effectively infinite) numbers of juveniles, among which N are sampled randomly to form the next adult generation (note that hermaphroditic haploid individuals occur in some species of mosses, ferns and algae).

Throughout the paper, fitness W denotes the overall fecundity of an individual and depends on the values of n quantitative phenotypic traits under stabilizing selection, represented by the vector z = (z 1 , z 2 , . . . z n ). In the following, we use greek letters α, β, γ... to denote phenotypic traits, while latin letters i, j, k... will denote loci. We assume that each phenotypic trait can be decomposed into a genetic and an environmental component:

z α = g α + e α ( 1
)
where g α is the individual's genetic contribution to trait α ("breeding value"), and where the environmental effect e α is independent of the genotype of the individual and is sampled from a Gaussian distribution with mean 0 and variance V e (the same for all traits). Average phenotypes and breeding values in the population are denoted z α and g α (with z α ≈ g α when the population is sufficiently large). As we assume no genotype × environment interaction, the variance of trait α is given by:

V α = V g, α + V e (2)
where V g, α is the genetic variance for trait α (variance of g α ). The genetic covariance between traits α and β (covariance between g α and g β ) will be denoted C g, αβ . Finally, X will denote the expected value of the quantity X at mutation-selection-drift equilibrium: for example, V g, α is the average genetic variance for trait α at equilibrium.

As we will see, some of our analytical results on the selective forces acting on the rate of sex do not depend on the specific shape of the fitness function. However, our simulation programs and some of our approximations assume an isotropic, Gaussianshaped fitness function around the phenotypic optimum, located at z = (0, 0, . . . 0):

W = exp - n α=1 z 2 α 2ω 2 , ( 3
)
where ω 2 represents the strength of selection. The mean fitness associated with a given genotype (obtained by averaging over the distribution of environmental effects e α ) is

given by:

W g = W g,max exp - n α=1 g 2 α 2V s (4)
where V s = ω 2 + V e , and where W g,max = (ω 2 /V s ) n/2 is the mean fitness of an optimal genotype (e.g., Lande, 1976a).

Genetic architecture of traits and mutational bias. We assume that selected traits are coded by loci with additive effects, so that

g α = j=1 g αj (5)
where g αj is the contribution of the allele at locus j on trait α. Loci are assumed biallelic (although some of our results on the selective forces acting on sex are valid under more general architectures), the alleles at each locus being denoted 0 and 1.

Assuming that an individual carrying allele 0 at all loci is at the phenotypic optimum, the contribution of locus j on trait α can be written as:

g αj = r αj X j , (6)
where X j is an indicator variable equal to 1 if the individual carries allele 1 at locus j (while X j = 0 otherwise), and r αj is the effect of allele 1 at locus j on trait α (note that r αj may be negative). The frequency of allele 1 at locus j is denoted p j , while q j = 1 -p j . At each locus, we assume that mutation occurs at the same rate u in both directions (from 0 to 1 and from 1 to 0), while U = u denotes the mutation rate on the whole set of loci affecting selected traits. As in previous works [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]Roze and Blanckaert, 2014), we introduce a parameter m measuring the degree of pleiotropy of mutations: each locus only affects a subset m (sampled randomly and independently for each locus) of the n traits under selection.

We assume that the distribution of r αj over all loci affecting trait α has average b and standard deviation a -the same for all traits -without any covariance between mutational effects on the different traits. From equation 4, the average deleterious effect of mutations on log W g (in an optimal genotype) is given by:

s d = 1 j=1 n α=1 r 2 αj 2V s = m (a 2 + b 2 ) 2V s . (7)
The parameter b represents the degree of mutational bias, since mutation tends to displace mean phenotypes away from the optimum when b = 0. In the following, mutational bias will be measured using a scaled parameter θ, defined as θ = b 2 / (a 2 + b 2 ) and varying between 0 and 1. For a given value of s d (mean fitness effect of mutations), the degree of mutational bias θ (between 0 and 1) while keeping the average fitness effect of mutations s d constant. This is equivalent to the approach used by [START_REF] Zhang | The anomalous effects of biased mutation revisited: mean-optimum deviation and apparent directional selection under stabilizing selection[END_REF], in which the variance of mutational effects decreases as the degree of mutational bias increases in order to maintain a constant mutational variance V M , defined as the per generation increase in phenotypic variance due to mutation (in our

model, V M = m n U (a 2 + b 2 ) = 2V s s d U/n).
Finally, we can note that while the average coefficient of epistasis (on fitness) between mutations is zero in the absence of bias (e.g., [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF], it becomes negative when θ > 0. Indeed, defining epistasis e as a deviation from additivity of mutational effects on log W g , we have (assuming that the number of loci is large):

e = - 2 ( -1) j =k n α=1 r αj r αk 2V s = -2ρ θ s d (8)
with ρ = m/n. In the extreme case when θ = 1 and ρ = 1 (all mutations have exactly the same phenotypic effect), epistasis becomes constant for all pairs of mutations and equals -2s d .

Change in phenotypic basis. Due to the symmetry of our model, average trait values and genetic variances at equilibrium should be the same for all traits, while mutational bias (θ > 0) will tend to displace the mean phenotype of the population in the direction of the (1, 1, . . . 1) vector. Although the effects of mutations on the different traits are not correlated, mutation generates genetic covariances C g, αβ between traits in the population (of the same magnitude for all pairs of traits), since individuals carrying more 1 alleles in their genomes tend to lie further in the direction of the (1, 1, . . . 1) vector. For analytical derivations, it is useful to define a new phenotypic basis in which the average mutational bias lies along the first axis and in which the genetic variance-covariance matrix is diagonal, thus eliminating covariances between traits (see Figure 1). This can be done by defining new breeding values g α as: .g., p. 380 in Anton, 2005). The fitnesses of genotypes in the new basis are still given by equation 4, replacing g α by g α . The average effect of mutations on g1 = g 1 / √ 2V s is given by: b1 = 1

g 1 = 1 √ n n β=1 g β (9) g α = 1 (α -1) α (α -1) g α - α-1 β=1 g β , α > 1 (10) (e
j=1 r1j = 1 j=1 1 √ n n β=1 r βj √ 2V s , (11) yielding
: b1 = ρ θ s d ( 12
)
where again ρ = m/n (note that equation 8 may thus be written as e = -2 b1 2 ). Due to the mutational bias, z 1 will tend to be positive, while the genetic variance along the first axis (V g,1 ) will be larger than along the other axes (see Figure 1).

Genetic control of the rate of sex. In order to explore the selective forces acting on reproductive mode, we will assume that a given individual may invest proportions s and 1 -s of its resources in sexual and asexual reproduction (respectively), and that genetic variation for s exists in the population. As in previous papers [START_REF] Michod | Deleterious mutations and selection for sex in finite, diploid populations[END_REF][START_REF] Otto | Differential selection between the sexes and selection for sex[END_REF]Roze, 2014), we introduce a direct cost of sex c by assuming that the probabilities that an individual is the maternal parent of a juvenile through asexual and sexual reproduction are proportional to 1 -s and s/c, respectively (c = 1 in the absence of cost, while c = 2 corresponds to a twofold cost of sex). This cost may be caused by anisogamy (cost of males): for example c = 2 when half of the resources invested in sex are used to produce male gametes, assuming that the same amount of resources is needed to produce a female gamete and an asexual spore. Alternatively, the cost may result from the failure of gametes to find a partner (assuming that a proportion 1 -1/c of gametes are lost), or to extra energetic costs associated with gamete production compared with asexual spore production. The rate of sex σ of an individual is defined as the proportion of sexually produced individuals among its maternally produced offspring, given by:

σ = s c (1 -s) + s (13)
(σ = s in the absence of cost). We assume that, like the other traits, investment in sex can be decomposed into an additive genetic and an environmental component:

s = s + g s + e s ( 14
)
where s is the average investment in sex in the population, g s = i g si (g si being the effect of the allele present at locus i on investment in sex) while e s is sampled from a centered Gaussian distribution with variance V e, s . These equations assume that the distribution of values of s in the population is not too close to 0 or 1 (otherwise the assumption of additivity may not hold, as s cannot be lower than 0 or higher than 1).

As above, the variance in s in the population is given by V g, s + V e, s (where V g, s is the variance in g s ). Throughout the paper, we will assume that loci affecting investment in sex do not affect the traits under stabilizing selection.

Assuming that the variance in s in the population is sufficiently small, the rate of sex σ may also be decomposed into an additive genetic and an environmental component:

σ = σ + g σ + e σ ( 15
)
where σ is the mean rate of sex. From equation 13 and 14 (and assuming that g s and e s are small, of order ), we have:

σ ≈ s c (1 -s) + s , V g, σ ≈ c 2 [c (1 -s) + s] 4 V g, s (16)
(to leading order in ) where V g, σ is the genetic variance for the rate of sex (variance of g σ ).

Simulation programs. Our individual-based simulation programs (written in C++) are available from Dryad, and described in Supplementary File S1. The genome of each individual consists in a single linear chromosome with map length R (average number of cross-overs at meiosis). The biallelic loci affecting the n traits under stabilizing selection are equally spaced along the chromosome, each of these loci affecting a subset of m randomly sampled traits as described above. Investment in sexual re-production s is coded by s multiallelic loci (with an infinite number of possible alleles per locus), which are also equally spaced along the chromosome (see Figure 2); assuming multiallelic loci ensures that all possible rates of sex between 0 and 1 may be achieved even when the number of loci affecting investment in sex is low. Mutational effects at these loci are sampled from a centered Gaussian distribution with variance a 2 s (the mutational effect being added to the value coded by the allele before mutation). Investment in sex s is obtained by summing allelic effects at all these s loci, and adding an environmental component drawn from a centered Gaussian distribution with variance V e, s (if the value obtained is lower than 0 or higher than 1, it is then set to 0 or 1). In a different version of the program the s multiallelic loci do not affect investment in sex (which is fixed), but correspond to neutral loci which are used to estimate the effective population size N e . For this, diversity at each of these neutral loci is computed as D = 1i p 2 i (where p i is the frequency of allele i), and the effective population size is estimated by N e ≈ D/ 2µ 1 -D , where D is the average diversity over neutral loci and generations, and µ the mutation rate at each neutral locus (generally fixed to 10 -3 ). Simulations with a fixed rate of sex generally lasted 10 5 generations, while simulations in which investment in sex was free to evolve lasted 2 × 10 6 generations (however the rate of sex generally reached an equilibrium within the first 5 × 10 5 generations).

EFFECT OF MUTATIONAL BIAS ON MEAN FITNESS

In this section we assume that the rate of sex σ is fixed, and explore the effect of σ and of mutational bias on mean fitness. The mutation load L measures the reduction in mean fitness of the population due to the presence of deleterious alleles, and is defined as:

L = 1 - W W g,max , ( 17
)
where W is mean fitness and W g,max the fitness of an optimal genotype. Throughout this section, we assume an isotropic, Gaussian-shaped fitness function (equation 3).

Assuming that the variance in log-fitness in the population is small and that population size is large, we have (see Supplementary File S2):

L ≈ 1 -exp - 1 2V s n α=1 V g, α + g α 2 . ( 18
)
In the absence of mutational bias, the effect of deviations of mean phenotypes from their optimal values (the term in g α 2 in equation 18) is proportional to 1/N e , and should thus remain small when N e is sufficiently large (Lande, 1976b;Charlesworth, 2013b). However, in the presence of mutational bias, drift may cause substantial deviations of 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]. Simple approximations for the load can be obtained in the regime where selection is negligible relative to drift at all loci. Assuming that the variance of g α due to drift is small, we have g α 2 ≈ g α 2 , while g α = j=i r αj p j in our biallelic model. Using equation 12, and the fact that p j = 1/2 under symmetric mutation when the effect of selection at locus j is neglected, one obtains:

1 2V s n α=1 g α 2 ≈ 1 4 b1 2 (19)
where b1 = √ ρ θ s d is the (scaled) magnitude of mutational bias (along the z 1 axis).

Furthermore, linkage disequilibria between loci should be close to zero on average when selection is sufficiently weak, in which case the genetic variance for trait α is given by:

V g, α ≈ j=1 r 2 αj p j q j (20)
(e.g., Lynch and Walsh, 1998). Given that p j q j ≈ N u/ (1 + 4N u) at mutation-drift balance, one obtains from equation 20:

1 2V s n α=1 V g, α ≈ s d N U 1 + 4N u , (21)
finally giving:

L ≈ 1 -exp -s d N U 1 + 4N u - 1 4 b1 2 (22)
Equation 22 is equivalent to equation 8 in Roze and Blanckaert (2014) that selection may have significant effects on mean trait values even when p i q i at each locus is mainly controlled by mutation and drift [START_REF] Robertson | A theory of limits in artificial selection[END_REF]Campbell, 1984;Barton, 1989;Charlesworth, 2013a). Based on this, it is possible to derive a better approximation for low s d by taking the effect of selection on g α into account, while still neglecting the effect of selection on genetic variance (and neglecting linkage disequilibria). This yields (see Supplementary File S2 for derivation):

L ≈ 1 -exp    -s d N U 1 + 4N u - b1 2 4 1 + s d n [1 + θ (m -1)] N 1+4N u 2    . ( 23
)
Approximations for the regime where genetic variances are significantly affected by mutation, selection and drift are more difficult to obtain. Under very strong selection against mutant alleles (so that the contribution to future generations of individuals deviating from the optimum can be neglected), the mutation load becomes (for both sexual and asexual populations):

L ≈ 1 -e -U (24)
(e.g., [START_REF] Kimura | The mutational load with epistatic gene interactions in fitness[END_REF]. Under sexual reproduction, equation 24 also holds under weaker selection in the absence of mutational bias (θ = 0), as long as drift and linkage disequilibria may be neglected (e.g., Bürger, 1998, Supplementary File S2).

In the case of an asexual population, an expression for the load at mutation-selection balance (still in the absence of mutational bias, and neglecting drift) can be obtained assuming a Gaussian distribution of trait values in the population: (Lande, 1980a;Roze and Blanckaert, 2014). Generalizing these expressions to introduce mutational bias is not straightforward in the context of our biallelic model, as the degree of mutational bias changes depending on the position of mean phenotypes;

L ≈ 1 -exp - n 2 U s d (25)
however, previous studies have shown that the effect of mutational bias is generally small in regimes where drift is negligible [START_REF] Waxman | The anomalous effects of biased mutation[END_REF][START_REF] Zhang | The anomalous effects of biased mutation revisited: mean-optimum deviation and apparent directional selection under stabilizing selection[END_REF]. In Supplementary File S2, we show that a deterministic approximation for the load in a sexual population under the maximum level of bias (θ = 1) is given by:

L ≈ 1 -exp - 4ρU -s d + s d (8ρU + s d ) 8ρ (26)
(see Supplementary File S2 for the same expression in terms of b1 , s d and U ).

Figure 3 shows the equilibrium mutation load as a function of s d , for different rates of sex and levels of mutational bias. In the absence of mutational bias (θ = 0), the load is generally higher in sexual (σ = 1) than in asexual (σ = 0) populations, due to the fact that recombination breaks combinations of alleles with compensatory effects (recombination load). This pattern reverses for high values of s d (s d = 0.1 in Figure 3), as the frequency of deleterious alleles is increased by Hill-Robertson effects in asexual populations. While the effect of mutational bias (with θ = 0.1) on the mean fitness of sexual populations stays modest, it greatly increases the load of asexual populations for small values of s d (between 10 -5 and 10 -3 ) -see Supplementary

Figure S1 for results under stronger bias (θ = 0.5). Supplementary Figure S2 shows that this increase in L is caused by deviations of mean phenotypes from the optimum, due to the combined effects of mutational bias and drift. Indeed, Figure 3 shows that the effective population size of asexual populations (estimated from the average diversity at neutral loci, see Methods) is greatly reduced by background selection effects.

As shown by Figure 4, equation 22 correctly predicts the increase in load caused by mutational bias at very low values of s d , but rapidly overestimates L as s d increases, as it neglects the effects of selection (see Supplementary Figure S3 for the relative effects of genetic variance and of deviations of mean phenotypes from the optimum).

In the case of sexual populations, equation 23 provides better predictions (dotted curves in Figure 4) but still fails when s d is not very small, as it neglects the effect of selection on genetic variances. In agreement with previous results [START_REF] Waxman | The anomalous effects of biased mutation[END_REF][START_REF] Zhang | The anomalous effects of biased mutation revisited: mean-optimum deviation and apparent directional selection under stabilizing selection[END_REF], we find that in sexual populations, the effect of mutational bias stays rather small in the deterministic regime (N s d 1). Very strong levels of bias (θ = 0.5, 1) decrease the load in this regime, this effect being correctly predicted by our deterministic approximation for θ = 1 (equation 26): this is due to the fact that mutational bias generates negative epistasis (on average) between deleterious alleles (equation 8), reducing the mutation load of sexual 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]. Figure 5 shows that the effect of mutational bias increases as the number of loci increases (allowing stronger deviations from the fitness optimum) and as population size decreases (see Supplementary Figures S4 -S6 for results under stronger bias and for the relative effects of genetic variance and of deviations of mean phenotypes from the optimum on the load). The effects of the degree of pleiotropy of mutations m and of the total number of selected traits n are shown on Figure 6. The mutation load increases with the degree of pleiotropy (Figure 6, top panels): indeed, the magnitude of mutational bias b1 increases with ρ θ (with ρ = m/n, equation 12). Increasing m while keeping ρ θ constant has only little effect on the load (Figure 6, middle panels), confirming that m mostly affects the load through its effect on b1 . Finally, Figure 6 shows that increasing n while keeping m/n (and thus b1 ) constant has little effect on the load in sexual populations, while it increases the load of asexual populations due to stronger deviations of mean phenotypes from the optimum (see Supplementary Figures S7). Indeed, increasing the dimensionality n of the fitness landscape reduces the chances that a deleterious allele can be compensated by mutations at other loci, and thus enhances the effect of mutational bias in asexuals.

Overall, these results show that the combined action of mutational bias and genetic drift may greatly reduce the mean fitness of asexual populations when the average fitness effect of mutations is small to moderate, this increase in load being maximized for intermediate strengths of selection against deleterious alleles s d , higher values of pleiotropy m/n, number of selected traits n and number of loci , and for lower values of population size N . In the next section, we will see how this translates into selection on modifier genes affecting the rate of sex of individuals.

EVOLUTION OF SEX

Analytical approximations. Expressions for the effect of selection on the rate of sex are derived in Supplementary File S3, assuming weak selection, a Gaussian distribution of traits affecting fecundity (the z α 's) and a low variance for the rate of sex in the population. Under these assumptions, the change in the mean rate of sex over one generation (∆σ) can be decomposed into two terms, representing the effect of the cost of sex (direct selection), and indirect selection caused by the effect of sex on genetic associations between loci affecting fecundity:

∆σ = ∆ cost g σ + ∆ ind g σ . ( 27
)
As shown in Supplementary File S3, ∆ cost g σ ≈ β cost V g, σ , where

β cost ≈ - c -1 1 + (c -1) σ (28)
represents the direct selection gradient (selecting against sex when c > 1). Indirect selection in turn decomposes into two terms, sometimes called the "short-term" and "long-term" effect of breaking genetic associations (e.g., Agrawal, 2006):

∆ ind g σ = ∆ short g σ + ∆ long g σ . (29)
The short-term effect is due to the fact that, in the presence of epistatic interactions, breaking genetic associations between loci affects the mean fitness of offspring. Under our isotropic fitness function (equation 3), and assuming that phenotypes are measured in a basis that eliminates covariances between traits (the basis defined by equations 9 and 10), ∆ short g σ is given by:

∆ short g σ ≈ n α=1 ∂ ln W ∂V g, α M g,σαα (30)
where M g,σαα is the third moment E (g σ -g σ ) (g α -g α ) 2 (where E stands for the average over all individuals). A more general expression for arbitrary fitness function is given in the Appendix (see Supplementary File S3 for derivation). Under stabilizing selection, ∂ ln W /∂V g, α is negative (mean fitness decreases as the genetic variance for selected traits increases). Furthermore, selection tends to generate associations (linkage disequilibria) between alleles at different loci with compensatory effects on selected traits, thereby reducing V g, α . By breaking these associations, sex increases the genetic variance among offspring: therefore, the genetic variance tends to be higher among individuals that engage more in sex (i.e., with higher values of g σ ) than among individuals that engage less in sex, translating into a positive value of M g,σαα . The term representing the short-term effect (∆ short g σ ) is thus negative, corresponding to the short-term cost of breaking genetic associations that have been generated by selection -one can show that this term is equivalent to the term in δV g in Charlesworth's (1993) recombination modifier model, see also Appendix 2 of [START_REF] Barton | A general model for the evolution of recombination[END_REF].

The long-term effect stems from the fact that increasing the genetic variance among offspring allows a better response to directional selection, and can be written as:

∆ long g σ ≈ n α=1 ∂ ln W ∂z α C g,σα (31)
where

C g,σα = E [(g σ -g σ ) (g α -g α )]
is the genetic covariance between the rate of sex (1993). Charlesworth (1993) and [START_REF] Barton | A general model for the evolution of recombination[END_REF] showed how the short-term and longterm effect can be expressed in terms of mean trait values and genetic variances for selected traits in a recombination modifier model, neglecting the effects of genetic drift on genetic associations and using a quasi linkage equilibrium (QLE) approximation.

An equivalent derivation for the case of the present model is given in Supplementary File S3, the main results being summarized in the Appendix. For this, we assume an infinite population size, large number of loci affecting fecundity, weak selection and low variance for the rate of sex in the population; we also assume that the rate of sex is not too low (for the QLE approximation to hold). One obtains that ∆

short g σ ≈ β short V g, σ , ∆ long g σ ≈ β long V g, σ ,
where the short and long-term selection gradients are given by:

β short ≈ - 1 2V s 2 r h,1 σ 2 n α=1 V g, α 2 - z 1 2 V g, 1 2 V s , (32)
β long ≈ 1 r h,2 σ - 1 r h,1 1 σ 2 z 1 2 V g, 1 2 V s 3 . ( 33
)
Equations 32 and 33 assume that traits are measured in the phenotypic basis given by equations 9 and 10, so that only the first phenotypic trait (with average z 1 and genetic variance V g, 1 2) is under directional selection. The terms r h,1 and r h,2 that appear in the denominators of β short and β long correspond to harmonic mean recombination rates among loci. Defining r ijk as the probability that at least one recombination event occurs at meiosis between a locus i affecting investment in sex and loci j and k affecting selected traits, r h,1 is the harmonic average of r ijk over all possible triplets of loci i, j

and k, while r h,2 is the harmonic average of r ij r ijk , where r ij is the recombination rate between loci i and j. The maximum possible values of r h,1 and r h,2 (obtained for the case of freely recombining loci) are thus 3/4 and 3/8, respectively.

Equations 32 and 33 indicate that both the short-term and long-term selection gradients increase as the mean rate of sex in the population σ decreases, β long increasing more rapidly (due to the term in 1/σ 3 ). However, both expressions diverge as σ tends to zero, due to the QLE approximation. Equation 33 also shows that the long-term effect vanishes in the absence of mutational bias (z 1 = 0). The genetic architecture of investment in sex affects β long and β long through r h,1 and r h,2 . Provided that the number of loci affecting fecundity is large and that their distribution over the genome is relatively uniform, the harmonic averages of r ijk and r ij r ijk over all j and k should be similar for all loci i affecting investment in sex, and the indirect selection gradient should thus be little affected by the number of loci coding for the rate of sex.

As we have seen in the previous section, it is difficult to obtain general analytical expressions for mean trait values (z 1 ) and genetic variances (V g, α ) at mutationselection-drift equilibrium under mutational bias, for arbitrary values of s d and σ, and we were thus not able to express the mean rate of sex in the population at equilibrium in terms of the different parameters of the model. One can note, however, that the approximations above for the short and long term selection gradients can be expressed in terms of the effect of sex on the average and variance in fitness among offspring, that could (at least in principle) be measured from an experimental population [START_REF] Barton | A general model for the evolution of recombination[END_REF]Charlesworth and Barton, 1996). Indeed, denoting W sex and W asex the mean fitness of sexually and asexually produced offspring (respectively), and Var A, sex (ln W ),

Var A, asex (ln W ) the additive variance in log fitness among sexually and asexually produced offspring, we have (see Supplementary File S3):

β short ≈ ∆ 1 r h,1 σ , β long ≈ 1 r h,2 σ - 1 r h,1 ∆ 2 σ ( 34
)
with:

∆ 1 = ln W sex -ln W asex , (35)
∆ 2 = Var A, sex (ln W ) -Var A, asex (ln W ) . ( 36
)
Equations 34 -36 are valid in principle for any shape of the fitness function (not necessarily Gaussian), as long as selection is sufficiently weak and the number of selected loci is sufficiently large. However, as the previous results, they assume that genetic associations remain small (QLE approximation), causing them to diverge as the mean rate of sex in the population tends to zero.

Simulation results. Figure 7 shows that, in agreement with the discussion above, the number of loci affecting investment in sex has very little effect on the mean rate of sex in the population (σ) at equilibrium (the numbers 9, 99 and 999 were chosen so that the number of loci affecting fecundity between two loci affecting sex is 1000, 100 and 10, respectively -see Figure 2). As shown by Figures 789, the population evolves towards asexuality in the absence of mutational bias (θ = 0), while increasing the magnitude of mutational bias b1 (by increasing either θ or m/n, see equation 12)

increases the equilibrium rate of sex. Higher rates of sex evolve under higher values of the mutation rate U , larger numbers of selected loci and lower values of population size N , due to stronger effects of mutational bias (Figures 89). Similarly, increasing the dimensionality of the fitness landscape n while keeping m/n constant (so that b1 stays constant) enhances the effect of mutational bias in asexuals (Figures 6, S7), favoring higher rates of sex (Figure 8). The mean fitness effect of deleterious alleles has a non-monotonic effect on selection for sex, the equilibrium rate of sex being maximized for intermediate values of s d (Figure 9). The genome map length R also has a non-monotonic effect on the equilibrium rate of sex (Figure 9): up to a certain point, increasing linkage favors sex since the long-term benefit of sex increases faster than the short-term cost as linkage becomes tighter (as can be seen from equations 32 and 33, and the fact that r h,2 decreases faster than r h,1 as recombination rates decrease).

However, indirect selection vanishes when R tends to zero (since sex becomes genetically equivalent to asexual reproduction), in which case the rate of sex evolves towards zero when sex is costly -Figure 9 shows that low rates of sex may be maintained in the population, probably due to hitchhiking effects between loci affecting investment into sex and loci affecting selected traits. Finally, Figure 8 shows that higher rates of sex are maintained in the absence of a direct cost of sex (c = 1), although the rate of sex still evolves towards zero when mutational bias is absent (θ = 0).

Our simulation program was modified in order to test the validity of the QLE approximations shown above (equations 32 -34) for different values of σ. In this modified version, we introduce genetic variation for investment in sex but constrain σ to stay in a given range by sampling the value of alleles at loci affecting sex after mutation from a uniform distribution with variance a 2 s , without adding the value of the allele before mutation. The short and long-term selection gradients were estimated from equations 30 and 31 (divided by V g, σ ), using equations A2 and A3 and measuring the moments z α , V g, α , M g,σαα and C g,σα for all traits α. For this, the value of g σ was estimated for each individual from the average rate of sex σ of 100 clonally produced offspring (all with different environmental components of investment in sex e s ), given by equation 13. The terms ∆ 1 and ∆ 2 of equation 34 were also measured every 100 generations by producing a pool of offspring by sexual reproduction and another pool by asexual reproduction, and measuring the mean fitness and additive variance in log fitness within each pool of offspring. The additive variance in log fitness was estimated from the covariance in log fitness Cov(ln W ) between sexually (or asexually) produced offspring and their own sexually produced offspring, using Var A (ln W ) = 4Cov(ln W ) -Var(ln W ) (Lynch and Walsh, 1998, Supplementary File S3). Figure 10 shows that the QLE approximation provides correct predictions of the indirect selection gradients when selection is sufficiently weak (s d = 10 -4 , for the parameter values used in Figure 10) and for intermediate rates of sex (while the QLE expressions diverge as σ approaches zero). Discrepancies appear for s d = 10 -3 , however, and become more important for s d = 10 -2 . These discrepancies are probably due to a breakdown of the different assumptions used to derive equations 32 -36 (e.g., weak genetic associations, negligible effect of associations involving more than 2 or 3 loci, distribution of breeding values close to a Gaussian distribution), and possibly also to the effect of drift on genetic associations (through the Hill-Robertson effect), which is not taken into account in our analysis.

DISCUSSION

Epistasis and drift are the two major sources of genetic associations that have been considered in theoretical studies on the benefits of sex and recombination. Epistasis may favor recombination when it is negative on average, that is, when the fitness effect of a deleterious allele is increased by the presence of other deleterious alleles at other loci, or conversely when the fitness effect of a beneficial allele is decreased by the presence of other beneficial alleles in the genome. However, epistatic interactions also generate a short-term cost for recombination (since recombinant offspring tend to have a lower mean fitness than their parents in a constant environment), so that high rates of recombination can only be favored when epistasis is weak relative to the strength of selection, and not too variable across loci [START_REF] Barton | A general model for the evolution of recombination[END_REF][START_REF] Otto | Deleterious mutations, variable epistatic interactions, and the evolution of recombination[END_REF].

Epistatic interactions (on fitness) arise naturally in models of selection acting on quantitative phenotypic traits. In agreement with the results mentioned above, Gaussian (or quadratic) stabilizing selection around a fixed optimum in an infinite population is expected to disfavor recombination in the absence of mutational bias (Charlesworth, 1993). Indeed, at equilibrium the mean phenotype of the population is centered on the optimum, in which case epistasis between deleterious alleles is zero on average, with a given variance [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF] -epistasis between two alleles displacing the phenotype in the same direction is negative (due to the negative curvature of the fitness function), while epistasis between alleles having opposite (compensatory) effects on the phenotype is positive. Away from the optimum, epistasis between deleterious alleles is negative on average (while epistasis between beneficial alleles is also negative, e.g., [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF], generating a deterministic advantage for recombination (Charlesworth, 1993).

Our simulation results confirm that, in the absence of mutational bias on phenotypic traits, populations evolve towards obligate asexuality when the phenotypic optimum remains constant over time, even when population size is finite. This stands in contrast with previous simulation results assuming fixed epistasis across loci (always negative or always positive), that found only minor effects of epistasis compared with the stochastic (Hill-Robertson) effects that favor recombination in initially asexual (or non recombining) populations [START_REF] Keightley | Interference among deleterious mutations favours sex and recombination in finite populations[END_REF]. When mutational bias is included in the model, however, positive rates of sex are maintained in the population at equilibrium. Indeed, mutational bias tends to displace mean phenotypes away from the optimum (thereby increasing the mutation load), this effect being stronger in asexual populations in which the variance in fitness may be greatly lowered by negative associations between loci, reducing their ability to respond to directional selection. Extending [START_REF] Barton | A general model for the evolution of recombination[END_REF] QLE analysis to our model, we obtained deterministic approximations for the short and long-term indirect selection gradients acting on sex in terms of mean trait values and genetic variances, and showed that these approximations provide reasonable predictions when selection acting at the different loci is sufficiently weak and when the rate of sex is not too low (Figure 10). This implies that, in this parameter range, selection for sex is mainly driven by negative linkage disequilibria caused by epistasis (although drift may play a significant role by increasing the distance between the mean phenotype and the optimum, and therefore the magnitude of directional selection). The Hill-Robertson effect may become more important in parameter ranges where the QLE approximation fails (strong selection and/or low rate of sex); however, the lack of suitable analytical method to cover such regimes makes it difficult to assess its relative effect.

Recent experimental evolution studies showed that higher rates of sex or outcrossing may evolve in populations adapting to a new environment or coevolving with a pathogen, possibly through the generation of advantageous genotypes by recombination and segregation [START_REF] Becks | Higher rates of sex evolve in spatially heterogeneous environments[END_REF]Agrawal, 2010, 2012;[START_REF] Morran | Running with the Red Queen: host-parasite coevolution selects for biparental sex[END_REF][START_REF] Luijckx | Higher rates of sex evolve during adaptation to more complex environments[END_REF]. In adapting populations of monogonont rotifers, Becks and Agrawal (2012) showed that sexually produced offspring tend to have a lower mean fitness and a higher variance in fitness than asexually produced offspring, in agreement with predictions from models with concave fitness functions such as the one used in this paper. However, how to relate the effect of sex on the mean and variance in fitness of offspring with the strength of indirect selection for sex is not immediately obvious. Transposing [START_REF] Barton | A general model for the evolution of recombination[END_REF] and Charlesworth and Barton's (1996) analysis of recombination modifier models to our sex modifier model, we showed that simple relations exist between the short and long-term selection gradients for sex and the effect of sex on the fitness of offspring (equations 34 -36). However, several important caveats must be noted:

(i) these relations only hold in the QLE regime, and thus break down when the rate of sex in the population is low; (ii) they depend on average recombination rates between loci affecting fitness and loci affecting the rate of sex (through r h,1 and r h,2 in equation 34), which are generally unknown (although lower bounds for selection gradients can be obtained by replacing these terms by their values under free recombination, i.e. r h,1 = 3/4 and r h,2 = 3/8); (iii) the long-term selection gradient is expressed in terms of the effect of sex on the additive variance in fitness among offspring, which will generally be more difficult to measure than the variance in fitness. Nevertheless, estimations of the effect of sex on the mean and variance in fitness among offspring still convey important information on the existence and sign of short and long-term selection gradients on sex (e.g., [START_REF] Peters | Liberating genetic variance through sex[END_REF][START_REF] Sharp | Evolution of sex: using experimental genomics to select among competing theories[END_REF].

For a given genomic mutation rate U , our model predicts that increasing the dimensionality of the fitness landscape n increases selection for sex (Figure 8). Indeed, the variance of epistasis between mutations decreases as n increases [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF], epistasis vanishing as n tends to infinity, since mutations become orthogonal in this limit (without any possible compensatory effect). In other words, strong epistatic interactions (in particular, compensatory effects between deleterious alleles) are more likely to occur when the dimensionality of the fitness landscape is low, and these strong interactions tend to favor asexual reproduction (that can maintain coadapted multilocus genotypes). However, we can note that our model assumes that all loci have the same probability of affecting any trait: under a more modular genetic architecture where different sets of loci affect different sets of traits (modular pleiotropy, e.g., [START_REF] Welch | Modularity and the cost of complexity[END_REF][START_REF] Chevin | Fisher's model and the genomics of adaptation: restricted pleiotropy, heterogenous mutation, and parallel evolution[END_REF][START_REF] Chebib | What affects the predictability of evolutionary constraints using a G-matrix? The relative effects of modular pleiotropy and mutational correlation[END_REF], the magnitude of epistatic interactions may be more dependent on the average number of traits coded by a given module than on the total number of selected traits, which may lead to different results. In general, the range of realistic values for the dimensionality of fitness landscapes remains difficult to assess: while a large number of traits in an organism may be under selection, many of those traits are probably correlated, reducing the effective dimensionality of the landscape (Martin and Lenormand, 2006). In VSV and φX174 viruses, the effective number of selected traits was estimated to be around 10 and 45 (respectively) based on predictions from Fisher's geometric model on the relation between N e and population mean fitness [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], but this number may be much higher in multicellular eukaryotes.

As we have seen, mutational bias is required for sex to be favored in a constant environment. Some evidence for mutational bias on quantitative traits has been obtained from Drosophila and Caenorhabditis elegans (e.g., [START_REF] Santiago | The distribution of spontaneous mutations on quantitative traits and fitness in Drosophila melanogaster[END_REF][START_REF] Lyman | Effects of single P -element insertions on bristle number and viability in Drosophila melangaster[END_REF][START_REF] Keightley | EMS-induced polygenic mutation rates for nine quantitative characters in Drosophila melanogaster[END_REF][START_REF] Ostrow | Mutational bias for body size in rhabditid nematodes[END_REF][START_REF] García-Dorado | Properties of spontaneous mutations affecting quantitative traits[END_REF]; however, how to relate these data with the parameter θ measuring bias in our model is not immediately obvious. In particular, a downward mutational bias is often observed on traits that may be seen as fitness components, but such a bias is expected in our model at the optimum even when θ = 0 (since fitness can only decrease at an optimum). Traits that have a less direct relation with fitness sometimes show mutational bias (e.g., metabolite pool size, [START_REF] Davies | The mutational structure of metabolism in Caenorhabditis elegans[END_REF], sometimes not (e.g., mitotic spindle traits, [START_REF] Farhadifar | Mutation is a sufficient and robust predictor of genetic variation for mitotic spindle traits in Caenorhabditis elegans[END_REF] but it is again difficult to relate such measures to θ, since the relation between these traits and fitness is generally poorly known.

Information on θ may rather be obtained from the distribution of fitness effects of mutations. Indeed, bias causes mutation to push phenotypic traits in a given direction away from the optimum, so that the proportion of beneficial mutations should always stay below 0.5, even for small-effect mutations occurring in a non-optimal genotype.

By contrast, in the absence of bias the proportion of beneficial mutations tends to 0.5 as one moves away from the optimum, the convergence to 0.5 being faster for smallereffect mutations. As a consequence of this high rate of compensatory mutations, drift load generally stays mild in the absence of bias unless population size is very small (the load being roughly proportional to n/N , e.g., Lande, 1980b;[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], while it may reach much higher values when mutational bias is present, as shown in the present paper. Compensatory mutations has been best studied in model organisms such as bacteriophages, bacteria, nematodes and yeasts where they were shown to be common (e.g., [START_REF] Levin | Compensatory mutations, antibiotic resistance and the population genetics of adaptive evolution in bacteria[END_REF][START_REF] Poon | The rate of compensatory mutation in the DNA bacteriophage φX174[END_REF][START_REF] Estes | Fitness recovery and compensatory evolution in natural mutant lines of C. elegans[END_REF][START_REF] Szamecz | The genomic landscape of compensatory evolution[END_REF]. However, more work is needed to better understand how the rate of compensatory mutations changes with the degree of maladaptation of individuals, in order to gain more insights on realistic levels of mutational bias (as modeled here).

Finally, we can note that the equilibrium rate of sex in the population generally stays small when the cost of sex is moderate to strong (Figures 8,9), the highest rates of sex being always achieved under complete bias (θ = 1), that is, when compensatory mutations are not possible. Similarly, low levels of costly sex are also maintained in most cases in models on the evolution of sex due to deleterious mutations without epistasis [START_REF] Michod | Deleterious mutations and selection for sex in finite, diploid populations[END_REF][START_REF] Otto | Differential selection between the sexes and selection for sex[END_REF]Roze, 2014). Exploring to what extent higher levels of sex may be maintained in models including environmental change would thus be of interest, and will be the subject of future work.

APPENDIX: GENERAL QLE RESULTS

Assuming that the distribution of phenotypic traits affecting fecundity in the population is approximately Gaussian and that selection is weak, a general expression for indirect selection on the rate of sex is given by (see Supplementary File S3 for derivation):

∆ ind g σ ≈ n α=1 ∂ ln W ∂z α C g,σα + α≤β ∂ ln W ∂C g,αβ M g,σαβ (A1)
where the second sum is over all possible pairs of selected traits, including α = β. Equation A1 is equivalent to Charlesworth's (1993) decomposition of the selection gradient for a recombination modifier allele into two terms (equation A10 in Charlesworth, 1993, see also Appendix 2 of [START_REF] Barton | A general model for the evolution of recombination[END_REF]. The first term of equation A1 (equivalent to the term in δz in Charlesworth, 1993) represents indirect selection caused by the effect of sex on mean phenotypes. With our Gaussian, isotropic fitness function (equation 3),

we have:

∂ ln W ∂z α = - z α V g, α + V s , ( A2
)
which is approximately -z α /V s when selection is weak (V g, α V s ). In our model, directional selection occurs along the axis corresponding to the direction of the mutational bias, and therefore only the first term of the sum (for α = 1) will contribute when phenotypes are measured in the basis defined by equations 9 and 10. The second term of equation A1 (equivalent to the term in δV g in Charlesworth, 1993) represents indirect selection caused by the effect of sex on the genetic variance-covariance matrix:

∂ ln W /∂C g,αβ describes how mean fitness is affected by the genetic covariance between traits α and β, while the third moment

M g,σαβ = E [(g σ -g σ ) (g α -g α ) (g β -g β )]
(where E stands for the average over all individuals) describes to what extent the genetic covariance between traits α and β differs between subsets of the populations with different rates of sex. As shown in Supplementary File S3, under an isotropic, Gaussian fitness function and measuring phenotypes in the basis defined by equations 9 and 10, only the terms with α = β differ from zero, and the second term of equation A1 thus becomes n α=1 ∂ ln W /∂V g, α M g,σαα . The selection gradient ∂ ln W /∂V g, α measures the strength of stabilizing selection on trait α, and is given by:

∂ ln W /∂V g, α = - 1 2 (V g, α + V s ) + 1 2 ∂ ln W ∂z α 2 ≈ - 1 2V s 1 - z α 2 V s (A3)
where again the term z α 2 will differ from zero only for the first phenotypic trait in the basis defined by equations 9 and 10. Note that the second term of equation A1 does not appear in classic expressions describing the effect of selection on correlated characters [START_REF] Lande | Quantitative genetic analysis of multivariate evolution, applied to brain:body size allometry[END_REF], as these assume a multivariate Gaussian distribution of phenotypic traits. Here we cannot assume that the joint distribution of the rate of sex σ and of the traits affecting fecundity is multivariate Gaussian: in particular, sex tends to increase V g, α by breaking negative genetic associations (linkage disequilibria between alleles with compensatory effects on trait α), generating a positive third moment M g,σαα .

Following Charlesworth (1993) and [START_REF] Barton | A general model for the evolution of recombination[END_REF], the moments C g,σα and M g,σαβ that appear in equation A1 may be expressed in terms of the genetic variance for the rate of sex V g, σ and genetic variances (and covariances) for selected traits using a QLE argument. The derivation (shown in Supplementary File S3) supposes that selection is weak relative to effective recombination rates between loci (and thus that the rate of sex is not too low), so that linkage disequilibria remain small. Furthermore, it neglects the effects of genetic associations involving more than three loci.

Under these assumptions, one obtains for M g,σαβ : that eliminates covariances between traits, we have (assuming that the number of loci affecting selected traits is large):

M g,σαβ ≈ - ∆ sel D αβ r h,1 σ 2 V g, σ (
∆ sel D αβ ≈ (1 + I αβ ) ∂ ln W ∂C g,αβ - ∂ ln W ∂z α ∂ ln W ∂z β V g, α V g, β (A5)
where I αβ equals 1 if α = β, and 0 otherwise. Under an isotropic, Gaussian fitness function, it is possible to show that the term between brackets in equation A5 equals 0 when α = β, while it is approximately -1/V s when α = β (Supplementary File S3, equation A3). In this case, equations A3 -A5 yield equation 32 in the main text.

The QLE expression for the genetic covariance C g,σα that appears in the first term of equation A1 writes (see Supplementary File S3 for derivation):

C g,σα ≈ - 1 r h,2 σ - 1 r h,1 1 σ 2 n β=1 ∂ ln W ∂z β (∆ sel D αβ ) V g, σ (A6)
where r h,2 is defined in the main text. Under an isotropic, Gaussian fitness function, equations A2 -A5 yield equation 33 in the main text. Strength of stabilizing selection on phenotypic traits optimal values by the same amount, and generates a positive covariance among traits.

V s = ω 2 + V e Strength of
Traits g 1 and g 2 are defined by rotating the phenotypic basis (equations 9 and 10, red axes) so that the covariance between g 1 and g 2 is zero, while mutational bias displaces phenotypes along the g 1 axis.

……….

map length R locus affec2ng phenotypic traits z α (biallelic) locus affec2ng investment in sex s (mul2allelic) Investment in sex is coded by s multiallelic loci, which are also regularly spaced along the chromosome. When s = 1, the locus affecting investment in sex is located at the mid-point of the chromosome. of mutations m and numbers of selected traits n. In the middle panels, m θ is kept for n = 10, 20 and 50, respectively.

• σ = 0 • σ = 0.1 • σ = 1 1 10 -2 10 -4 s d 0.2 0.4 0.6 0.8 1.0 L θ = 0 • σ = 0 • σ = 0.1 • σ = 1 1 10 -2 10 -4 s d 0.2 0.4 0.6 0.8 1.0 L θ = 0.1 • σ = 0 • σ = 0.1 • σ = 1 1 10 -2 10 -4
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • θ = 0 • θ = 0.1 • θ = 0.25 • θ = 0.5 • θ = 1 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • θ = 0 • θ = 0.1 • θ = 0.25 • θ = 0.5 • θ = 1 1 10 -2 10 -4 10 -6
• • • • • • • • • • • • • • • • • • • • • • • • • • • ℓ = 100 • ℓ = 1000 • ℓ = 10000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • ℓ = 100 • ℓ = 1000 • ℓ = 10000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • N = 100 • N = 1000 • N = 5000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • N = 100 • N = 1000 • N = 5000
• • • • • • • • • • • • • • • • • • • • • • • • • • • m = 1 • m = 10 • m = 50 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • m = 1 • m = 10 • m = 50 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • m = 1, θ = 0.1 • m = 10, θ = 0.01 • m = 50, θ = 0.002 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • m = 1, θ = 0.1 • m = 10, θ = 0.01 • m = 50, θ = 0.002 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • n = 40 • n = 60 • n = 100 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • n = 40 • n = 60 • n = 100
• • • • • • • • • • • • • • • • • • • • • ℓ s = 1 • ℓ s = 9 • ℓ s = 99 • ℓ s =
• • • • • • • • • • • • • • • N = 500 • N = 5000 • N =
• • • • • • • • • • • • • • • • n = 10 • n = 20 • n =
σ • • • • • • • • • • • • • • • • m = 1 • m = 5 • m = 50 0.2 0.4 0.6 0.8 1.0 θ 0.02 0.04 0.06 0.08 0.10 σ • • • • • • • • • • • • • • • • c = 1 • c = 1.2 • c
• • • • • • • • • • • • • • • • ℓ = 10 3 • ℓ = 10 4 • ℓ =
σ • • • • • • • • • • • • • • • • • • • • • s d = 0.1 • s d = 0.01 • s d = 0.001 • s d = 0.
σ • • • • • • • • • • • • • • • • • • • • • • • • • • • • • R = 0 • R = 0.5 • R = 2 • R = 5 • R = 10 0.2 0.4 0.6 0.8 1.0 θ 0.05 0.10 0.15 σ • • • • • • • • • • • • • • • • U = 0.1 • U = 0.5 • U = 1 0.2 0.
• • • • • • • • • • • • • • • • • • • •                      •, σ = 0 •, σ = 0.1 •, σ = 1 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g θ = 0 • • • • • • • • • • • • • • • • • • • • •                      1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g θ = 0.1 • • • • • • • • • • • • • • • • • • • • •                      1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g θ = 0.5 Figure S2. Terms W g = exp [-n α=1 g α 2 / (2V s )], W Vg = exp [-n α=1 V g, α / (2V s )]
representing the effect of departures of mean phenotypes from the optimum (W g , circles, solid lines) and the effect of genetic variance (W Vg , squares, dashed lines) on the mutation load (L ≈ 1 -W g W Vg , see Supplementary File S2), for different values of s d and θ. Parameter values are as in Figures 3 andS1. Note that the lines simply connect simulation results and do not correspond to analytical approximations.

Mutational bias causes an increase in load through W g .

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •                                              • θ = 0 • θ = 0.1 • θ = 0.25 • θ = 0.5 • θ = 1 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •                                              •, θ = 0 •, θ = 0.1 •, θ = 0.25 •, θ = 0.5 •, θ = 1 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 1 Figure S3
. Same as Figure 4 in the main text, showing

W g = exp [-n α=1 g α 2 / (2V s )]
(circles, solid lines) and

W Vg = exp [-n α=1 V g, α / (2V s )] (squares, dashed lines). • • • • • • • • • • • • • • • • • • • • • • • • • • • • ℓ = 100 • ℓ = 1000 • ℓ = 10000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • ℓ = 100 • ℓ = 1000 • ℓ = 10000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • N = 100 • N = 1000 • N = 5000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • N = 100 • N = 1000 • N = 5000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 L σ = 1 Figure S4
. Same as Figure 5 in the main text, with θ = 0.5. (circles, solid lines) and The fitness of each individual is then computed according to equation 3 in the main text, where ω 2 is fixed to 10; however, as noted above, the values of ω 2 and V e should have little effect on the results (for given values of s d and θ), since V s = ω 2 + V e may be considered as a scaling factor.

• • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, ℓ = 100 •, ℓ = 1000 •, ℓ = 10000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, ℓ = 100 •, ℓ = 1000 •, ℓ = 10000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, N = 100 •, N = 1000 •, N = 5000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, N = 100 •, N = 1000 •, N = 5000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 1
W Vg = exp [-n α=1 V g, α / (2V s )] (squares, dashed lines). • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, ℓ = 100 •, ℓ = 1000 •, ℓ = 10000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, ℓ = 100 •, ℓ = 1000 •, ℓ = 10000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, N = 100 •, N = 1000 •, N = 5000 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 1 Figure S6. Same as Figure S4, showing W g = exp [-n α=1 g α 2 / (2V s )] (circles, solid lines) and W Vg = exp [-n α=1 V g, α / (2V s )] (squares, dashed lines). • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, m = 1 •, m = 10 •, m = 50 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, m = 1 •, m = 10 •, m = 50 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, m = 1, θ = 0.1 •, m = 10, θ = 0.01 •, m = 50, θ = 0.002 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, m = 1, θ = 0.1 •, m = 10, θ = 0.01 •, m = 50, θ = 0.002 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 1 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, n = 40 •, n = 60 •, n = 100 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 0 • • • • • • • • • • • • • • • • • • • • • • • • • • •                            •, n = 40 •, n = 60 •, n = 100 1 10 -2 10 -4 10 -6 s d 0.2 0.4 0.6 0.8 1.0 W g _ , W V g σ = 1
Investment in sexual reproduction s is coded by s loci, which are also equally spaced along the chromosome. These loci are multiallelic, investment in sex being given by:

s = s init + s i=1 g si + e s (1)
where s init is the initial investment in sex, g si the effect of locus i on s, and e s an environmental component drawn from a centered Gaussian distribution with variance V e, s . If the value of s obtained from equation 1 is lower than 0 or higher than 1, it is then set to 0 or 1 (respectively). During a number of preliminary generations, g si is fixed to zero for all loci affecting the rate of sex. Then, mutation occurs at a rate U s per generation on the whole set of s loci. When a mutation occurs at locus i, a quantity drawn from a centered Gaussian distribution with variance a 2 s is added to g si .

For each of the N individuals of the next generation, a maternal parent is sampled with a probability proportional to W 1 -s + s c , where W is its fitness and s its investment into sex. With probability 1-σ (where σ is given by equation 13 in the main text), the new individual is produced asexually and carries the same genotype as its mother, except for new mutations (the number of mutations on biallelic loci affecting the traits under stabilizing selection is drawn from a Poisson distribution with parameter U ). With the complementary probability, the new individual is produced sexually; in this case a paternal parent is sampled with a probability proportional to W s, and a recombinant offspring is produced (the number of cross-overs occurring at meiosis is sampled from a Poisson distribution with parameter R, and the position of each cross-over is drawn from a uniform distribution along the chromosome). Every 100 generations, the mean investment in sex, mean rate of sex, mean fitness, mean trait values, genetic variances and covariances among traits and some higher moments of phenotypic distributions are recorded by the program.

Genetic associations and decomposition of the genetic variance. Using the notation of Barton and Turelli (1991) and [START_REF] Kirkpatrick | General models of multilocus evolution[END_REF], we denote p i = E [X i ] the frequency of allele 1 at locus i and define ζ i as:

ζ i = X i -p i . (3)
Furthermore, products of ζ i variables are denoted:

ζ U = i∈U ζ i (4)
where U represents a set of loci. For example, for U = {i, j}, we have:

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

D U = E [ζ U ] (6)
In particular, D ij is the linkage disequilibrium between loci i and j. As we will see, associations involving repeated indices (such as D iij = E (X i -p i ) 2 (X j -p j ) ) sometimes appear in the computations. Using the fact that X i 2 = X i (since X i equals 0 or 1), repeated indices can be eliminated using the relation:

D Uii = p i q i D U + (1 -2p i ) D Ui (7)
with q i = 1 -p i (e.g., equation 5 in [START_REF] Kirkpatrick | General models of multilocus evolution[END_REF]. In particular (and because D j = E [X j -p j ] = 0), we have

D iij = (1 -2p i ) D ij . Similarly, D iijj = p i q i p j q j + (1 -2p i ) (1 -2p j ) D ij , while D ii = p i q i .
The genetic variance for trait α in the population is given by:

V g, α = E (g α -g α ) 2 (8)
From equations 5 and 6 in the main text:

g α = E i=1 r αi X i = i=1 r αi p i (9)
so that g α = i r αi p i . Using the definitions above, we have:

V g, α = E   i=1 r αi (X i -p i ) 2   = E   i=1 r αi ζ i 2   = E i,j r αi r αj ζ i ζ j (10)
where the last sum is over all i and j (including i = j). Using equations 4 and 6, one obtains:

V g, α = i=1 r αi 2 p i q i + i =j r αi r αj D ij . (11)
In the following, we assume that linkage disequilibria remain negligible, so that

V g, α ≈ i=1 r αi 2 p i q i .
Neglecting the effects of selection on mean trait values and genetic variance.

Simple approximations for g α and V g, α are obtained for the regime where s d is so low that selection has negligible effects on p i and p i q i , compared with the effects of mutation and drift. Because drift does not change expected allele frequencies, the change in p i over one generation is given by (neglecting the effect of selection):

p i t+1 = u + (1 -2u) p i t . (12)
yielding p i = 1/2 at equilibrium. Using the change in phenotypic basis given by equations 9 and 10 in the main text, we have:

1 2V s n α=1 g α 2 = 1 2V s g 1 2 (13)
since only the first trait in the new basis (along which the mutational bias occurs)

should differ from zero, on average, at equilibrium. Equation 9 and equation 11 in the main text then yield:

1 2V s n α=1 g α 2 = 1 4 b1 2 ( 14
)
where b1 is the scaled magnitude of mutational bias.

Neglecting the effects of selection yields the following recursion for p i q i :

p i q i t+1 ≈ 1 - 1 N [u + (1 -4u) p i q i t ] (15)
so that p i q i ≈ N u/ (1 + 4N u) at equilibrium (assuming large N and small u). Noting that α i r αi 2 = 2V s s d (see equation 7 in the main text), one obtains:

1 2V s n α=1 V g, α ≈ s d N U 1 + 4N u . (16)
Equations 2, 14 and 16 yield the following approximation for the load (assuming

g α 2 ≈ g α 2 ): L ≈ 1 -exp -s d N U 1 + 4N u - 1 4 b1 2 , (17)
equivalent to equation 8 in Roze and Blanckaert (2014) in the absence of mutational bias ( b1 = 0). Comparisons with individual-based simulations show that equation 17 does indeed provide correct predictions when s d is very low (see Figures 234in the main text). As s d increases, g α and V g, α depart more and more from equations 14 and 16; however, simulations indicate that equation 16 stays valid over a wider range of values of s than equation 14, in agreement with previous observations that selection may have significant effects on mean trait values even when p i q i at each locus is mainly controlled by mutation and drift (Campbell, 1984;Barton, 1989;Charlesworth, 2013a). Based on this, we can derive a better approximation for low s d by taking the effect of selection on g α into account, but still neglect the effect of selection on V g, α ,

as shown in the next subsection.

Effect of selection on mean trait values in the low s d regime. From equation 12, we have:

g α t+1 = u i r αi + (1 -2u) ( g α t + ∆ sel g α t ) (18)
where ∆ sel g α is the change in g α due to selection, given by:

∆ sel g α = E W g W g α -g α = E W g W (g α -g α ) . (19)
Assuming weak selection (V g, α /V s small, of order ), we have from equation 4 in the main text (to the first order in ):

W g W g,max ≈ 1 - 1 2V s n α=1 g α 2 = 1 - 1 2V s n α=1 g α 2 + 2 (g α -g α ) g α + (g α -g α ) 2 , (20) yielding:
W W g,max ≈ 1 - 1 2V s n α=1 g α 2 + V g, α , (21)
and thus:

W g W ≈ 1 - 1 V s n α=1 g α (g α -g α ) - 1 2V s n α=1 (g α -g α ) 2 -V g, α . (22)
From equations 19 and 22, one obtains:

∆ sel g α ≈ - 1 V s n β=1 g β C g, αβ - 1 2V s n β=1 M g, αββ (23)
where M g, αββ is the third moment E (g α -g α ) (g β -g β ) 2 . Assuming that the distribution of phenotypes in the population stays close to a Gaussian distribution, M g, αββ should be close to zero. Furthermore, assuming that fluctuations in g β and C g, αβ due to drift remain small, g β C g, αβ ≈ g β C g, αβ . Measuring traits in the phenotypic basis defined by equations 9 and 10 in the main text (so that C g, αβ = 0 for α = β), one obtains the following expression for g 1 at equilibrium (from equations 18 and 23):

g 1 ≈ U b 1 1 -(1 -2u) 1 -V g, 1 Vs . ( 24
)
Neglecting the effects of selection on V g, 1 ≈ i=1 (r 1i ) 2 p i q i and noting that

i=1 (r 1i ) 2 = m n (a 2 + m b 2 ) = 2V s s d 1 n [1 + θ (m - 1
)] (using equation 9 in the main text), we have:

V g, 1 ≈ 2V s s d 1 n [1 + θ (m -1)] N U 1 + 4N u . ( 25
)
From equations 24 and 25, one obtains (assuming u and s d are small):

g 1 ≈ b 1 2 1 + s d n [1 + θ (m -1)] N 1+4N u (26)
Equations 2, 16 and 26 yield (assuming g α 2 ≈ g α 2 ):

L ≈ 1 -exp    -s d N U 1 + 4N u - b1 2 4 1 + s d n [1 + θ (m -1)] N 1+4N u 2    . (27)
Comparisons with individual-based simulations confirm that equation 27 provides better predictions than equation 17 in the case of sexual populations, as long as s d is sufficiently small (see Figures 345in the main text). Equation 27 fails when s d is not very small, however, as selection affects genetic variances at equilibrium. Unfortunately, we could not obtain any simple expression for the genetic variance (and mean fitness) in this regime for arbitrary θ, although an approximation can be obtained for θ = 1, as shown in the next subsection.

Effect of selection on genetic variance and approximations for the mutationselection regime. Neglecting linkage disequilibria, genetic variances can be expressed in terms of the genetic diversities p i q i at the different loci (equation 11). Extending equation 15 to include selection yields:

p i q i t+1 ≈ 1 - 1 N u + (1 -4u) p sel i q sel i t . (28)
Furthermore, noting that

p i q i = D ii = E [ζ ii ],
we have, to the first order in :

p sel i q sel i = E W W ζ ii . (29)
Decomposing g α , g α and V g, α as sums over loci (using equations 9 and 11) and introducing the centered variables ζ i = X i -p i , we have from equation 22:

W g W = 1 + i=1 a i ζ i + i,j a ij (ζ ij -D ij ) (30)
with a i = -1

Vs n α=1 g α r αi and a ij = -n α=1 (r αi r αj ) / (2V s ), both of order (Barton and Turelli, 1991). Using equations 29 and 30, and neglecting linkage disequilibria, one obtains:

p sel i q sel i = p i q i - 1 V s n α=1 r αi z α (1 -2p i ) p i q i - 1 2V s n α=1 r αi 2 (1 -2p i ) 2 p i q i . (31)
Equations 28 and 31 lead to the following recursion for the genetic variance:

V g, α t+1 ≈ 1 - 1 N 2V s s d U n + (1 -4u) V g, α t - 1 V s n β=1 z β C ααβ t - 1 2V s i=1 r αi 2 n β=1 r βi 2 (1 -2p i ) 2 p i q i t . (32)
It is not possible to derive an expression for V g, α at equilibrium under mutation, selection and drift from equation 32 -one may assume that z β C ααβ ≈ z β C ααβ and that C ααβ ≈ 0, and assume that p i q i is small at each locus so that (1

-2p i ) 2 ≈ 1
(and the second sum on the second line of equation 32 becomes β V g, β ), but the resulting approximation does not work well when s d is small (as p i q i may not be small), nor when s d is large and in the presence of mutational bias (as C ααβ = 0, results not shown). Neglecting drift, and in the absence of mutational bias (θ = 0), the change in p i q i over one generation is (from equations 28 and 31):

∆ (p i q i ) ≈ u (1 -2p i ) 2 - 1 2V s n α=1 r αi 2 (1 -2p i ) 2 p i q i (33) so that either p i = 1/2 or p i q i = 2V s u/ ( n α=1 r αi 2 ) at equilibrium. When s d u,
most loci should be at the second equilibrium, in which case V g, α ≈ 2V s U/n, and

L ≈ 1 -e -U .
Another approximation can be obtained for the case where θ = 1 (no variance of mutational effects) and when drift is negligible. Indeed, in this case alleles 1 are deleterious, and p i should thus be small at equilibrium. To the first order in p i , equations 28 and 31 give for the change in p i over one generation:

∆p i ≈ u - 1 V s n α=1 r αi g α p i - 1 2V s n α=1 r αi 2 p i (34)
which may also be written as (using the change in phenotypic basis given by equations 9 and 10 in the main text):

∆p i ≈ u - 1 V s b 1 g 1 p i -s d p i . (35)
From this, the change in g 1 ≈ i r 1i p i is:

∆g 1 ≈ U b 1 - 1 V s b 1 g 1 2 -s d g 1 (36)
yielding, at equilibrium:

g 1 √ 2V s ≈ 8 b1 2 U + s d 2 -s d 4 b1 . (37)
Finally, we have α V g, α ≈ i α r αi 2 p i , which is also 2V s s d p i (as p i should be the same at all loci when θ = 1). Noting that g 1 ≈ b 1 p i , we thus have:

n α=1 V g, α 2V s ≈ s d b1 g 1 √ 2V s . (38)
Equations 2, 37 and 38 finally lead to:

L ≈ 1 -exp   - 4 b1 2 U + s d 8 b1 2 U + s d 2 -s d 8 b1 2   (39)
or, in terms of s d , ρ and U :

L ≈ 1 -exp - 4ρU -s d + s d (8ρU + s d ) 8ρ . (40)
Simulations confirm that equation 40 provides accurate predictions for θ = 1 (in sexual populations), when s d is sufficiently high (see Figure 4 in the main text).

FILE S3: QLE MODEL FOR THE EVOLUTION OF SEX

We derive here expressions for the change in mean rate of sex in the limit of an infinite population, using a quasi-linkage equilibrium (QLE) argument. For this, we use Turelli and Barton's (1990) method (see also [START_REF] Barton | A general model for the evolution of recombination[END_REF] to express the effect of selection on genetic associations in terms of partial derivatives of ln W with respect to mean trait values and genetic variances/covariances. Note that the derivations given below are in principle valid for any number of possible alleles at each locus (not necessarily biallelic loci) and any fitness function (not necessarily Gaussian), as long as the distribution of phenotypes affecting fecundity stays approximately Gaussian.

Definitions. Extending the notation of Turelli and Barton (1990) to multiple traits, we define the centered variable ζ αj as:

ζ αj = g αj -g αj (1)
(where again g αj is the effect of the allele present at locus j on trait α in a given individual, and g αj its average over all individuals). Genetic associations are defined as

C U = E [ζ U ] (2)
where E stands for the average over all individuals, and with ζ U = x ζ x , each x bearing two elements, the trait α and the locus j. For example, C αj αj = E (g αj -g αj ) 2 while

C αj αk βk = E [(g αj -g αj ) (g αk -g αk ) (g βk -g βk )]
. Using these definitions, the genetic variance for trait α can be written as:

V g, α = E   j ζ αj 2   = j=1 C αj αj + j =k C αj αk = V 0 g, α + D αα (3)
where V 0 g, α = j C αj αj is the "genic variance" for trait α (genetic variance in a population with the same allele frequencies, at linkage equilibrium), and

D αα = j =k C αj αk
is the effect of linkage disequilibria on the variance. Similarly, the genetic covariance between traits α and β can be decomposed as:

C g, αβ = j=1 C αj βj + j =k C αj βk = C 0 g, αβ + D αβ . (4)
As explained in the main text, we assume that investment in sexual reproduction s is also a polygenic trait with independent genetic and environmental contributions:

s = s + g s + e s (5)
where e s is sampled in a Gaussian distribution with mean 0 and variance V e, s , and assuming additive effects of loci affecting s:

g s = j g sj (6)
where g sj is the effect of the allele at locus j on investment in sex (we assume that loci affecting the rate of sex do not affect the other traits). Assuming that the variance in s in the population is sufficiently small (g s , g e small, of order η), the rate of sex σ = s/ [c (1 -s) + s] of an individual can also be decomposed into an additive genetic and an environmental component:

σ = σ + g σ + e σ (7)
with (to leading order in η):

σ ≈ s c (1 -s) + s , g σ ≈ c [c (1 -s) + s] 2 g s , V g, σ ≈ c 2 [c (1 -s) + s] 4 V g, s . (8)
The expected change in σ over one generation (denoted ∆σ) corresponds to the change in g σ in the parental generation due to differences in fecundities among individuals and to the cost of sex. In the following, we derive deterministic approximations for ∆σ, assuming that phenotypic traits affecting fecundity are normally distributed (this implies that the number of loci affecting each of these traits is sufficiently large, each locus having a sufficiently small effect on the trait). We will also use a quasi-linkage equilibrium approximation, assuming that rates of sex and recombination are not too small, so that genetic associations between loci are small and equilibrate fast relative to change in allele frequencies. Finally, we will assume that the genetic variance in the rate of sex in the population (V g, σ ) is small (however, we do not make any assumption on the number of loci affecting the rate of sex). For this, we will decompose a generation into two steps: the first ("selection") corresponds to the differential reproduction of individuals due to differences in fecundity (according to the values of their phenotypes z 1 , . . . z n ), while the second ("reproduction") corresponds to the effect of the cost of sex and of sexual recombination (strictly, this second step also involves selection when c > 1, since individuals investing more in sex are disfavored). In the next sections, we derive expressions for changes in mean breeding values during these two steps.

Effect of selection on mean breeding values. The effect of selection on g α can be written as:

∆ sel g α = E W g W (g α -g α ) (9)
where W g is the mean fecundity of individuals with breeding values g = (g 1 , g 2 , . . .)

and W the mean fecundity of the whole population. Following [START_REF] Barton | A general model for the evolution of recombination[END_REF], we assume that selection is weak and approximate W g /W by:

W g W ≈ 1 + α (g α -g α ) ∂ ln W ∂z α + α≤β [(g α -g α ) (g β -g β ) -C g, αβ ] ∂ ln W ∂C g, αβ . (10)
(see Appendix A), where the last sum includes the terms for α = β, which involve partial derivatives of ln W with respect to V g, α . From equations 9 and 10, and assuming a Gaussian distribution of breeding values, we recover the classic expression: [START_REF] Lande | Quantitative genetic analysis of multivariate evolution, applied to brain:body size allometry[END_REF]. The change in g σ is obtained similarly:

∆ sel g α = β ∂ ln W ∂z β C g, αβ (11)
∆ sel g σ = E W g W (g σ -g σ ) . (12)
However, we can no longer assume that the joint distribution of investment into sex σ and of traits affecting fecundity is multivariate normal (indeed, genetic variances and covariances may differ between subgroups of individuals differing in their values of g σ , due to the effect of sexual recombination on genetic associations). From equations 10 and 12, one obtains:

∆ sel g σ ≈ α ∂ ln W ∂z α C g, σα + α≤β ∂ ln W ∂C g, αβ M g, σαβ (13)
where M g, σαβ is the moment E [(g σ -g σ ) (g α -g α ) (g β -g β )]. Equation 13 is equivalent to Charlesworth (1993)'s decomposition of the selection gradient for a recombination modifier allele into two terms (equation A10 in Charlesworth, 1993, see also Appendix 2 of [START_REF] Barton | A general model for the evolution of recombination[END_REF]. The first part of equation 13 shows that under directional selection acting on trait α, a covariance between g α and g σ generates indirect selection on σ (this is equivalent to the term in δz in Charlesworth, 1993). The second part of equation 13 (equivalent to the term in δV g in Charlesworth, 1993) corresponds to indirect selection on σ due to different genetic variances and covariances for selected traits among subgroups of individuals with different rates of sex. For example, ∂ ln W /∂V g, α < 0 under stabilizing selection acting on a single trait α (the immediate effect of increasing genetic variance is to decrease mean fitness) and in this situation we also expect that higher rates of sex tend to increase genetic variance, so that E (g σ -g σ ) (g α -g α ) 2 > 0, and the second term of equation 13 selects against sex.

This term is equivalent to the "short-term effect" in models for the evolution of sex (or recombination) with epistasis (e.g., Agrawal, 2006). Now, if the population mean phenotype z = (z 1 , . . . z n ) is displaced from the optimum, the higher genetic variance associated with sex will increase the efficiency of selection, generating associations between higher values of g σ and values of g α closer to the optimum, that in turn favor sex. This effect is represented by the first term of equation 13 and corresponds to the "long term effect" (favoring sex due to an increased efficiency of selection).

Selection gradients ∂ ln W /∂z α and ∂ ln W /∂C g, αβ take simple forms in the case of a fully isotropic model with Gaussian stabilizing selection:

W = exp -α (z α -θ α ) 2 2ω 2 (14)
where ω 2 represents the strength of selection (the same for all traits), and where the phenotypic optimum is located at θ = (θ 1 , θ 2 , . . . θ n ). A general expression for mean fitness under Gaussian stabilizing selection (and when the maximal possible fitness is 1, as implied by equation 14) is given by: [START_REF] Gomulkiewicz | Demographic and genetic constraints on evolution[END_REF], where det(A) is the determinant of matrix A, T stands for matrix/vector transpose, S is a matrix determining the pattern of multivariate stabilizing selection, and P is the phenotypic variance-covariance matrix.

W = det (S + P) -1 S exp - 1 2 (z -θ) T (S + P) -1 (z -θ) (15)
When fitness is given by equation 14, S = ω 2 I (where I is the identity matrix). Furthermore, our assumption of independent, identically distributed environmental effects yields P = G + V e I, where G is the genetic variance-covariance matrix (whose diagonal elements are genetic variances V g, α , and off-diagonal elements genetic covariances

C g, αβ )
. It is always possible to find an orthonormal basis in which the G matrix is diagonal, that is, to define new phenotypic traits as linear combinations of the "true" phenotypic traits so that the new traits are independent, and fitness is still given by equation 14. Assuming that phenotypes are measured in this new basis, we show in Appendix B that:

∂ ln W ∂z α = - z α -θ α V g, α + V s ( 16
)
where V s = ω 2 + V e , while:

∂ ln W ∂V g, α = - 1 2 (V g, α + V s ) + 1 2 ∂ ln W ∂z α 2 (17) ∂ ln W ∂C g, αβ = ∂ ln W ∂z α ∂ ln W ∂z β , (18)
for α = β.

Change in mean rate of sex during reproduction. To compute the change in g σ during reproduction (due to the cost of sex), we first compute the change in g s .

We have:

∆ rep g s = E c (1 -s) c (1 -s ) + s (g s -g s ) + E s ♀ s ♂ s [c (1 -s ) + s ] (g s,♀ -g s ) + (g s,♂ -g s ) 2 (19)
where the primes denote averages among individuals after selection (that is, weighting each individual by its relative fecundity), and where the average on the second line is over all possible pairs of female and male parents, s ♀ and s ♂ being the investments in sex of these parents, and g s,♀ , g s,♂ their value of g s . Equation 19 can be understood as follows. The term on the first line is the proportion of asexually produced offspring -which is (1 -s ) / 1 -s + s c -multiplied by the change in the mean value of g s among those offspring relative to the parents: for this, each parent is weighted by its relative contribution to the pool of asexually produced offspring, which is (1 -s) / (1 -s ). The term on the second line is the proportion of sexually produced offspring -which is s c / 1 -s + s c -multiplied by the change in the mean value of g s among those offspring. On average, the mean value of g s among the offspring of a given female and male is (g s,♀ + g s,♂ ) /2, where g s,♀ and g s,♂ are the values of g s in the parents. Furthermore, the relative contributions of these parents to the pool of sexually produced offspring are s ♀ /s and s ♂ /s . Replacing s by s + g s -g s + g s + e s (and similarly for g s,♀ , g s,♂ ) in equation 19, and using s = s + g s finally yields:

∆ rep g s = - c -1 c (1 -s ) + s V g, s . (20)
Equation 20 represents the effect of direct selection against sex (whenever c > 1), and is equivalent to the expression derived in Roze (2014) in the case of a single biallelic sex modifier locus. Strictly, s and V g, s in equation 20 are the mean and genetic variance for investment in sex after selection (weighting each individual by its relative fecundity). However, taking into account the change in s and V g, s due to selection would introduce terms in V g, s 2 in equation 20; neglecting those terms, s and V g, s in equation 20 can thus be replaced by their values s and V g, s at the start of the generation (before selection). From equations 8 and 20, one then obtains:

∆ rep g σ ≈ - c -1 1 + (c -1) σ V g, σ . (21)
Assuming no mutational bias on σ, the change in the mean rate of sex over one generation is given by:

∆σ = ∆ sel g σ + ∆ rep g σ . (22)
In the following, we derive approximate expressions for the moments C g, sα and M g, sαβ that appear in the expression of ∆ sel g σ (equation 13). However before that, we will compute an expression for the contribution of linkage disequilibria to the genetic variances and covariances between traits affecting fecundity (D αβ = j =k C αj βk , equations 3 and 4), at quasi-linkage equilibrium.

Genetic associations between selected loci. Neglecting genetic variance in the rate of sex, C αj βk at the next generation is given by:

C αj βk = (1 -σ) C αj βk + σ (1 -r jk ) C αj βk (23)
where the double prime denotes variables measured at the next generation (after reproduction), and r jk is the recombination rate between loci j and k. The first term of equation 23 is the proportion of asexually produced offspring, multiplied by the genetic association among those offspring, which is the same as among parents. The second term is the proportion of sexually produced offspring, in which C αj βk is decreased by a factor 1 -r jk due to recombination. Equation 23 can be written under the simpler form:

C αj βk = (1 -ρ jk ) C αj βk (24)
where ρ jk = σ r jk is the "effective" recombination rate between loci j and k.

The effect of selection on C αj βk can be computed as follows (Turelli and Barton, 1990;[START_REF] Barton | A general model for the evolution of recombination[END_REF]. We have:

C αj βk = E [(g αj -g αj ) (g βk -g βk )] (25)
where again the prime denotes averages after selection (weighting each individual by its relative fecundity). Equation 25can also be written:

C αj βk = E [(g αj -g αj -∆ sel g αj ) (g βk -g βk -∆ sel g βk )] . (26)
In the following, we use the notation C sel U for genetic associations measured after selection, but using as "reference values" (the g αj in equation 1) allelic averages before selection: in particular, C sel αj βk = E [(g αj -g αj ) (g βk -g βk )]. Expanding equation 26 and noting that C sel αj = E [g αj -g αj ] = ∆ sel g αj , one obtains:

C αj βk = C sel αj βk -(∆ sel g αj ) (∆ sel g βk ) . (27)
Furthermore, we have:

C sel αj βk = E W g W ζ αj βk , ∆ sel g αj = E W g W ζ αj . (28)
From equation 10, and noting that g α -g α = j ζ αj , while (g α -g α ) (g β -g β )-C g, αβ = j,k (ζ αj βk -C αj βk ) (where the last sum is over all pairs of loci j and k including j = k), we have:

W g W ≈ 1 + α ∂ ln W ∂z α j ζ αj + α≤β ∂ ln W ∂C g, αβ j,k (ζ αj βk -C αj βk ) . (29)
Equations 28 and 29 yield:

C sel αj βk = C αj βk + γ ∂ ln W ∂z γ i C γi αj βk + γ≤δ ∂ ln W ∂C g,γδ h,i (C γh δi αj βk -C γh δi C αj βk ) . (30)
Equation 30 shows that C αj αk is affected by higher-order associations (involving 3 or 4 loci). These associations are in turn generated by the effect of selection, and eroded by recombination. In the following we assume that selection is sufficiently weak relative to recombination, so that between-locus associations remain small (Turelli and Barton, 1990). Ignoring terms involving between-locus associations in the sums on the righthand-side of equation 30, only the terms for h = j, i = k and h = k, i = j in the last sum remain, giving:

C sel αj αk ≈ C αj αk + γ≤δ ∂ ln W ∂C g,γδ (C αj γj C βk δk + C αj δj C βk γk ) . (31)
Equations 28 and 29 also yield (neglecting between-locus associations):

∆ sel g αj = γ ∂ ln W ∂z γ i C αj γj + γ≤δ ∂ ln W ∂C g,γδ h,i C αj γj δj . (32)
C αj βk at QLE is obtained by solving C αj βk = C αj βk . From equations 24, 27, 31 and 32, this yields:

C αj βk = 1 ρ jk -1 ∆ sel C αj βk (33) with ∆ sel C αj βk = γ,δ (1 + I γδ ) ∂ ln W ∂C g,γδ C αj γj C βk δk -(∆ sel g αj ) (∆ sel g αk ) . (34)
where I γδ equals 1 if γ = δ, and 0 otherwise. Summing over all loci, one obtains for

D αβ = j =k C αj βk : D αβ ≈ 1 ρ h -1 ∆ sel D αβ (35)
where ρ h is the harmonic mean of ρ jk over all pairs of loci affecting fecundity, and with

∆ sel D αβ ≈ γ,δ (1 + I γδ ) ∂ ln W ∂C g,γδ - ∂ ln W ∂z γ ∂ ln W ∂z δ C 0 g, αγ C 0 g, βδ . (36
D αα = j =k C αj αk ≈ 1 ρ h -1 2 ∂ ln W ∂V g, α - ∂ ln W ∂z α 2 V g, α 2 . ( 37
)
When the fitness function is given by equation 14, this simplifies to (using equations 16 and 17):

D αα ≈ - 1 ρ h -1 V g, α 2 V g, α + V s (38)
corresponding to the result obtained by [START_REF] Bulmer | The Mathematical Theory of Quantitative Genetics[END_REF] under the assumption of exchangeable loci (equations A3c and A6a in Charlesworth, 1993).

Indirect selection for sex: "short-term effect". As discussed earlier, the "shortterm effect" is represented by the term on the second line of equation 13, which depends on moments M g, σαβ = E [(g σ -g σ ) (g α -g α ) (g β -g β )] (for all traits α, β affecting fecundity). From equation 8, we have:

M g, σαβ ≈ c [c (1 -s) + s] 2 M g, sαβ (39)
with

M g, sαβ = E [(g s -g s ) (g α -g α ) (g β -g β )]
. Furthermore, using our definition of genetic associations, M g, sαβ can be decomposed as:

M g, sαβ = i,j,k C si αj βk (40)
where the sum is over all loci i affecting investment sex, and over all pairs of loci j and k affecting traits α, β. A QLE approximation for C si αj βk can be obtained as follows.

At the next generation, we have:

C si αj βk = E [(g si -g si ) (g αj -g αj ) (g βk -g βk )] (41)
where again the double primes denote averages over individuals of the next generation (after reproduction). Equation 41 can also be written:

C si αj βk = E (g si -g si -∆ rep g si ) (g αj -g αj -∆ rep g αj ) × (g βk -g βk -∆ rep g βk ) (42)
where g si , g αj , g βk are the averages of g si , g αj , g βk among selected parents (weighting each parent by its relative fecundity), and ∆ rep g si = g si -g si the change in g si during reproduction, due to the cost of sex (and similarly for ∆ rep g αj and ∆ rep g βk ).

In

C si αj βk = C rep si αj βk -(∆ rep g si ) C rep αj βk -(∆ rep g αj ) C rep si βk -(∆ rep g βk ) C rep si αj + 2 (∆ rep g si ) (∆ rep g αj ) (∆ rep g βk ) . (43)
The change in g αj during reproduction is generated by the cost of sex and by genetic associations between locus j and loci affecting investment in sex, and is thus proportional to V g, s (the same is true for ∆ rep g βk ). Furthermore, the sum over all i and j of C si βk is the genetic covariance between trait β and investment in sex s, which is also proportional to V g, s . As a consequence, the last three terms of equation 43 will generate terms in O(V g, s 2 ), and will thus be ignored, so that:

C si αj βk ≈ C rep si αj βk -(∆ rep g si ) C rep αj βk . (44)
Using a similar reasoning as when deriving equation 19, an expression for C rep si αj βk is given by (using s ≈ s, and for j = k):

C rep si αj βk = E c (1 -s) c (1 -s) + s ζ si αj βk + E s ♀ s ♂ s [c (1 -s) + s] r ijk,∅ ζ si αj βk,♀ + r ∅,ijk ζ si αj βk,♂ + r i,jk ζ si,♀ ζ αj βk,♂ + r jk,i ζ αj βk,♀ ζ si,♂ + r ij,k ζ si αj,♀ ζ βk,♂ + r k,ij ζ βk,♀ ζ si αj,♂ + r ik,j ζ si βk,♀ ζ αj,♂ + r j,ik ζ αj,♀ ζ si βk,♂ . (45)
In equation 45, r S,T is the probability that a meiotic product inherits the set S of loci from the maternal genome, and the set T of loci from the paternal genome, while

C rep si αj βk = 1 - s c (1 -s) + s (1 -r ijk,∅ -r ∅,ijk ) C si αj βk + 1 c (1 -s) + s -(c -r ijk,∅ -r ∅,ijk ) h C sh si αj βk + (r i,jk + r jk,i ) C αj βk h C sh si + l C sl αj βk h C sh si + (r ij,k + r k,ij ) C si αj h C sh βk + l C sl si αj h C sh βk + (r ik,j + r j,ik ) C si βk h C sh αj + l C sl si βk h C sh αj . (46)
Many of the terms of equation 46 may be neglected when V g, s is small, using the fact that sums over all loci of associations involving one or several "s" indices are proportional to V g, s : therefore, the terms on the last two lines of equation 46 and the last term of the third line will generate terms in V g, s 2 . Furthermore, we will neglect linkage disequilibria between loci affecting investment in sex, so that only the terms for h = i remain in the sums above, and equation 46 simplifies to:

C rep si αj βk ≈ (1 -ρ ijk ) C si αj βk - c -r ijk,∅ -r ∅,ijk c (1 -s) + s C si si αj βk + r i,jk + r jk,i c (1 -s) + s C si si C αj βk (47)
with ρ ijk = σ r ijk , r ijk = 1 -r ijk,∅ -r ∅,ijk being the probability that at least one recombination event occurs between loci i, j and k. Because C si si αj βk ≈ C si si C αj βk to leading order, equation 47 further simplifies to:

C rep si αj βk ≈ (1 -ρ ijk ) C si αj βk - c -1 + r jk c (1 -s) + s C si si C αj βk (48)
The term ∆ rep g si in equation 44 is given by:

∆ rep g si = E c (1 -s) c (1 -s) + s ζ si + E s ♀ s ♂ s [c (1 -s) + s] ζ si,♀ + ζ si,♂ 2 . ( 49
)
Neglecting linkage disequilibria between loci affecting the rate of sex, this yields:

∆ rep g si ≈ - c -1 c (1 -s) + s C si si . ( 50
)
From equation 24, C rep αj βk = (1 -ρ jk ) C αj βk to leading order (that is, neglecting genetic variation for the rate of sex), so that:

(∆ rep g si ) C rep αj βk ≈ - c -1 c (1 -s) + s (1 -ρ jk ) C si si C αj βk . (51)
Putting everything together, one obtains from equations 44, 48 and 51:

C si αj βk ≈ (1 -ρ ijk ) C si αj βk - c r jk [c (1 -s) + s] 2 C si si C αj βk . ( 52
)
The change in C si si due to selection may be neglected under our assumption that V g, s is small (as it would generate terms in V g, s 2 ). Furthermore, the effect of selection on C si αj βk can be neglected when selection is weak, as it involves higher-order associations between loci i, j, k and other loci, which are themselves generated by the effect of selection at these loci. Using these approximations, equation 52 becomes:

C si αj βk ≈ (1 -ρ ijk ) C si αj βk - c r jk [c (1 -s) + s] 2 C si si C αj βk (53)
giving at QLE:

C si αj βk ≈ - 1 ρ ijk c r jk C αj βk [c (1 -s) + s] 2 C si si . ( 54
)
From the results of the preceding subsection, C αj βk ≈ (∆ sel C αj βk ) /ρ jk (where ∆ sel C αj βk is given by equation 34), so that r jk C αj βk ≈ ∆ sel C αj βk /σ. Equation 54 thus simplifies to:

C si αj βk ≈ - 1 r ijk c s 2 (∆ sel C αj βk ) C si si . ( 55
)
The same reasoning as above can be used to compute C si αj βj , which is found to be negligible. Summing over all loci, one thus obtains:

M g, sαβ ≈ - 1 r h,1 c s 2 (∆ sel D αβ ) V g, s (56)
where r h,1 is the harmonic mean of r ijk over all triplets of loci involving one locus affecting sex and two loci affecting fecundity. From this, we have (using equations 8 and 39):

M g, σαβ ≈ - ∆ sel D αβ r h,1 σ 2 V g, σ . (57)
Indirect selection for sex due to the short-term effect (the second term of equation 13)

can thus be written approximately as:

- 1 r h,1 σ 2 α≤β ∂ ln W ∂C g, αβ ∆ sel D αβ V g, σ . (58)
We will see later that the term between parentheses can be expressed in terms of the effect of sex on the mean fitness of offspring.

Indirect selection for sex: "long-term effect". The long-term effect depends on genetic covariances between the rate of sex and traits affecting fecundity (C g, σα ). From equation 8, we have:

C g, σα ≈ c [c (1 -s) + s] 2 C g, sα (59)
with C g, sα = E [(g s -g s ) (g α -g α )], which can be decomposed as:

C g, sα = i,j C si αj . ( 60
)
Using the same approach as above, one obtains for the effect of reproduction on C si αj :

C si αj = C rep si αj -(∆ rep g si ) (∆ rep g αj ) ≈ C rep si αj (61)
since the term (∆ rep g si ) (∆ rep g αj ) will generate terms in V g, s 2 . Neglecting linkage disequilibria between loci affecting the rate of sex and other terms of order V g, s 2 , we have:

C rep si αj ≈ (1 -ρ ij ) C si αj - c -1 + r ij c (1 -s) + s C si si αj . ( 62
)
To the first order in V g, s ,

C
while the term (∆ sel g si ) (∆ sel g αj ) is of higher order in the strength of selection, and may thus be neglected. Finally, using the same method as in the previous subsection shows that associations C si αj βk γl (that appear on the second line of equation 67)

are proportional to C si si C αj βk γl . However, 3-locus associations C αj βk γl are of higher order in the strength of selection than 2-locus associations, and we will assume that the sum over all loci of these associations is negligible relative to the sum of pairwise associations C αj βk . This leaves us with the following recursion for C si αj :

C si αj ≈ (1 -ρ ij ) C si αj + β ∂ ln W ∂z β k C si αj βk . ( 68
)
At QLE, and using equation 55, we thus have:

C si αj ≈ - 1 ρ ij -1 c s 2 k 1 r ijk β ∂ ln W ∂z β (∆ sel C αj βk ) C si si , (69)
and summing over all loci:

C g, sα ≈ - 1 r h,2 σ - 1 r h,1 c s 2 β ∂ ln W ∂z β (∆ sel D αβ ) V g, s (70)
where r h,2 is the harmonic mean of r ij r ijk over all triplets of loci i, j and k, where i affects investment in sex while j and k affect fecundity. Equations 8 and 59 then yield:

C g, σα ≈ - 1 r h,2 σ - 1 r h,1 1 σ 2 β ∂ ln W ∂z β (∆ sel D αβ ) V g, σ (71)
and indirect selection for sex due to the long term effect (first term of equation 13) is thus approximately:

- 1 r h,2 σ - 1 r h,1 1 σ 2 α,β ∂ ln W ∂z α ∂ ln W ∂z β ∆ sel D αβ V g, σ . (72)
Note that the term in the first parenthesis of equation 72 is positive, 1/r h,1 becoming negligible compared with 1/ (r h,2 σ) as the rate of sex decreases. due to the effect of sexual recombination. From equation 33, we have ∆ sel D αβ ≈ j =k ρ jk C αj βk , so that:

∆ sex/asex C g, αβ ≈ - 1 σ ∆ sel D αβ . (75)
Therefore, from equation 73, the difference in ln W between sexually and asexually produced offspring is given by:

∆ 1 = ln W sex -ln W asex ≈ - 1 σ α≤β ∂ ln W ∂C g, αβ ∆ sel D αβ (76)
and indirect selection for sex due to the short-term effect (equation 58) thus becomes:

∆ 1 r h,1 σ V g, σ . (77)
Following [START_REF] Barton | A general model for the evolution of recombination[END_REF] and Charlesworth and Barton (1996), selection for sex due to the long-term effect can be expressed in terms of the effect of sex on the variance in log-fitness among offspring. From equation 10 we have, to leading order in selection gradients: (79)

ln W g -ln W ≈ α (g α -g α ) ∂ ln W ∂z α + α≤β [(g α -g α ) (g β -g β ) -C g, αβ ] ∂ ln W ∂C g, αβ (78)
Equation 79 is approximately equivalent to the first two lines of equation A3b in Charlesworth and Barton (1996), corresponding to the additive and epistatic components of the variance in log fitness (denoted hereafter V A and V AA ). Using equations 75 and 79, the sum appearing in the expression for the strength of selection for sex due to the long term effect (equation 72) can be expressed as -σ (V A, sex -V A, asex ), where

V A, sex and V A, asex are the additive components of the variance in log fitness (first term of equation 79) among offspring produced by sexual and asexual reproduction, respectively. Selection for sex due to the long term effect thus becomes:

1 r h,2 σ - 1 r h,1 ∆ 2 σ V g, σ (80)
with ∆ 2 = V A, sex -V A, asex . Assuming that epistasis is weak relative to directional selection, Charlesworth and Barton (1996) show that the effect of recombination on V AA may be neglected, in which case the long term effect can be expressed in terms of the effect of recombination on Var [ln W g ]. However, in situations where epistatic interactions may be of the same order of magnitude as directional selection (as in the present model), the additive component of Var [ln W g ] should be estimated, for example from the covariance between parents and offspring (e.g., Lynch and Walsh, 1998).

Indeed, under the assumption of a sufficiently large number of loci with weak effects so that the joint distribution of trait values in parents and offspring is approximately multivariate Gaussian, the covariance in log fitness between parents and offspring is:

Cov PO [ln W g ] ≈ α,β
where C PO g, αβ is the covariance between g α in the parents and g β in their offspring. Using series:

W * ≈ W + α (g α * -g α ) ∂W ∂z α + α≤β C * g, αβ -C g, αβ ∂W ∂C g, αβ (A4)
(note that each (α, β) pair is counted only once in the last sum). From equations A3

and A4, we have

∂W g ∂g α ≈ ∂W ∂z α , ∂ 2 W g ∂g α 2 ≈ 2 ∂W ∂V g, α , ∂ 2 W g ∂g α ∂g β ≈ ∂W ∂C g, αβ (α = β) (A5)
and equation A2 and A5 yield (after dividing both sides by W ):

W g W ≈ 1 + α (g α -g α ) ∂ ln W ∂z α + α≤β [(g α -g α ) (g β -g β ) -C g, αβ ]
∂ ln W ∂C g, αβ .

(A6)

Turelli, M. and N. H. Barton. 1990. Dynamics of polygenic characters under selection.

Theor. Popul. Biol. 38:1-57.

---. 1994. Genetic and statistical analyses of strong selection on polygenic traits: what, me normal? Genetics 138:913-941.

θ

  will thus allow us to explore a continuum between two extreme situations corresponding to two classical models: θ = 0 corresponds to Fisher's geometrical model without mutational bias, with a variance a 2 of mutational effects and possible compensatory effects among different mutations, while θ = 1 corresponds to a situation where all mutations have the same fitness effect (a 2 = 0) and selection thus becomes directional (alleles 1 are disfavored), without any possible compensatory effect among mutations. Note that a 2 and b 2 are simply expressed in terms of s d and θ, as a 2 = 2V s (1 -θ) s d /m and b 2 = 2V s θ s d /m. Furthermore, equation 4 indicates that the parameters a, b, ω 2 and V e should only affect changes in genotype frequencies through the scaled parameters ã = a/ √ 2V s and b = b/ √ 2V s , since genotypic fitnesses become independent of V s when expressed in terms of the scaled phenotypic traits gα = g α / √ 2V s . For a given choice of s d and θ, the results should thus not depend on ω 2 and V e . Using the parameters s d and θ (instead of ã and b) will allow us to change

σ

  and trait α. Equation 31 corresponds to the classical expression describing the effect of selection on correlated characters[START_REF] Lande | Quantitative genetic analysis of multivariate evolution, applied to brain:body size allometry[END_REF]: if selection favors higher values of trait α (∂ ln W /∂z α > 0), a positive genetic covariance between traits α and σ will lead to the evolution of higher values of σ. In our model, directional selection is caused by mutational bias displacing mean phenotypes from the optimum, and thus occurs along the first phenotypic axis of the basis defined by equations 9 and 10 (∂ ln W /∂z α = 0 along all other axes). Because sex increases the response to directional selection by increasing the genetic variance among offspring, trait values tend to be closer to the optimum in individuals that engage more in sex: C g,σα has the same sign as ∂ ln W /∂z α , and ∆ long g σ is thus positive -this term is equivalent to the term in δz in Charlesworth

  A4) where r h,1 is defined in the main text. The term ∆ sel D αβ in the numerator of equation A4 measures the change in D αβ (per generation) due to selection, where D αβ is the contribution of linkage disequilibria to the genetic covariance between traits α and β. As shown in Supplementary File S3, when phenotypes are measured in a basis

sFigure 1 .

 1 Figure 1. Effect of mutational bias, illustrated for n = 2 (the fitness optimum

Figure 2 .

 2 Figure 2. Simulated genetic architecture. Traits affecting fecundity are coded by

Figure 3 .

 3 Figure 3. Top: average mutation load as a function of the mean fitness effect of

Figure 4 .

 4 Figure 4. Average mutation load in asexual (left) and sexual (right) populations as a

Figure 5 .

 5 Figure 5. Average mutation load in asexual (left) and sexual (right) populations as

Figure 6 .

 6 Figure 6. Average mutation load in asexual (left) and sexual (right) populations as

Figure 7 .

 7 Figure 7. Mean rate of sex in the population at equilibrium as a function of the

Figure 8 .

 8 Figure 8. Mean rate of sex at equilibrium as a function of the degree of mutational

Figure 9 .Figure 10 .Figure S1 .

 910S1 Figure 9. Mean rate of sex at equilibrium as a function of the degree of mutational

Figure S5 .

 S5 Figure S5. Same as Figure 5 in the main text, showing W g = exp [-n α=1 g α 2 / (2V s )]

Figure S7 .

 S7 Figure S7. Same as Figure 6 in the main text, showing W g = exp [-n α=1 g α 2 / (2V s )](circles, solid lines) and W Vg = exp [-n α=1 V g, α / (2V s )] (squares, dashed lines).

FILE

  S1: DESCRIPTION OF SIMULATION PROGRAMSThe genome of each individual consists in a single linear chromosome with map length R (average number of cross-overs at meiosis). The loci affecting the n traits under stabilizing selection are biallelic and equally spaced along the chromosome, the genome of an individual at these loci being represented by a set of bits (0 or 1). At the beginning of the simulation, the effects of allele 1 at each locus on the different phenotypes are drawn and stored in a table: as explained above, each locus only affects a subset of m randomly sampled traits, the effect on each of these traits being drawn from a Gaussian distribution with standard deviation a = 2V s (1 -θ) s d /m and average b = 2V s θ s d /m. At the start of each generation, genetic components g α are computed for each individual given its genotype, and environmental components e α are drawn from a Gaussian distribution with mean 0 and variance V e , fixed to 1/n to avoid that fitness reaches very low values when the number of selected traits is large.

  si si αj ≈ C sel si si αj -(∆ sel g αj ) C si si(63) while from equation 29:C sel si si αj ≈ C si si αj + β βγ k,l (C si si αj βk γl -C si si αj C βk γl) 65, and using the fact that C si si αj βk ≈ C si si C αj βk , C si si αj βk γl ≈ C si si C αj βk γl to the first order in V g, s , one obtains that the effect of selection on C si si αj is negligible, which finally leads to C si si αj ≈ C si si αj ≈ 0 at QLE.The effect of selection on C si αj is given by:C si αj ≈ C sel si αj -(∆ sel g si ) (∆ sel g αj ) , βγ k,l (C si αj βk γl -C si αj C βk γl ) ,

  Expressing indirect selection in terms of the effect of sex on the fitness of offspring. The terms between parentheses in equation 58 and 72 (involving ∆ sel D αβ ) provide intuitive understanding of the mechanisms generating indirect selection for sex, but would be difficult to measure in a real population. However, using our hypothesis of weak selection and Gaussian distribution of traits affecting fecundity, these can be expressed in terms of the effect of sex on the fecundity of offspring, that could (at least in principle) be measured in an experimental population. Indeed, a Taylor series on ln W provides the following approximation for the effect of a change in mean breeding values and/or in the genetic variance-covariance matrix on ln W : imagine an experiment where we sample a sufficiently large number of individuals from the population (so that genetic associations within this pool of individuals are representative of associations in the whole population) and let them produce a pool of offspring by sexual reproduction and another pool by asexual reproduction, both pools should have the same mean breeding values (on average), while genetic variances and covariances (measured separately within each pool of offspring) should differ by an amount: ∆ sex/asex C g, αβ = C sex g, αβ -C asex g, αβ =j =k r jk C αj βk (74)

Table 1 :

 1 Parameters and variables of the model.

	N	Population size
	n	Number of selected traits
	m	Degree of pleiotropy of mutations
	ρ = m/n	Scaled pleiotropy
	V e	Environmental variance (on selected traits)
	ω 2

  the following, we use the notation C rep

U

for genetic associations measured after reproduction, but using as "reference values" (the g αj in equation 1) allelic averages after selection (g αj ): for example,

C rep αj βk = E [(g αj -g αj ) (g βk -g βk )]

. Expanding equation 42 and noting that E [g si -g si ] = ∆ rep g si , one obtains:

  so that the variance in ln W g among individuals is:Var [ln W g ] ≈

	α,β	C g, αβ	∂ ln W ∂z α	∂ ln W ∂z β
	+
	α≤β γ≤δ

(C g, αγ C g, βδ + C g, αδ C g, βγ ) ∂ ln W ∂C g, αβ ∂ ln W ∂C g, γδ .

We thank Brian Charlesworth for suggesting considering the effect of mutational bias, two anonymous reviewers for helpful 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 SexChange, ANR-14-CE02-0001). analyzed the model, EV and DR wrote the article.

data is 10.5061/dryad.50q65s0.

constant by decreasing θ as m increases, while m/n is kept constant in the bottom panels by increasing m as n increases (i.e., m = 4, 6 and 10 when n = 40, 60 and 100, respectively). The horizontal dashed lines correspond to equation 24 (1 -e -U ), the dashed curves to equation 22, and the dotted curves to equation 23. Parameter values are as in Figure 3 with θ = 0.1 and n = 50 unless specified otherwise.

FILE S2: APPROXIMATIONS FOR THE EFFECT OF MUTATIONAL BIAS ON MEAN FITNESS

We explain here the derivation of the approximations given in the main text for the effect of mutational bias on the load, assuming a Gaussian fitness function and biallelic loci. Throughout the following, the notation X (also denoted E [X]) stands for the average of the quantity X over all individuals of the population, while the notation Y stands for the expected value of quantity Y in the population at mutation-selection-drift equilibrium. In particular, W is the expected value of the population mean fitness. Assuming that the variance in log-fitness among individuals remains small, we have W ≈ e ln W ; furthermore, assuming that the variance in ln W due to drift is small yields:

(1)

From equations 4 and 18 in the main text, this yields:

Equation 1 shows that the load can be decomposed into the two terms,

representing the decrease in mean fitness due to genetic variance maintained in the population, and to deviations of the mean phenotype from the optimum, respectively. If population size is sufficiently large, the variance of mean phenotypes due to drift should remain small (Lande, 1976;Charlesworth, 2013b), so that g α 2 ≈ g α 2 ; this is confirmed by simulations (results not shown). In the following, we thus derive approximations for L by computing expressions for g α and V g, α in different limit cases.

C PO g, αβ = C g, αβ /2, equation 81 becomes:

Assuming that selection is weak (meaning that the variance in W g is small), we can approximate W g by a Taylor series around g = (g 1 , g 2 , . . .):

where the partial derivatives are taken in g, and the last sum includes α = β. Averaging over all individuals yields W ≈ W g (g) + 1 2 α,β C g, αβ ∂ 2 W g / (∂g α ∂g β ), so that equation A1 can also be written as:

The derivatives of W g in equation can be expressed in terms of derivatives of W (Barton and Turelli, 1991;Turelli and Barton, 1994). Consider the effect of a slight change in the distribution of breeding values g on mean fitness: g α and C g, αβ change to g α * and C * g, αβ , causing mean fitness to change from W to W * . Replacing g α -g α by g α -g α * + g α * -g α in equation A2 and averaging over the new state of the population yields:

Note that terms (g α * -g α ) (g β * -g β ) appearing in the second sum have been neglected, as we assume that g α * -g α is small for all α. Another expression for W * can be obtained by developing W (which is a function of g α = z α and C g, αβ for all α, β) as a Taylor APPENDIX B: SELECTION GRADIENTS WITH ISOTROPIC, GAUSSIAN

FITNESS FUNCTION

From equation 15, we have:

so that:

and

Tr (S + P) -1 ∂ (S + P) ∂G + 1 2 (z -θ) T (S + P) -1 ∂ (S + P) ∂G (S + P) -1 (z -θ)

where Tr stands for the trace of a matrix. If phenotypes are measured in a basis that eliminates covariances among traits, (S + P) -1 is a diagonal matrix, with elements 1/ (V g, α + V s ) on its diagonal (with V s = ω 2 + V e ). In that case, equations B2 and B3 yield equations 16 -18.