# Mutation Rate Evolution in Partially Selfing and Partially Asexual Organisms

^{*}Centre National de la Recherche Scientifique (CNRS), UMI 3614, Evolutionary Biology and Ecology of Algae, Roscoff, 29688 France^{†}Sorbonne Universités, UPMC Université Paris VI, Roscoff, 29680 France

- 1Corresponding author: Centre National de la Recherche Scientifique, Station Biologique de Roscoff, Place Georges Teissier, CS 90074, 29688 Roscoff Cedex, France. E-mail: roze{at}sb-roscoff.fr

## Abstract

Different factors can influence the evolution of the mutation rate of a species: costs associated with DNA replication fidelity, indirect selection caused by the mutations produced (that should generally favor lower mutation rates, given that most mutations affecting fitness are deleterious), and genetic drift, which may render selection acting on weak mutators inefficient. In this paper, we use a two-locus model to compute the strength of indirect selection acting on a modifier locus that affects the mutation rate toward a deleterious allele at a second, linked, locus, in a population undergoing partial selfing or partial clonality. The results show that uniparental reproduction increases the effect of indirect selection for lower mutation rates. Extrapolating to the case of a whole genome with many deleterious alleles, and introducing a direct cost to DNA replication fidelity, the results can be used to compute the evolutionarily stable mutation rate, *U*. In the absence of mutational bias toward higher *U*, the analytical prediction fits well with individual-based, multilocus simulation results. When such a bias is added into the simulations, however, genetic drift may lead to the maintenance of higher mutation rates, and this effect may be amplified in highly selfing or highly clonal populations due to their reduced effective population size.

- clonality
- deleterious mutation
- modifier model
- multilocus population genetics
- self-fertilization

RATES of spontaneous mutation per nucleotide and per cell division span several orders of magnitudes within eukaryotes (*e.g.*, Sung *et al.* 2012; Lynch *et al.* 2016), providing evidence that mutation rates evolve over long timescales. Furthermore, mutation rate variation within a species has been documented in various groups of organisms including bacteria, green algae, and fruit flies (*e.g.*, Demerec 1937; Woodruff *et al.* 1984; Miller 1996; Haag-Liautard *et al.* 2007; Ness *et al.* 2015), suggesting that mutation rates may possibly change rapidly if natural selection can act upon this genetic variation. Changes in mutation rates have indeed been observed during evolution experiments (Sniegowski *et al.* 2000): for example, increased mutation rate in evolving populations of bacteria due to the fixation of mutator genotypes (*e.g.*, Sniegowski *et al.* 1997), or decreased mutation rate in populations of *Drosophila* exposed to X-irradiation during several generations (Nöthel 1987).

Starting with Sturtevant (1937), a number of evolutionary forces that may influence the evolution of mutation rates have been identified (Drake *et al.* 1998; Sniegowski *et al.* 2000; Baer *et al.* 2007; Lynch 2010). Because most mutations affecting fitness are deleterious (Eyre-Walker and Keightley 2007), alleles coding for higher mutation rates should be associated with less fit genetic backgrounds, thus favoring reduced mutation rates. Using a two-locus modifier model in which one locus affects the mutation rate between alleles at a linked locus directly affecting fitness, Kimura (1967) showed that the strength of selection to reduce mutation in a panmictic, diploid population is approximately where is the change in mutation rate caused by the modifier locus, *r* the recombination rate between the two loci, and the heterozygous effect of the deleterious allele, assumed different from zero—see Karlin and McGregor (1974) for the case of a fully recessive deleterious allele. This result was later generalized by Kondrashov (1995), Dawson (1998, 1999), Johnson (1999a), and Lynch (2008) to the case of modifiers changing the deleterious mutation rate over a whole genome. In asexual populations, this effect may be compensated by the higher rate of production of beneficial alleles by mutator genotypes, which may lead to transient increases in mutation rate when mutators hitchhike with the beneficial alleles they created (*e.g.*, Leigh 1970; Eshel 1973; Taddei *et al.* 1997; Tenaillon *et al.* 1999; André and Godelle 2006). In sexual populations, however, recombination destroys the association between mutators and beneficial alleles, and selection for reduced mutation due to the effect of deleterious alleles should generally prevail (Leigh 1970; Johnson 1999b).

The maintenance of nonzero mutation rates is often considered as the result of two opposing forces: selection for reduced mutation rates due to the deleterious effect of most mutations, and the intrinsic cost of DNA replication fidelity (*e.g.*, Kimura 1967; Drake *et al.* 1998; Sniegowski *et al.* 2000; Baer *et al.* 2007). More recently, Lynch (2008, 2011) proposed that the equilibrium value of the mutation rate may instead result from a balance between indirect selection and genetic drift: indeed, once the mutation rate has decreased to a very low level, the strength of selection for further increases in replication fidelity may become weaker than genetic drift. The mutation rate would thus reach higher values in populations with lower effective population size, due to less efficient selection acting on modifier alleles reducing mutation: this agrees with the observation that the mutation rate is lower in species with larger estimated (Lynch 2010; Sung *et al.* 2012; Lynch *et al.* 2016).

Based on the above-mentioned results of Kimura (1967), reproductive systems that reduce effective recombination rates (such as selfing or clonality) should increase the strength of selection for lower mutation rates (as mutators tend to stay longer associated with the deleterious alleles they produce). In the extreme case of full selfing or full clonality, the strength of selection against a mutator allele becomes equivalent to the increase in mutation load that it causes: neglecting drift, this corresponds to the increase in mutation rate caused by the mutator (*e.g.*, Sturtevant 1937; Drake *et al.* 1998). Using multilocus simulations incorporating a cost of replication fidelity, Sloan and Panjeti (2010) showed that the equilibrium deleterious mutation rate is indeed lower in asexual than in sexual populations, generating an indirect benefit for asexuality. With selfing, selection for lower mutation rates should be further enhanced by the increased fitness effect of deleterious alleles due to increased homozygosity. However, background selection may strongly reduce the effective size of highly selfing or clonally reproducing populations (Nordborg 1997; Glémin and Ronfort 2013; Agrawal and Hartfield 2016; Roze 2016), which, according to Lynch’s (2010) hypothesis mentioned above, may possibly increase the equilibrium mutation rate. The overall effect of selfing or clonality on the evolution of mutation rates thus remains unclear, and has been little explored.

In this paper, we extend Kimura’s (1967) two-locus model to compute the strength of indirect selection acting on a mutation modifier locus in a partially selfing or partially clonal diploid population. The results confirm that uniparental reproduction increases selection against mutator alleles due to stronger associations with deleterious alleles. Under partial selfing, the strength of indirect selection generated by closely linked loci can be approximated by replacing *r* and *h* in Kimura (1967)’s result by effective recombination and dominance coefficients and (*e.g.*, Glémin and Ronfort 2013; Hartfield and Glémin 2016; Roze 2016). However, this approximation underestimates the effect of more distant loci, which may become important when the selfing rate is high. We then extrapolate from this two-locus model to derive expressions for the genomic deleterious mutation rate at equilibrium between indirect selection generated by deleterious alleles and the cost of replication fidelity, and show that these expressions correctly predict the outcome of individual-based multilocus simulations. Finally, using two different simulation models with different assumptions on the genetic architecture of the mutation rate, we show that consistent with Lynch’s (2010) hypothesis, populations with lower effective size may maintain higher mutation rates, provided that mutations increasing replication fidelity (antimutator alleles) occur less frequently than those decreasing it (mutator alleles). In some cases, intermediate rates of outcrossing lead to lower mutation rates than obligate outcrossing or obligate selfing/clonality, due to strong reductions in the effective size of highly selfing or clonal populations caused by background selection.

## Methods

### Two-locus model

Our analytical model represents a very large (effectively infinite) population of diploid individuals with discrete generations. As in Kimura (1967), we consider the evolution of a locus (denoted *M*) affecting the mutation rate at a second locus (denoted *A*), which directly affects fitness. Two alleles (denoted 0 and 1) segregate at each locus; we assume that allele 1 at locus *A* is deleterious, reducing fitness by a factor in heterozygotes and in homozygotes. For simplicity, we assume additivity at the mutation modifier locus (locus *M*), the mutation rate at locus *A* being and in individuals with genotype 00, 01, and 11 at locus *M*, respectively. We assume that mutations from 0 to 1 and from 1 to 0 occur at the same rate; however, this hypothesis should not significantly affect the results, as the effect of back mutations will be negligible as long as the deleterious allele stays at low frequency in the population. We also introduce an intrinsic cost of DNA replication fidelity, so that individuals with lower mutation rates pay a fitness cost. For this, we will assume that the fitness of an individual is multiplied by a function, , that increases with the mutation rate, and will consider different forms of cost function. Individuals contribute to the next generation in proportion to their fitness; under partial selfing a proportion *α* of juveniles is produced by selfing, while under partial asexuality a proportion *γ* is produced clonally (the remaining proportion or being produced by outcrossing with random union of gametes). Finally, *r* measures the recombination rate between the two loci. We assume that mutation occurs after selection, before recombination; however, assuming that mutation occurs just after recombination yields the same results (as long as the mutation rate depends on the genotype of the diploid parent).

Following previous works (*e.g.*, Barton and Turelli 1991; Kirkpatrick *et al.* 2002; Roze 2015, 2016), genetic associations within and between loci are defined as follows. We define and as indicator variables that equal 1 if a given individual carries allele 1 at locus *i* on its first or second haplotype, respectively, and 0 otherwise. The frequency of allele 1 at locus *i* in the whole population is thus given by where stands for the average over all individuals. Defining the centered variables and as(1)the genetic association between the sets and of loci present on the two haplotypes of the same individual is given by: ,(2)where(3)(note that ), and where sets and may be the empty set ∅, *M*, *A* or Associations between genes present on the same haplotype of an individual () will be simply denoted For example, is a measure of the departure from Hardy-Weinberg equilibrium at locus *M*, while represents the linkage disequilibrium (LD) between loci *M* and *A* (genetic association between alleles present on the same haplotype). Similarly, measures the association between alleles at loci *M* and *A* present on different haplotypes of the same individual.

In the following, we assume that both loci have weak effects (*s*, small) and derive an expression for the change in (the frequency of allele 1 at locus *M*) to the first order in *s* and We will see that this expression includes different forms of genetic associations. Assuming that the effective recombination rate is large relative to we will then use a quasi-linkage equilibrium approximation (QLE) to express these associations in terms of allele frequencies, and of the different parameters of the model. Finally, the result will be extrapolated to compute the overall strength of selection on a modifier allele affecting the mutation rate at a large number of selected loci, assuming that genetic associations between those loci can be neglected.

### Multilocus simulations

Our simulation program (written in C++ and available from Dryad) is modified from Roze (2015, 2016), and represents a finite population of *N* diploid individuals whose genome consists in a linear chromosome along which deleterious mutations occur every generation. For simplicity, all mutations have the same selection and dominance coefficients (*s*, *h*). A mutation modifier locus is located at the midpoint of the chromosome, and controls the deleterious mutation rate (the mutation rate of an individual being the average of the values coded by its two modifier alleles). At the start of each generation, the fitness of every individual is computed as(4)where *U* is the deleterious mutation rate of the individual (per haploid genome); the function representing the cost of replication fidelity; and *i*, *j* are the number of heterozygous and homozygous deleterious alleles present in the genome of the individual. In general, we will use the cost function but different functions will also be considered (as explained in the *Results* section). To form each of the *N* juveniles of the next generation, an individual is sampled randomly to serve as a maternal parent. If the fitness of the individual (divided by the maximal fitness in the population) is higher than a random number sampled from a uniform distribution between 0 and 1, the individual is retained, otherwise another individual is sampled until the test is satisfied. Under partial selfing, the mother self-fertilizes with probability *α*, in which case the new individual is formed by two recombinant chromosomes from the same parent. Under partial asexuality, the mother reproduces clonally with probability *γ*, in which case the genome of the new individual is a copy of the maternal genome. If the mother reproduces by outcrossing (with probability or ), a second individual is sampled using the same procedure as above to serve as a father, and the genome of the new individual is generated from recombinant chromosomes from both parents. During meiosis, the number of cross-overs is sampled from a Poisson distribution with parameter *R* (genome map length, in Morgans), and the position of each cross-over is sampled from a uniform distribution. The parameter *R* will typically take large values () in order to mimic a whole genome with multiple chromosomes. Deleterious mutations occur once the parents have been selected, before recombination (note that different offspring from the same parent will carry different new mutations). The number of new deleterious mutations on each chromosome is sampled from a Poisson distribution whose parameter corresponds to the mutation rate of the parent, and the position of each new mutation is sampled from a uniform distribution along the chromosome. Back mutations do not occur, and any deleterious allele that has reached fixation is removed from the population in order to increase execution speed.

During a number of preliminary generations (usually 2000), the deleterious mutation rate of each individual (per haploid genome) is set to Then, new alleles coding for different mutation rates can appear at the modifier locus (at rate per generation, where *U* is the deleterious mutation rate of the individual). When a mutation occurs at the modifier locus, the mutation rate coded by the new allele is sampled from a Gaussian distribution centered on the value of the allele before mutation, with variance (if the new value is negative, it is set to zero). As explained in the *Results* section, different mutational models were also considered, including a bias toward higher values of *U* and scaling of with *U*. The program generally runs for generations, the equilibrium mutation rate being computed by averaging over the last generations.

A second simulation program considers a different genetic architecture for the mutation rate: instead of being coded by a single locus, *U* depends on 1000 biallelic loci evenly spaced along the chromosome. Alleles at each of these loci are denoted 0 and 1; genotypes carrying alleles 0 at all loci have mutation rate At each locus, allele 1 increases the mutation rate, by an amount that is sampled (independently for each locus) from an exponential distribution with parameter *λ*. The effects of alleles 1 at the same or at different loci are additive. During the first 2000 generations, all loci affecting the mutation rate are fixed for allele 0; then, during generations mutations occur at rate at each of these loci (mutations and back mutations occur at the same rate). Selection and recombination are implemented as in the previous program.

### Data availability

The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article. Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.1j6b0.

## Results

### Change in frequency at the mutation modifier locus

In the Appendix, we show that an expression for the change in frequency of allele 1 at locus *M*, to the first order in and *s* is given by:(5)The term on the first line of Equation 5 represents the effect of the cost of replication fidelity, favoring alleles that increase the mutation rate. This direct selective pressure increases with the rate of inbreeding ( factor), due to increased homozygosity at locus *M*. The term on the second line represents the effect of indirect selection disfavoring mutator alleles, as these alleles tend to be more often associated with the deleterious allele at the second locus. Indeed, and represent the association between allele 1 at locus *M* and the deleterious allele at locus *A* on the same or on the other haplotype, whereas represents the association between allele 1 at locus *M* and homozygotes for the deleterious allele at locus *A*. We now derive expressions for these associations at QLE.

### Expressions for genetic associations

In the following, and denote genetic associations measured after selection, mutation, and recombination/segregation (respectively). Recursions for genetic associations over each step of the life cycle are computed to the first order in *s* and We will also assume that the deleterious allele stays at low frequency and neglect terms in Under high effective recombination, it is sufficient to express associations to the first order in neglecting terms in *s*. However, these expressions diverge when the effective recombination rate tends to zero [due to terms or in their denominators]. To obtain more accurate expressions for the case where the effective recombination rate is of order *s* (assuming that the deleterious allele stays at mutation-selection balance and that is sufficiently small, so that the QLE approximation holds), we include terms in *s* in the recursions for genetic associations, by computing the effect of selection on these associations.

#### Selection:

To leading order, the effect of selection on can be written as: .(6)Note that Equation 6 is an approximation, as it neglects the change due to selection of allele frequencies and that appear in and (see Equation 1). However, taking this change in allele frequencies into account would introduce a term of order (*e.g.*, Equation 23 in Barton and Turelli 1991) that can be neglected here. From Equations 6 and A6, neglecting terms in and terms proportional to one obtains:(7)Similarly, one arrives at:,(8)(9)Equations 7–9 show that selection tends to decrease genetic associations between the two loci, as it reduces the frequency of the deleterious allele in the population.

#### Mutation:

The effect of genotype-dependent mutation on genetic associations can be computed as follows. Mutation changes the frequency of allele 1 at locus *A* to (where is the average mutation rate at locus *A*), while in a given individual, changes to with probability *u* (the mutation rate of the individual at locus *A*), and remains unchanged with probability Therefore, after mutation is given by:(10)where is the average over all individuals before mutation. Replacing by and *u* by Equation 10 yields [after neglecting terms of order and terms in ]:(11)Similarly, one obtains:(12)(13)(indeed, one can show that the term in of is of order ). Equations 11 and 12 show that the modifier effect generates an association between the allele increasing mutation and the deleterious allele at the other locus, either on the same or on the other haplotype. The buildup of and is stronger under inbreeding (term in Equations 11 and 12), due to increased homozygosity at the modifier locus: the probability that a deleterious allele is present in the same genome as an allele coding for higher mutation is increased when the high-mutation allele is also present on the other haplotype.

#### Recombination/segregation:

The effects of recombination and segregation depend on the reproductive system. Under partial selfing, we have:(14)(15)(16)Equation 16 assuming that is small. Under partial asexuality:

(17)(18)#### QLE expressions:

The expressions given above can be used to compute solutions for and at QLE, corresponding to the equilibrium values of genetic associations for given values of allele frequencies (under the assumption that associations equilibrate fast relative to changes in allele frequencies). For the case of a partially selfing population, one obtains (assuming and thus neglecting terms in ):(19)(20)(21)with and Equations 19–21 can be used to obtain the strength of indirect selection acting on the modifier locus, given by (from Equation 5):(22)In the absence of selfing (), simplifies to:(23)which agrees with the result obtained by Kimura when the modifier effect is weak (Equation 2 in Kimura 1967). Under complete selfing (), simplifies to

When linkage is tight (small *r*), a separation of timescales argument yields (Nordborg 1997; Roze 2016): this may also be seen from Equations 19 and 20, neglecting terms in *r* in the numerators. From this, one obtains that is equivalent to the result obtained under random mating (Equation 23), replacing by *h* by and *r* by This can also be obtained directly by neglecting the term in *r* in the numerator of Equation 22, and the term in in the denominator. Figure 1A compares the prediction from Equation 22 and the prediction obtained by replacing *h*, and *r* in Equation 23 by effective parameters: both yield undistinguishable results when linkage is sufficiently tight, but discrepancies appear when loci are loosely linked.

Under partial asexuality, genetic associations at QLE are given by:(24)(25)As shown by Figure 1B, the strength of indirect selection is roughly equivalent to the expression obtained under random mating, replacing *r* by as long as the rate of clonal reproduction *γ* stays moderate. This is not the case for higher values of *γ*, however, due to the extra contribution of the association (mutations generated on the other haplotype, that remain associated to the modifier allele due to clonal reproduction). As under complete selfing, one obtains that under full asexuality (). In both cases, the strength of indirect selection in the absence of recombination corresponds to the increase in mutation load caused by the modifier (as the modifier allele stays associated with the deleterious alleles it generates). Under selfing, the increase in load equals the increase in haploid mutation rate (since under full selfing), which is as individuals quickly become homozygous at the modifier locus. Under clonality, the increase in load is twice the increase in haploid mutation rate (since assuming *h* is significantly >0), which is as individuals stay heterozygous at the modifier locus.

### Multilocus extrapolation and simulation results

These two-locus results can be extrapolated to the case of a modifier affecting the mutation rate in the whole genome by integrating over the genetic map. If the map length is sufficiently large, the overall effect can be approximated by assuming free recombination among all loci, replacing *r* by and by (the increase in haploid mutation rate caused by allele 1 at the modifier locus) in the expressions above. The evolutionarily stable mutation rate (at which indirect selection to reduce the mutation rate exactly balances the cost of replication fidelity) can then be obtained by solving for *U*, where the strength of direct selection is given by assuming is small (see Equation 5). Under partial selfing, and assuming free recombination, this yields(26)when the cost function is given by so that Equation 26 simplifies to when and to when The equivalent expression for partial asexuality is given by:(27)simplifying to under full asexuality ().

Figure 2 shows that the predicted value for the evolutionarily stable mutation rate obtained by integrating Equations 22, 24, and 25 over the genetic map (see *Mathematica* notebook in Supplemental Material, File S2 for the integration) generally fits well with the multilocus simulation results (using our first simulation program, with a single modifier locus). With a genome map length of 20M, the simpler expressions obtained assuming free recombination (Equations 26 and 27) also provide accurate predictions: discrepancies appear for lower values of *s*, due to the fact that deleterious alleles segregate at higher numbers of loci, increasing the number of deleterious alleles that are closely linked to the mutation modifier locus. Figure 2, C and D show that integrating Equations 22, 24, and 25 over the genetic map also provides accurate predictions for lower values of map length *R*. As shown by Figure 2E, discrepancies between analytical and simulation results appear for low values of *h* and intermediate selfing rates ( 0.4 in Figure 2E): these discrepancies are possibly generated by identity disequilibria between selected loci (correlations in homozygosity), which are neglected in the analytical model (the discrepancies observed for in Figure 2A may also be caused by identity disequilibria).

Figure S1 in File S1 shows the approximations obtained for when replacing *h* and *r* by the effective parameters and (under partial selfing) and *r* by (under partial clonality) in the expression for indirect selection under random mating (Equation 23). Although these approximations tend to overestimate by underestimating the strength of indirect selection generated by distant loci, they often stay relatively close to the more exact expressions given above, the discrepancy being stronger for intermediate selfing or clonality rates, and for weaker strength of selection against deleterious alleles.

### Effects of population size and mutational bias

Figure 3, A and B show that changing population size *N* from to or to has little effect on the average mutation rate at equilibrium (although the variance of *U* around its average value increases as *N* decreases). This may seem at odds with the prediction of Lynch (2010) mentioned in the Introduction, which states that the mutation rate should be lower in populations with larger in regimes where a substantial proportion of mutations changing *U* are significantly affected by genetic drift. Indeed, averaging over the distribution of mutational effects at the modifier locus, the mean value of (corresponding to the average strength of indirect selection acting on a new modifier allele under random mating and free recombination) is close to in the simulations, and thus of the same order of magnitude as the strength of genetic drift (at least for and ). Nevertheless, decreasing *N* from to does not significantly affect the average value of *U* at equilibrium (see Figure 3). Similar results are obtained when using different forms of cost function In Figure 3, C and D, where *c* is set to so that under random mating for the parameter values used in Figure 3, according to our approximations: as illustrated by Figure 4, the selection gradient obtained () is less steep around than with the exponential cost function used in Figure 2 and Figure 3, A and B. Finally, in Figure 3, E and F, yielding a linear selection gradient (). Parameters *a* and *b* were set to and so that under random mating for the parameter values used in Figure 3, while the slope of the selection gradient at is the same as with the cost function used in Figure 2 and Figure 3, A and B (see Figure 4). Because the mutation rate evolves to very low levels for these parameter values at sufficiently high selfing or clonality rates, we maintained a minimum mutation rate of at the modifier locus in the simulation program to prevent that the population remains stuck in the absorbing state of perfect replication fidelity (). Figure S2 in File S1 shows the same simulation results as Figure 3, E and F, displayed on a log scale.

The reason for the limited effect of drift on the average value of *U* observed in Figure 3 is the absence of mutational bias at the modifier locus: indeed, mutations increase or decrease *U* with the same probability. Drift may have stronger effects when mutations affecting *U* tend to occur more often in a particular direction (Zhang and Hill 2008; Lynch 2011; Charlesworth 2013): most likely in the direction of increased values of *U*, as it should be easier to impair DNA replication fidelity than to improve it. Indeed, when such a mutational bias is added in the simulation program (by introducing a parameter *β* such that a fraction of mutations at the modifier locus tend to increase *U*), the mutation rate evolves toward higher values when population size is sufficiently small. When this is the case, *U* keeps increasing unless one assumes that the average size of mutational steps is proportional to *U*, so that fewer mutations may fix by drift as *U* reaches higher values (*e.g.*, Lynch 2011). Figure 5 shows the results of simulations in which is sampled in a half-Gaussian distribution with SD (for both and ), where *U* is the mutation rate coded by the modifier allele before mutation. For the parameter values used in Figure 5, the equilibrium mutation rate is slightly higher for than for but generally remains close to the evolutionarily stable strategy (ESS) value (see Figure S3 in File S1 for the same results shown on a log scale). Genetic drift has a much stronger effect on the equilibrium value of *U* when however. As expected, increasing the degree of mutational bias (by increasing *β*) or decreasing the average size of mutational steps (by reducing ) amplifies the effects of drift, causing higher values of *U* to evolve. Intermediate selfing or clonality rates bring *U* closer to its ESS value (by increasing the strength of indirect selection), but *U* may increase again as *α* or *γ* approach 1, due to background selection amplifying the effect of drift. In some simulations, background selection caused a runaway process in which the reduction in leads to elevated mutation rate, further reducing When this is the case, *U* reaches very high values, and the program has to be stopped manually: this happened for clonality rates higher than the right-most points in Figure 6, B, D, and F, and for and in Figure 6E (while an equilibrium was reached for ). As shown by Figure S4, Figure S5, Figure S6, and Figure S7 in File S1, qualitatively similar results were obtained under the different cost functions shown in Figure 4.

Similar results were also obtained using our second simulation program, representing a more realistic genetic architecture in which *U* is controlled by *L* biallelic loci. Since the minimal mutation rate (corresponding to the mutation rate of a genotype carrying allele 0 at each mutation modifier locus) is and since the heterozygous effect of each modifier locus is sampled from an exponential distribution with parameter *λ* (whose average is ), the average value of *U* should thus be close to in regimes where the evolution of *U* is mainly controlled by drift (assuming large *L*, and additivity within and between modifier loci). Below mutations at modifier loci thus tend to increase *U* (mutational bias). Figure 6 shows simulation results for and (so that ). As can be seen on the figure, *U* becomes closer to its mutation-drift equilibrium value as *N* decreases. Again, increasing the selfing rate or the clonality rate tends to reduce *U* by increasing the strength of indirect selection; however, above a given threshold for *α* or *γ* (that depends on population size), *U* increases as the selfing or clonality rate increases (due to background selection effects). Note that we could not obtain estimates for the equilibrium mutation rate under high clonality rates (Figure 6B, on the right of the right-most points) because deleterious alleles tend to accumulate in the heterozygous state in the population, causing the program to become increasingly slow.

## Discussion

It is widely accepted that mutation rates are maintained at low levels to avoid producing an overly strong burden of deleterious alleles. In this paper, we confirm that this deterministic force favoring lower mutation rates is increased by uniparental reproduction, and compute the strength of this effect in populations with intermediate selfing rates or clonality rates. In agreement with previous separation of timescales arguments (Nordborg 1997; Roze 2016), when linkage between loci is sufficiently tight, the result obtained under partial selfing becomes equivalent to the expression obtained under random mating (Kimura 1967), replacing the dominance coefficient of deleterious alleles and recombination rates by the effective parameters and however, this expression underestimates the strength of indirect selection generated by loosely linked loci. Introducing a direct fitness cost associated with DNA replication fidelity, we could obtain simple approximations for the evolutionarily stable mutation rate, which were confirmed by multilocus, individual-based simulations, in the absence of mutational bias at the mutation modifier locus.

When a mutational bias toward lower fidelity of DNA replication (*i.e.*, higher mutation rate) is added into the model, the average value of the mutation rate at equilibrium becomes more sensitive to genetic drift, in agreement with general results on the evolution of quantitative traits under mutation, selection and drift (Zhang and Hill 2008; Charlesworth 2013). In that case, the mutation rate *U* may reach high values when mutations affecting *U* have a weak effect, so that the effect of indirect selection acting on these mutations becomes weaker than the strength of genetic drift: this is the essence of the argument proposed by Lynch to explain the observed negative correlation between estimated effective population size and the mutation rate (*e.g.*, Lynch 2010; Sung *et al.* 2012; Lynch *et al.* 2016). However, the mutational bias by itself is not a sufficient condition for *U* to stabilize around a value that depends on This will occur under the extra assumption that the average size of mutational steps at mutation modifier loci increases with the mutation rate, so that the relative effect of drift at these loci becomes weaker as *U* increases (Lynch 2011, Figure 5 of the present article). Alternatively, when *U* is affected by a sufficiently large number of loci with a distribution of effects, the fraction of loci at which indirect selection is weaker than drift will decrease as increases, which may also generate a negative relation between the average value of *U* and (Figure 6). After adding these ingredients into our simulation programs, we observed three possible types of outcomes: either drift has only a limited effect, and *U* stays close to its deterministic equilibrium (at which indirect selection generated by deleterious alleles exactly balances the cost of replication fidelity), or *U* stabilizes around a higher value that depends on or a runaway process occurs, under which drift causes the evolution of higher *U*, in turn reducing through background selection effects, causing further increase in *U* (which should eventually lead to population extinction through mutational meltdown).

Given that per-nucleotide mutation rates are very low in most species, the hypothesis of a mutational bias toward higher mutation rates seems reasonable. Very little is known on the distribution of the effects of mutations affecting *U*, however, or on how this distribution may change as *U* evolves. Interactions between mutations affecting DNA repair pathways have been demonstrated in bacteria and yeast. In some cases, positive epistasis (on the mutation rate) has been shown between mutations in genes with partially redundant effects, such as *MutM* and *MutY* in *Escherichia coli* (Michaels *et al.* 1992; Fowler *et al.* 2003), or *MSH3* and *MSH6* in *Saccharomyces cerevisiae* (Marsischky *et al.* 1996). These examples provide possible scenarios under which the effect of a modifier allele would increase with the baseline mutation rate. However, examples of negative epistasis also exist, for example when a repair pathway involves the combined activity of two proteins (such as *MSH2* and *MSH6* in yeast, Marsischky *et al.* 1996). Overall, we still lack a clear picture of how the average effect of mutator/antimutator alleles should change with the baseline mutation rate. Furthermore, these studies generally focus on mutators with large effects, which may not be representative of the majority of mutations affecting *U*. Obtaining more detailed information on the genetic architecture of mutation rate variation within natural populations would represent an important progress, but remains a formidable task.

While our results show that the effects of partial selfing and partial clonality are very similar in regimes where drift has only a limited effect, differences appear when drift is stronger and may lead to mutation accumulation. Muller’s ratchet occurred in some of our simulations with partial or complete selfing: for and in Figure 2, and in different cases with in Figure 5 (for in Figure 5A, and 1 in Figure 5B, 0.3 and in Figure 5E). It also occurred for 2000, 3000, and high *α* () in Figure 6. Because our simulation program removes fixed mutations, it could still continue to run and *U* generally stabilized, except for 0.3, and in Figure 5E, where *U* increased to very high values and the program had to be stopped when However, in all these cases, any real population would eventually reach extinction due to mutation accumulation. By contrast, at high clonality rates mutations tend to accumulate in the heterozygous state (even when ), an effect already observed in previous studies by Charlesworth *et al.* (1993a,b) and Roze and Michod (2010). Because these mutations could in principle still be removed from the population by rare segregation events, the simulation program does not eliminate them and becomes increasingly slow. This occurred in nearly all cases with in Figure 2 and Figure 3, although the mutation rate reached an equilibrium before the program had to be stopped. It also occurred for values of *γ* higher than the right-most points in Figure 5 and Figure 6 (at a faster rate as *γ* increased), in which case the program had to be stopped before *U* had reached equilibrium. Again, in all these situations, the population would eventually go extinct by mutational meltdown. These results outline two important differences between partial selfing and partial clonality: (1) the mutation accumulation regime is reached sooner under partial clonality than under partial selfing as the rate of uniparental reproduction increases (due to the absence of segregation), and (2) a runaway process leading to very high mutation rates may occur at moderate selfing rates ( 0.3 in Figure 5E, see also Figure S2 and Figure S3 in File S1), while it does not occur at higher selfing rates. This last effect is not observed under partial clonality, and could possibly be due to identity disequilibria between selected loci reducing the efficiency of selection (*e.g.*, Roze 2015).

Provided that mutation rate polymorphism exists within populations, the evolution of *U* could in principle interact with the evolution of different aspects of reproductive systems, such as reproductive modes or mating systems. In regimes where the effect of drift on the evolution of *U* stays negligible, this should favor the evolution of uniparental reproduction, as it should be associated with lower mutation rates (Sloan and Panjeti 2010). However, the evolution of *U* may have opposite effects in regimes where it is more strongly affected by drift, and where high rates of uniparental reproduction may trigger the evolution of higher mutation rates, through stronger background selection effects. Whether the evolution of the mutation rate would occur on a sufficiently fast timescale to significantly affect the evolution of reproductive systems should in principle depend on the genetic architecture of *U* and of the reproductive system; however, this should be explored more rigorously using theoretical approaches. Even if mutation rate evolution is not fast enough to have a significant impact on evolutionary transitions between mating systems, *U* may change in response to a switch in reproductive system. This may affect the species-level selection component acting on the evolution of reproductive systems, for example, by accelerating the extinction of selfing or asexual lineages by mutational meltdown if *U* reaches higher levels due to stronger drift. More generally, it may affect the long-term evolutionary potential of selfing or asexual species, or the relation between the selfing rate and level of inbreeding depression across species. It would thus be of particular interest to obtain mutation rate estimates from different pairs of closely related species with contrasted reproductive systems, in order to see if a general pattern emerges.

Finally, our model makes a number of assumptions on selection against deleterious alleles: in particular, all deleterious alleles have the same selection and dominance coefficient, while drift has no significant effect on their equilibrium frequency. Given the concave shape of the relation between the selection coefficient of deleterious alleles *s* and the strength of indirect selection acting on the mutation rate, introducing variability in *s* across loci (while keeping the average constant) should in principle reduce the overall strength of selection for lower mutation rates. However, our deterministic model cannot be used to predict the effect of deleterious alleles for which (whose frequency is significantly affected by drift), while our infinite sites simulation program cannot deal with very low *s* values as individuals then carry very large numbers of mutations, causing the program to become extremely slow. Given that an important proportion of mutations may possibly fall in the parameter region (*e.g.*, Eyre-Walker and Keightley 2007), it would be important to explore the effect of such weakly selected deleterious alleles on the evolution of mutation rate modifiers. This may introduce new effects of reproductive systems on mutation rate evolution, as the reproductive system may affect the proportion of mutations on which natural selection is effective (by affecting ). Our model also assumes that all mutations affecting fitness are deleterious: while previous theoretical work has shown that beneficial mutations should only have a minor role on the evolution of mutation rates in sexual, outcrossing, populations (Leigh 1970; Johnson 1999b), their effect should become more important in populations undergoing high rates of selfing or clonality (since mutator alleles can stay associated with the beneficial alleles they produced), and may increase the equilibrium mutation rate in such populations. Last, as in most theoretical studies of mutation rate evolution, we have neglected epistatic interactions between selected mutations: in particular, our model does not take into account possible compensatory effects between deleterious alleles (*e.g.*, reciprocal sign epistasis). Allowing the sign of the fitness effect of mutations to depend on the genetic background (which typically occurs in models of directional selection acting on quantitative traits) may affect the selective forces acting on mutation modifier loci. Exploring the evolution of the mutation rate under more realistic assumptions on the genetic architecture of fitness will be the subject of future work.

## Acknowledgments

We thank Jean-Nicolas Jasmin and two anonymous reviewers for helpful discussions and comments, and the bioinformatics and computing service of Roscoff’s Biological Station (Abims platform) for computing time. This work was supported by the French Agence Nationale de la Recherche (project TRANS, ANR-11-BSV7-013 and project Clonix, ANR-11-BSV7-007), and by a post-doctoral fellowship from the Conseil Général du Finistère to C.G.

## Appendix

### Change in Frequency at the Modifier Locus

Using the definitions given in the main text, the fitness of an individual can be written as:(A1)where represents the cost of replication fidelity and *u* is the mutation rate of the individual, given by:(A2)Replacing by and noting that is the average mutation rate at locus *A*, we have and a Taylor series of to the first order in yields:(A3)Since when Equation A4 may also be written as:(A4)Replacing by in Equation A1, one then obtains (to the first order in ):(A5)with Denoting the average fitness in the population, this yields (to the first order in and *s*):(A6)The change in frequency of allele 1 at locus *M* (over one generation) is given by:(A7)where is the average over all individuals just before selection. From Equations 2 and A6, and using the fact that repeated indices appearing in genetic associations can be eliminated using the relation [with *e.g.*, Equation 5 in Kirkpatrick *et al.* (2002)], one arrives at:(A8)The first term of Equation A8 represents the effect of direct selection acting at locus *M* (due to the cost of replication fidelity), while the other terms (involving genetic associations) correspond to indirect selection. Associations and are of order and the term on the second line of Equation A8 is thus of order Furthermore, one can show that the association is of order *s* (*e.g.*, Roze 2015), and the third line of Equation A8 can thus be neglected [terms of order and ]. Finally, the term that appears on the last line of Equation A8 can also be written where is the identity disequilibrium between loci *M* and *A*, generated by partial selfing (Roze 2015). Similarly, one can show that the association (measuring the excess homozygosity at locus *M*) that appears on the first line of Equation A8 is affected by the identity disequilibrium through a term proportional to —see Equation 5 in Roze (2015). However, we show in the main text that contrarily to the expressions for and at QLE do not tend to zero when the frequency of the deleterious allele tends to zero. Therefore, assuming that the deleterious allele stays at low frequency ( small), we may neglect terms proportional to and thus neglect terms involving the identity disequilibrium. In this case, can be written as where *F* is the inbreeding coefficient (*e.g.*, Roze 2015), and Equation A8 simplifies to:

## Footnotes

Supplemental material is available online at www.genetics.org/lookup/suppl/doi:10.1534/genetics.117.300346/-/DC1.

*Communicating editor: J. Hermisson*

- Received April 19, 2017.
- Accepted September 28, 2017.

- Copyright © 2017 by the Genetics Society of America