The Degeneration of Asexual Haploid Populations and the Speed of Muller's Ratchet
 Isabel Gordo⇓ and
 Brian Charlesworth
 Institute of Cell, Animal and Population Biology, University of Edinburgh, Edinburgh EH9 3JT, United Kingdom
 Corresponding author: Isabel Gordo, Institute of Cell, Animal and Population Biology, Ashworth Laboratories, King's Bldgs., W. Mains Rd., Edinburgh EH9 3JT, United Kingdom. Email: i.gordo{at}ed.ac.uk
Abstract
The accumulation of deleterious mutations due to the process known as Muller's ratchet can lead to the degeneration of nonrecombining populations. We present an analytical approximation for the rate at which this process is expected to occur in a haploid population. The approximation is based on a diffusion equation and is valid when N exp(−u/s) ⪢ 1, where N is the population size, u is the rate at which deleterious mutations occur, and s is the effect of each mutation on fitness. Simulation results are presented to show that the approximation estimates the rate of the process better than previous approximations for values of mutation rates and selection coefficients that are compatible with the biological data. Under certain conditions, the ratchet can turn at a biologically significant rate when the deterministic equilibrium number of individuals free of mutations is substantially >100. The relevance of this process for the degeneration of Y or neoY chromosomes is discussed.
NEW mutations arise continuously within populations, the vast majority probably being slightly deleterious (Crow 1993). An effectively infinite asexual population subject to recurrent deleterious mutations will achieve an equilibrium resulting from the continuous appearance of new mutations opposed by selection against them: deterministic mutationselection balance. With independent and identical effects of each mutation, the equilibrium number of mutations in a haploid randomly mating population follows a Poisson distribution with mean u/s (Kimura and Maruyama 1966; Haigh 1978), where u is the per genome mutation rate and s is the selection coefficient against a deleterious mutation. At this deterministic equilibrium, the number of individuals free of mutations (the best or “leastloaded” class) in a large population of N breeding adults, is n_{0} = N exp(−u/s).
But random genetic drift plays a role in a finite population, and it may perturb this equilibrium, leading to the loss of the best class. In the absence of recombination, and with the reasonable assumption that back mutation is negligible for strongly selected mutations, the loss of this class is irreversible and mutations will continually accumulate in the population, leading to the decline of its mean fitness. This is the process known as Muller's ratchet (Muller 1964; Felsenstein 1974). Once the best class is lost, the new leastloaded class is now the one that has one mutation, but this is also subject to stochastic loss, so that a repetition of successive losses of the leastloaded class can be seen as successive clicks of the ratchet.
One of the important questions in this process concerns the rate or speed at which it operates, or: how much time does it take for the ratchet to click one notch? Because the degeneration of the Y chromosome (Charlesworth 1978, 1996; Rice 1994; Charlesworth and Charlesworth 1997, 1998) and the fate of asexual populations (Pamiloet al. 1987; Gabrielet al. 1993; Lynch et al. 1993, 1995) may involve the operation of Muller's ratchet, the quantification of its rate is of great biological importance.
Haigh (1978) suggested that the most important parameter for the ratchet mechanism would be n_{0}, because it is the loss of the best class that drives the process; the smaller the n_{0}, the faster the ratchet is likely to be. Although Haigh suggested an expression [Equation 9a in Haigh (1978)] for the time between clicks of the ratchet, this is basically a fit to his simulation results. Bell (1988) also suggested an expression based on a fit of the time to extinction of the best class as a function of n_{0} (roughly 10 n_{0}). Several attempts toward the quantification of the process have subsequently been made, either using a quantitative genetics approach for estimating the rate of change of the average number of mutations (Pamiloet al. 1987; Gabrielet al. 1993; Lynchet al. 1993; Higgs and Woodcock 1995; PrügelBennett 1997) or using diffusion theory to calculate the mean time to loss of the leastloaded class (Stephanet al. 1993; Charlesworth and Charlesworth 1997).
In this article, it is shown that the size of the best class is not sufficient to predict the speed of the ratchet. For the same value of n_{0}, the ratchet can turn at very different speeds. An approximation based on a diffusion equation, for the case n_{0} > 1, which almost certainly applies to the evolution of the Y chromosome (Charlesworth 1996), is presented, together with simulations of asexual haploid populations to check its validity. Gessler (1995) has derived an approximation that seems to work well for the case n_{0} < 1. Comparisons between the simulation results and the predictions from the analytical expression suggest that the formulas seem to predict the rate of the process better than the previous formulas for a region of parameter space that is of biological interest (Charlesworth and Charlesworth 1997).
APPROXIMATION BASED ON THE DIFFUSION EQUATION
We start with a haploid asexual population at equilibrium under mutationselection balance, with x_{0} = exp(−u/s) being the frequency of individuals in the leastloaded class. The existence of this equilibrium requires that the number of individuals in the best class (n_{0} = Nx_{0}) be >1. When n_{0} < 1, this equilibrium may not be approached in a finite population, because it requires the existence of individuals that have a very low probability of actually being present (Gessler 1995; Gessler and Xu 1999).
The way the frequency of the best class varies through time is dictated by mutation taking it below the equilibrium value, selection restoring it to that value, and by stochastic fluctuations due to drift. We wish to quantify how much time it takes for the frequency of the best class (from now on denoted by x), with initial value x_{0}, to reach the value zero. One way to do this is to use a diffusion equation for the density function of the time until absorption occurs, subject to the condition that x = 0 is the only absorbing state (Ewens 1979, Equations 4.39, 4.40, p. 123). To solve this equation, one has to evaluate the deterministic change (drift coefficient) and stochastic variance (diffusion coefficient) in x.
Assuming a WrightFisher population, the diffusion coefficient is just the variance due to binomial sampling of N individuals from the previous generation (Stephanet al. 1993; Charlesworth and Charlesworth 1997),
The drift coefficient, representing the expected change in x due to mutation and selection (Stephanet al. 1993; Charlesworth and Charlesworth 1997), is
Let us now make the simplifying assumption that the changes in mean fitness are sufficiently small that they can be approximated by small perturbations from the equilibrium value (Stephanet al. 1993). This implies that the system is close to its equilibrium state most of the time, as is supported by our simulations (see Figure 9). We assume that the perturbations in w are mostly due to small fluctuations in the leastloaded class. This assumption was employed by Stephan et al. (1993) and Charlesworth and Charlesworth (1997), and is justified in practice by the observation that, for large N, the distribution among the classes that are present at any time remains close to the Poisson distribution given by the deterministic equilibrium formula [see Table 1 and Figure 4 of Charlesworth and Charlesworth (1997)]. We may express the mean fitness close to equilibrium as a Taylor expansion in x/x_{0}:
Taking just the linear term in (2a), we obtain an approximation for the reduction in mean fitness below its equilibrium value as
When the frequency x is above its equilibrium value, the system responds with a reduction of mean fitness
In fact, the approach to the neighborhood of equilibrium takes some time (Haigh 1978), and just after a click the new best class is above its equilibrium value, so that
We now may write the drift coefficient as
Using these drift and diffusion coefficients, the time spent in the frequency interval [0, x_{0}] (Ewens 1979, Equation 4.39) is
Using expressions 3a and 3b, and evaluating the integrals numerically for a given population size, mutation rate, and selection coefficient, we obtain the expected time to loss of the leastloaded class as T(N, u, s) = T_{0,x0} + T_{x0,1}.
SIMULATION METHODS
For a given population size (N), genomic mutation rate to deleterious mutations (u), and selection coefficient against each mutation (s), haploid asexual populations were simulated starting at mutationselection equilibrium (Kimura and Maruyama 1966); i.e., the number of individuals in the class with m mutations is
Assuming that the sequence of events is mutation, reproduction, and selection, populations were then run for 100 generations. After this initial period, populations were run for >2000 generations and up to 100,000 generations for conditions under which the ratchet clicks slowly, so that the average time between clicks of the ratchet could be measured. Every generation, the number of mutations in every individual is counted and the number of individuals with the least number of mutations (leastloaded class) is recorded. If, at a given generation, the number of mutations in the leastloaded class increases, the ratchet has clicked. To form a new generation, individuals are sampled randomly from the previous generation, then subjected to the occurrence of mutations sampled from a Poisson distribution with mean u, and assigned probabilities of survival as (1 − s)^{k}, where k is the number of mutations that an individual carries. A new generation of N individuals is constructed by comparing the probability of survival of each individual with a pseudorandom number drawn from a uniform distribution in the interval [0, 1]. Each run was repeated several times; generally five replicates were performed to obtain the results presented in the next section.
Although this simulation procedure does not follow the fate of each mutation at a particular locus, which is extremely time consuming, it gives the same results as the multilocus stochastic simulations of Charlesworth and Charlesworth (1997) as far as the estimation of the time between clicks of the ratchet is concerned, for all parameter sets tested (results not shown).
SIMULATION RESULTS
If the loss of the leastloaded class is the determining factor in driving the ratchet (Haigh 1978), one would expect that the time between clicks of the ratchet would stay approximately constant over a range of parameter values that keep n_{0} constant. Figures 1, 2, and 3 show the simulation results for several parameter sets chosen such that n_{0} stays constant.
In Figure 1, s is kept constant at 0.015, and N changes with u to keep n_{0} constant (either 20.2 or 202). We observe that the time between clicks of the ratchet does not change significantly over an order of magnitude change in N. The increase in N seems to be compensated by the increase in u. In Figure 2, u is kept constant, at 0.1, and N changes together with s to keep n_{0} constant (with the same values as in Figure 1). Although for small values of n_{0} there seems to be no significant difference in the speed of the ratchet over an order of magnitude change in population size, for a higher value of n_{0} the difference is evident: as s becomes large (small N in the plot), the speed of the ratchet is greatly reduced. For example, for a population size of 705 individuals, with a selection coefficient of 0.08, we did not observe any click over 50,000 generations, but for a population of 19,000 individuals, with s = 0.022, the average time for a click is ~1560 generations. In Figure 3, N is kept constant and u changes with s (the mean equilibrium number of mutations, u/s, has the value 5), so that n_{0} is constant. We see that, for either small u and small s or for large u and large s, the speed of the ratchet is greatly reduced.
These results show that, as noted previously by Stephan et al. (1993), the size of the best class is not sufficient to predict the speed of the ratchet, because the ratchet can turn at very different speeds for a constant n_{0}. One observes that for the same n_{0}, keeping N constant and varying u and s so that u/s is constant, there is a value of u and s for which the time of the ratchet has a minimum, i.e., for the same n_{0}, increasing s can both slow down and speed up the ratchet (see the Ushaped curves in Figure 3). This shape is not predicted by any of the previous formulas. In contrast, such a minimum is predicted by the approximation presented here, although in the region of very small u and very small s the approximation gives lower times than the simulations. It also underestimates the time between clicks of the ratchet for small population size (or small n_{0}).
The possible reasons why this minimum is observed, and why the present approximation underestimates the time for small population size, deserve some comment. Other things being equal, decreasing s should speed up the ratchet and decreasing u should slow it down. In the case of Figure 3, both s and u are changing to keep u/s constant, so that a minimum may occur, due to the fact that the dependence of the time on the mutation rate is different from that on the selection coefficient. In the region where s is very small, so that each mutation has very little effect on fitness, u is also very small. This means that the probability of a mutation occurring is very small, and the force of mutation that drives individuals from the best class to the next class is greatly reduced, leading to a slower ratchet. In the region where the mutation rate is large, the selection coefficient is also large, so that although mutations keep appearing at a high rate, selection is so efficient in restoring the best class that a great number of individuals come from the leastloaded class, which leads to a slower ratchet.
One observes from the comparison of the theoretical formula and the simulation results that, as long as s is not extremely small or large (Figures 1, 2 and 3), the predictions seem to approximate the simulations reasonably well, especially when N is big (or n_{0} is large). If N is small, it is more difficult for the system to maintain itself close to equilibrium, because drift is dominating. In this case, an approximation based on small perturbations becomes inadequate.
We now ask how the speed of the ratchet changes with population size, for a given mutation rate and a constant selection coefficient. The simulation results and the expected times calculated with the various approximations discussed above are shown in Figures 4, 5 and 6, for different values of N, u, and s, as a function of n_{0}. The mutation rate and selection coefficient were chosen to lie in the parameter range that may be most relevant to the problem of the evolutionary degeneration of an incipient Y chromosome or neoY chromosome (Charlesworth 1978, 1996; Charlesworth and Charlesworth 1997). The values of u = 0.04 and s = 0.01 might be considered reasonable for the case of a Drosophila chromosome arm, if widely accepted estimates of mutational parameters are used (Charlesworth and Charlesworth 1998; Drakeet al. 1998). In each figure, the time between clicks of the ratchet is plotted against n_{0}, for a fixed value of u and s, so that an increase in n_{0} is solely due to an increase in N. The population size varies in the interval [1000; 13,500] in Figure 4, [1000; 10,000] in Figure 5, and [1000; 30,000] in Figure 6.
One observes that the approximation of Charlesworth and Charlesworth, contrary to what was previously thought (Charlesworth and Charlesworth 1997; Orr and Yuseob 1998), greatly underestimates the speed of the ratchet for large population sizes in the parameter range considered here. Although Stephan et al. could not establish exactly the range of validity of their two approximations, they suggested use of their Equation 8 for predicting the speed of the ratchet for the range of selection coefficients considered here. From comparison with the simulations presented, we see that their Equation 14 seems to describe the rate of change with N of the time between clicks of the ratchet better than their Equation 8, although it always underestimates the absolute time for this parameter range.
In simulations done to check the change in the ratchet's speed with different mutation rates (Figure 7), we can see that the range of parameters for which Stephan et al.'s Equation 14 gives a better quantification of the process than their Equation 8 is not only dependent on a large population size and strong selection, but is also modulated by the mutation rate. In simulations done to check the change of the ratchet's speed with different s (Figure 8) we see that, for a given N and u, as s increases, the estimated time given by their Equation 14 becomes an underestimate. In fact, Stephan et al. studied the process by dividing it in two separate phases: the establishment phase, during which the new leastloaded class reaches a value close to that of the deterministic equilibrium, and the extinction phase, during which the leastloaded class becomes extinct. Their Equation 14 is based on the assumption that the population size is large enough that the change in mean fitness of the population is sufficiently small to be approximated by small perturbations and that selection is strong enough that the process is mostly trapped in the extinction phase. This assumption is not supported by our simulations (see Figure 9). The time spent in the establishment phase is the one above the dashed line; the time spent in the extinction phase is the one below this line.
When deriving our approximation, we assumed that the net loss of mean fitness due to a click of the ratchet would be K ≈ 0.6se^{−u}, because the system did not recover its new equilibrium instantaneously. In Figure 10A, we plot the mean fitness of the population as a function of 1 − x/x_{0} for a set of simulation runs with parameters N = 10,000, u = 0.04, s = 0.01, after a click of the ratchet. As we assumed that the mean fitness of the population could be approximated as a linear function of 1 − x/x_{0}, the slope of the linear regression line plotted in the figure corresponds to the value of K in the theoretical approximation. The value assumed in the derivation (K = 5.76 × 10^{−3}) agrees well with the one from the regression. However, the agreement is not good for a population size of only 1000 (Figure 10B). This can be attributed at least to two factors: either the population size is so small that we cannot approximate the changes in mean fitness by small perturbations and/or the value of K is different from the one we are assuming. If we calculate the time by substituting the value of K from the linear fit in the theoretical expression, we find that the time obtained is still below the one measured in the simulations, so that an incorrect value for K is not the only source of the discrepancy.
DISCUSSION
Speed of the ratchet: The equilibrium size of the leastloaded class, n_{0}, is usually regarded as the chief parameter that determines the speed of the ratchet (Haigh 1978; Bell 1988; Gessler 1995; Charlesworth and Charlesworth 1998). The simulations presented here show that this parameter is not sufficient for this purpose, as noted previously by Stephan et al. (1993). For example, for an equilibrium size of the best class of ~200 we can get average times for one click of the ratchet varying from 900 to 8000 as a result of changes in selection coefficient (0.015–0.04) and mutation rate (0.075–0.2).
The parameter space considered here has been constrained to ensure n_{0} ⪢ 1, so that the Poisson distribution expected under mutationselection balance can be approximately attained. Gessler (1995) has shown that this balance will not be met for conditions under which n_{0} < 1, because selection is too weak to counter mutation pressure. In this case, the distribution of the number of mutations is not Poisson but is close to a shifted negative binomial distribution, whose parametrization allows an estimation of the rate of the ratchet.
For ⪢ 1, the approximation for the advance of the ratchet based on a diffusion equation presented here seems to make a better prediction of the time between clicks of the ratchet than the previous approximations, for moderate selection coefficients in a range that is compatible with the biological data, provided the population size is not too small, so that the stochastic fluctuations are not too violent. As noted before by Stephan et al. (1993), for intermediate selection coefficients the establishment phase and the extinction phase are blurred and a separate analysis of these two phases does not well predict the outcome. An approximation based on the assumption that the mean fitness of the population is affected solely by fluctuations of the leastloaded class
The observation that, for a constant n_{0} and constant s, we did not observe significant differences in the speed of the ratchet, over an order of magnitude change in the population size (Figure 1), suggests that n_{0}s is an important parameter, although the expression for the average time between clicks of the ratchet is not an explicit function of n_{0}s. For the parameter range considered here we observe that an increase of 10fold in n_{0}s caused a decrease of ~10fold in the speed of the ratchet.
Y and neoY chromosome degeneration: Because an incipient Y chromosome, or a neoY chromosome resulting from an autosome fusion or translocation, that fails to recombine with its homologue in the heterogametic sex is vulnerable to the ratchet, it is interesting to calculate the expected rate of its operation under the above approximation. The erosion of a protoY chromosome is very similar to the degeneration of a haploid asexual population if one replaces s by hs, where h is the dominance coefficient and s is the selection coefficient against homozygous mutations (Charlesworth and Charlesworth 1997). In the case of Drosophila, if one assumes an effective population size of males of 5 × 10^{5}, a deleterious mutation rate per Y chromosome of 0.04, and an average selection coefficient against a heterozygous mutation of 0.01, the present approximation gives a value of ~3 × 10^{25} generations for one click of the ratchet; if one sets 5 generations per year for Drosophila, this will correspond to 6 × 10^{24} years per click. However, if the mutation rate is slightly higher, say 0.06, this time would decay to 40,000 years per click. Our approximation suggests that an increase in the mean number of mutations from 4 (n_{0} = 9158) to 6 (n_{0} = 1239), just due to an increase in the mutation rate, increases the speed of the ratchet by 20 orders of magnitude. This comes from the fact that there is a more or less exponential increase of the time between clicks of the ratchet with an increase in n_{0} for a given selection coefficient. On the other hand, for a mutation rate of 0.04, if hs is ~0.005, n_{0} will be 168 and the estimated average time for a click is 697 generations. Because every click of the ratchet in a haploid population indirectly leads to the fixation of a deleterious mutation in the whole population (Charlesworth and Charlesworth 1997), this means that if the ratchet is clicking every 697 generations, in a period of 1 million years we would expect an incipient Drosophila Y chromosome to become fixed for ~7170 deleterious mutations (again assuming 5 generations per year).
The neoY chromosome system of Drosophila miranda constitutes an excellent clock to set the time scale over which the degeneration of a nonrecombining region is supposed to occur. The time estimated for the origin of the chromosomal rearrangement generating the neoY in D. miranda is ~1.25 million years ago (mya; S. Yi, personal communication). The neoY shows evidence for degeneration, and the neoX is partially dosage compensated (Steinemannet al. 1993; Steinemann and Steinemann 1998, 1999). These observations seem to suggest that, if there is a general process responsible for the degeneration of the nonrecombining segment of the genome such as Muller's ratchet, it is expected to show its signature over a time scale of the order of 10^{6} years.
In Table 1 we show the expected time for a click of the ratchet, under the approximation proposed here, for a population of half a million males for various values of u and hs. We observe that, for these values for the average effect of a mutation, the time for a click of the ratchet becomes biologically irrelevant when n_{0}hs goes above 15. If the average heterozygous effect of a nonlethal deleterious mutation is of the order of 1% (Charlesworth and Hughes 1999), under the above approximation, for n_{0} < 1500 the ratchet may play a role in the degeneration of the neoY, but its rate is probably too small to explain the degeneration observed if n_{0} > 1000. For the ratchet to be the main process causing the degeneration of the neoY, n_{0} has to be probably <500. In the case of lethal mutations, which probably occur at a rate of ~0.0025 for an incipient Y chromosome (Fryet al. 1999) and with hs of ~2% (Crow 1993), the time for a click of the ratchet in a population of half a million chromosomes is biologically irrelevant.
Muller's ratchet is a priori more likely to be an important force in driving the degeneration of mammalian Y chromosomes, given that their effective population size is around one order of magnitude less than that for Drosophila (Charlesworth and Charlesworth 1997). There is some evidence that the evolution of the mammalian Y chromosome has been punctuated by at least four events that suppressed recombination between the X and the Y (Lahn and Page 1999), the first event having occurred ~300 mya. If the deleterious mutation rate and average effect of mutations for a mammalian protoY chromosome were the same as those estimated for the Y in Drosophila, then with a population of 5 × 10^{4}, protoY chromosomes would degenerate due to this process at an average rate of 1 click every 40 thousand generations, for a deleterious mutation rate of 0.04. Our ignorance of the value of these parameters does not allow any final conclusion about the ratchet being a leading process in the degeneration of the mammalian Y, although it seems more likely than for the Drosophila case.
Because of the assumptions we have made to derive these results, some caution has to be taken considering their implications. First we must note that we have assumed, for simplicity, that all deleterious mutations have the same effect. However, recent work has suggested that an equal effect of mutations assumption does not fit the data from mutation accumulation experiments, which are designed to measure the mutation rate to deleterious mutations and the selection coefficients against those mutations (Keightley 1996; Davieset al. 1999; Fryet al. 1999). The occurrence of many mutations with small effects, and a few with large effects, seems to be more consistent with the results. If this is the case, the ratchet is expected to turn at a much higher speed than for a single selection coefficient of the order of 1% (Gessler and Xu 1999), but each turn will cause a very small decline in mean fitness if many mutations have very low selection coefficients. It is likely that the degeneration of the Y chromosome and the evolution of dosage compensation are driven by selection to increase the activity of the X relative to the Y in males, in response to the decline in mean fitness of the Y (Charlesworth 1996; Charlesworth and Charlesworth 1997). The strength of such selection is determined by the rate of this decline and will be very weak if it is small.
We have also assumed independence of mutational effects and it has been shown that epistasis slows down the speed of this process (Charlesworthet al. 1993; Kondrashov 1994). However, if there is in fact a distribution of mutational effects with a more or less exponential shape, epistasis will not stop the ratchet (Butcher 1995).
Acknowledgments
We thank D. Gessler and E. Baake for helpful comments. I. Gordo also thanks G. McVean for helpful discussions concerning this work. I. Gordo is supported by the Gulbenkian Foundation and Program PRAXIS XXI of Portugal, and B. Charlesworth by the Royal Society.
Footnotes

Communicating editor: W. Stephan
 Received August 23, 1999.
 Accepted November 22, 1999.
 Copyright © 2000 by the Genetics Society of America