Modes of Reproduction and the Accumulation of Deleterious Mutations With Multiplicative Fitness Effects
Patsy Haccou, Maria Victoria Schneider


Mutational load depends not only on the number and nature of mutations but also on the reproductive mode. Traditionally, only a few specific reproductive modes are considered in the search of explanations for the maintenance of sex. There are, however, many alternatives. Including these may give radically different conclusions. The theory on deterministic deleterious mutations states that in large populations segregation and recombination may lead to a lower load of deleterious mutations, provided that there are synergistic interactions. Empirical research suggests that effects of deleterious mutations are often multiplicative. Such situations have largely been ignored in the literature, since recombination and segregation have no effect on mutation load in the absence of epistasis. However, this is true only when clonal reproduction and sexual reproduction with equal male and female ploidy are considered. We consider several alternative reproductive modes that are all known to occur in insects: arrhenotoky, paternal genome elimination, apomictic thelytoky, and automictic thelytoky with different cytological mechanisms to restore diploidy. We give a method that is based on probability-generating functions, which provides analytical and numerical results on the distributions of deleterious mutations. Using this, we show that segregation and recombination do make a difference. Furthermore, we prove that a modified form of Haldane’s principle holds more generally for thelytokous reproduction. We discuss the implications of our results for evolutionary transitions between different reproductive modes in insects. Since the strength of Muller’s ratchet is reduced considerably for several forms of automictic thelytoky, many of our results are expected to be also valid for initially small populations.

ONE of the main dilemmas of evolutionary biology is that for all else being equal asexual populations have a twofold fitness advantage over their sexual counterparts (Williams 1975; Maynard Smith 1978). Thus, whenever the two reproductive strategies compete, the elimination of the sexual mode of reproduction is expected, unless there are factors that counterbalance its disadvantages. Nevertheless, most eukaryotes reproduce sexually (Bell 1982; Dybdahl and Lively 1995). Therefore, several theories have been developed to explain the maintenance of sexual reproduction (see Kondrashov 1993; Hurst and Peck 1996; Pecket al. 1998; Westet al. 1999).

One of the main theories that are currently used is the “deterministic deleterious mutation theory,” proposed by Kondrashov (1982, 1984). He showed that in large populations segregation and recombination may lead to a lower mutation load if there are synergistic interactions between deleterious alleles. This has led to much theoretical work on epistatic effects of such mutations (e.g., Charlesworth 1990; Otto and Feldman 1997; Barton and Charlesworth 1998). Empirical results, however, seem to indicate that it is doubtful whether such interactions indeed occur. It seems more likely that in many species effects of mildly deleterious mutations are more or less multiplicative (Willis 1993; Elena and Lenski 1997; Otto 1997; but see Riveroet al. 2003).

In the absence of epistasis, recombination has no effect on expected viability of females in sexual populations with equal male and female ploidy. Furthermore, their expected mutational load is then equal to that of clonally reproducing females. Probably because of this, the theoretical study of the evolutionary effects of deleterious mutations with multiplicative fitness effects was neglected. This may, however, have been premature. For instance, Siller (2001) showed that even without synergy, sexual reproduction does reduce the mutation load if there is a stronger selection against deleterious mutations in males than in females. This occurs, for instance, if females prefer to mate with males that have a lower mutation load. His results indicate that multiplicative fitness effects should not be ignored a priori.

In previous studies only two modes of reproduction were considered, i.e., diplodiploid (or haplohaploid) sexual and clonal asexual reproduction. However, many other modes of reproduction exist. An alternative sexual reproductive mode is, for instance, arrhenotoky, where males are produced parthenogenetically and are haploid, whereas females are diploid and produced from fertilized eggs. This occurs in many insects, especially in the Hymenoptera (e.g., Normark 2003). Since the difference in ploidy will lead to different selection pressures on males and females, similar effects as found by Siller (2001) can be expected to occur. Likewise, asexual reproduction is not always clonal, but involves meiosis and recombination in many organisms (so-called automictic thelytoky). Thus, depending on the cytological mechanism, recombination rate can affect expected viability in asexuals. As we show, this also occurs when fitness effects of deleterious alleles are multiplicative.

We present a method for calculating long-term expected viabilities when there are multiplicative effects, based on probability-generating functions, which can be used to derive analytical and numerical results on expected viabilities. Our method resembles the one introduced by Dawson (1999), based on cumulant-generating functions. However, he considers only one-dimensional generating functions, corresponding to systems where both males and females are effectively haploid and, as a consequence, segregation and recombination are equivalent.

We study asymptotic expected viabilities of populations with several sexual and asexual modes of reproduction, which are all known to occur in insects (e.g., Normark 2003), but are also used by other organisms (e.g., Suomalainenet al. 1987). This leads to surprising results. For instance, it appears that recombination is disfavored for some modes of sexual reproduction whereas it is favored in some asexual cases. Furthermore, we prove that a modified form of Haldane’s principle, which states that long-term mutation load depends only on the mutation rate and not on the selection pressure (Haldane 1937), holds for thelytokous reproduction in general.


We consider the following reproductive modes:

  • Diplodiploidy (amphimixis): This is the most well-known mode of sexual reproduction, where males and females are both diploid and produced from fertilized eggs.

  • Haplodiploidy: There are two different possibilities: arrhenotoky, where males are haploid and produced from unfertilized eggs and females are diploid and produced from fertilized eggs, or paternal genome elimination (PGE), where eggs are all fertilized but males do not pass their paternal genome to their offspring. In the latter case there are again two possibilities: males are somatically diploid but in sperm only the maternal genome is retained, or males are somatically haploid; i.e., their paternal genome is destroyed before embryogenesis.

  • Thelytoky: Females produce only daughters, from unfertilized eggs. There are two types of thelytokous reproduction: apomixis—clonal reproduction, which can be diploid or polyploid—and automixis, a meiosis that takes place after which diploidy is restored. Restoration of diploidy in the case of automixis can happen through six different cytological mechanisms, which are described by Suomalainen et al. (1987; see also White 1973). We use their terminology. The genetic consequences of four of the mechanisms are illustrated in Figure 1. In all cases the first three steps are the same and identical to the procedure in a normal meiosis: chromosome duplication is followed by recombination, which results in four chromosome sets. Although no real gametes are produced, we refer to these chromosome sets as “gametes” for convenience. The different chromatids are labeled with letters a′-d′ and the gametes with letters a-d in Figure 1. Gametes and chromatids are characterized by the origin of their centromere. The gametes a and b are homologous, and so are c and d. Zygotes are formed by duplication and/or combinations of the gametes. The mechanisms differ in the way this is done:

    • Gamete duplication: The cleavage nuclei fuse, or the halves of the divided chromosomes of the cleavage nuclei remain in the same nucleus. The chromosome sets are all duplicated and from the resulting pairs one is randomly selected to become the zygote.

    • Terminal fusion: The second polar nucleus fuses with the egg nucleus. One of the pairs of gametes (a, b) and (c, d) is randomly chosen to form the zygote.

    • Central fusion: The two central polar nuclei fuse. The pairs (a, b) and (c, d) act as two parents that produce the zygote together, so there are four possible combinations. One of those is randomly selected.

    • Second meiotic spindles fusion: This mechanism was described by Narbel-Hofstetter (1964) for Apterona helix. The result is that there are three possible zygote genotypes that occur with unequal probabilities. Suomalainen et al. (1987) refer to this as “mechanism D” and list chances 1/6, 4/6, 1/6 rather than the ones we use (1/4, 1/2, 1/4; cf. Figure 1) for the three zygote genotypes. This is, however, due to an erroneous interpretation of the cytological mechanism of the process described by Narbel-Hofstetter (L. W. Beukeboom and L. P. Pijnacker, personal communication).

    There are two additional automictic mechanisms, which are equivalent to apomixis with respect to their genetic results:

    • Gonoid thelytoky: The parthenogenetic egg undergoes two meiotic divisions, and the chromosomes pair at the second meiotic prophase.

    • Premeiotic doubling: A premeiotic doubling of chromosome number is reduced through meiosis.

Figure 1.

—Schematic of the genetic effects of thelytokous reproduction with four different automictic mechanisms. In all cases the first three steps are the same as in standard meiosis: chromosome duplication is followed by recombination. The mechanisms differ in the way the zygotes are formed from the resulting chromosome sets. As an example we show what may happen to one pair of homologous chromosomes (represented by blocks); the circles represent deleterious alleles. 1, parent genotype; 2, chromosome duplication (the chromatids are labeled a′-d′; a′ and b′ share a centromere, and so do c′ and d′); 3, recombination (the resulting chromosome sets are labeled a-d; the centromeres of a and b have the same origin, and so do those of c and d; and in this case recombination has led to an exchange of the deleterious allele between chromatids b′ and c′); 4, zygote genotypes and their expected frequencies for the different mechanisms.


Throughout we make the following assumptions:

  1. The population size is large, so that Muller’s ratchet does not operate.

  2. Generations are discrete, with the relevant steps in the life cycle occurring in the order: selection-mutation-recombination and reproduction.

  3. Fitness effects of deleterious mutations are multiplicative with a selection coefficient that depends on the ploidy of the individual.

  4. The number of mutatable loci is so large that the probability of homozygosity due to new mutations or segregation in sexual populations can be ignored.

  5. Mutation rate is constant. The numbers of new mutations per haploid genome are independent Poisson distributed with expectation λ/2.

  6. Loci are unlinked (see below).

  7. In sexual populations there is random mating.

Assumptions i, iv, v, and vii are as found in the standard literature (Kondrashov 1982, 1984; Siller 2001). The order of the steps in assumption ii is natural for haplodiploids as well as diploids if it is assumed that viability selection occurs. However, Kondrashov (1984) considers organisms with selection at the haploid stage, with the life cycle mutation-mating-recombination-selection. In his model recombination is synonymous to reproduction, as is also true for the models where haplotypes are considered.

Note that the formation of the gametes, exemplified in the top half of Figure 1, is the same in sexual females as in automictic ones. In sexual females a′ and b′ correspond to the maternally derived chromatids and c′ and d′ to the paternally derived ones. If we denote deleterious alleles by a “1” and wild-type alleles by a “0” (a notation that we use throughout), the configuration at a certain locus on the chromatids can be denoted by a set of four binary numbers, indicating the allele type at that locus on, respectively, a′, b′, c′, and d′. For instance, in the example of Figure 1 the initial chromatid state at the locus containing the deleterious allele is (0, 0, 1, 1). Similarly, after recombination has taken place the gamete configuration at a locus can be represented by such a set of binary numbers, in the example of Figure 1 by (0, 1, 0, 1). Table 1 shows all the six possible gamete states for loci at which the parent is heterozygous. We denote the total number of occurrences of gamete state i at loci with initial chromatid state (1, 1, 0, 0) by ni and those at loci with initial chromatid state (0, 0, 1, 1) by mi (i = 1,..., 6), so in Table 1 yi = ni + mi. Assumption vi implies that the numbers {n1,..., n6} and {m1,..., m6} are independent and multinomially distributed (see Equation B3). The recombination model determines the parameters of those distributions.

View this table:

Gamete configurations of deleterious mutations

There are many different recombination models (see, e.g., Lange 2002, for an overview), but some general considerations will hold for most of them. Let r be a measure of recombination rate, where r = 0 corresponds to absence of recombination. If the initial chromatid state at a certain locus is (1, 1, 0, 0) and there is no recombination, the gamete configuration at that locus will be of type 1. States 2-5 can all be produced by exchanging alleles between one pair of nonhomologous chromatids, whereas to reach state 6 two such changes are required. Thus, it is reasonable to assume that with this initial condition the probabilities of reaching gamete states 2-5 are all equal and initially (for low values of r) lower than the probability of ending up in state 1 but higher than the chance of gamete state 6. Furthermore, for r = 0 the probability that the gamete state is 1 equals one and as r becomes infinitely large all states have equal probability. We denote the probability of gamete state 1 if the initial chromatid state is 1 by π1(r), that of states 2-5 by π2(r), and that of state 6 by π3(r). It appears that the exact form of these probabilities does not affect the results very much. Here, we use the following expressions, π1(r)=16+13e32r+12er,π2(r)=1616e32r,π3(r)=1π1(r)4π2(r), (1) which were derived from a continuous-time Markov chain model for crossing over (see appendix a). If the initial chromatid state is (0, 0, 1, 1) the roles of states 1 and 6 are reversed, so in that case the probability of gamete state 1 is π3(r) and that of state 6 is π1(r).


For sexual diplodiploid populations or populations with (diploid) apomixis it is a well-known fact that the long-term expected viability depends only on the mutation rate, and when the mutation distribution is Poisson, as in our model, it equals e (e.g., Bürger 2000). The same result holds for the automictic mechanisms that have the same genetic consequences as apomixis, gonoid thelytoky, and premeiotic doubling. Furthermore, it can be shown to hold for PGE with male somatic diploidy and for polyploid apomixis as well.

To calculate expected viabilities for the other, more complicated, situations we use so-called probability-generating functions (PGFs) of the numbers of deleterious mutations. The generating function of a random variable x is defined by F(z)=E[zx], (2) where E[·] denotes the expectation over x, and z is a variable that can assume values between zero and one (see, e.g., Feller 1968). The generating function F(z) completely determines the distribution of x and can be used to calculate its moments, such as, for instance, expectation and variance. We derive recursion equations for the probability-generating functions of the numbers of deleterious mutations per genome in successive generations. Numerical iteration of such equations until the outcome is stable gives the PGFs of the stable distributions of these numbers, which can be used to calculate expected viabilities. We give an overview of the derivation of the recursion equations in the main text. Details of calculations are given in appendix b.

Haplodiploid reproduction: It follows from assumption iv that in females deleterious mutations occur only in heterozygous form. Genotypes of females are represented by two numbers (n, m) where n is the number of deleterious mutations on the maternal and m that on the paternal chromosomes. Genotypes of somatically haploid males are represented by one such number, denoted by k. It is easily seen that for our purposes arrhenotoky is equivalent to PGE with somatically haploid males. Thus the results in this section hold for both of these reproductive modes.

Populations are characterized by two probability-generating functions of the genotypes of females and males at the start of the tth generation, just before selection, Gt(z1,z2)=Et(z1nz2m),z1,z2[0,1],Ht(z)=Et(zk),z[0,1], (3) where Gt is a multivariate PGF (the two-dimensional generalization of Equation 2; see, e.g., Johnsonet al. 1997). Assumption iii implies that females of type (n, m) have survival chance (1 - hs)n+m, whereas males of type (k) survive with probability (1 -σs)k, where hs (the selection coefficient per heterozygous locus in females) and σs (the selection coefficient per deleterious allele in males) lie between 0 and 1. As a result the PGFs of female and male genotypes after selection become Gt(z1,z2)=Gt((1hs)z1,(1hs)z2)Gt((1hs),(1hs)),Ht(z)=Ht((1σs)z)Ht(1σs). (4)

Assumption v implies that mutation changes the PGFs to Gt(z1,z2)=e(λ2)(z1+z22)Gt(z1,z2),Ht(z)=e(λ2)(z1)Ht(z). (5)

Recombination occurs only in the diploid females. The possible configurations of deleterious alleles in the four types of gametes a-d shown in Table 1 can be used to derive the genotypes of the offspring. The number of deleterious mutations in a haploid son is simply the sum of the deleterious mutations on the female gamete. We assume that all gametes have equal probability of being used, so the PGF of the males in the next generation equals Ht+1(z)=14E[zy1+y3+y4]+14E[zy1+y2+y5]+14E[zy2+y4+y6]+14E[zy3+y5+y6]. (6)

In appendix b it is shown that this leads to Ht+1(z)=12(Gt((1ρ(r))z+ρ(r),ρ(r)z+(1ρ(r)))+Gt(ρ(r)z+(1ρ(r)),(1ρ(r))z+ρ(r))), (7) where ρ(r)=π3(r)+2π2(r). (8)

Since mating is random (assumption vii) and males are haploid, the PGF of the daughters is found by multiplying this function with the PGF of males, i.e., Gt+1(z1,z2)=12Ht(z2)×(Gt((1ρ(r))z1+ρ(r),ρ(r)z1+(1ρ(r)))+Gt(ρ(r)z1+(1ρ(r)),(1ρ(r))z1+ρ(r))). (9)

Combination of all the steps finally gives the recursion equations: Gt+1(z1,z2)=12e(λ2)(z1+z22)Ht((1σs)z2)Ht(1σs)Gt(1hs,1hs)×(Gt((1hs)(ρ(r)+(1ρ(r))z1),(1hs)(ρ(r)z1+(1ρ(r))))+Gt((1hs)(ρ(r)z1+(1ρ(r))),(1hs)(ρ(r)+(1ρ(r))z1))),Ht+1(z)=12e(λ2)(z1)1Gt(1hs,1hs)×(Gt((1hs)(ρ(r)+(1ρ(r))z),(1hs)(ρ(r)z+(1ρ(r))))+Gt((1hs)(ρ(r)z+(1ρ(r))),(1hs)(ρ(r)+(1ρ(r))z))). (10)

Automictic reproduction: Gamete duplication leads to homozygosis after one generation, and recombination does not have any effect on the form of the distribution of the numbers of deleterious mutations per individual. It can be shown easily that in this case the stable distribution of the number of mutations per individual is Poisson with mean λ/(2s), where λ is the mutation rate per diploid genome and s the selection coefficient per homozygous locus. As a consequence, the expected viability becomes e-λ/2.

To study the population genetic effects of the other mechanisms it is mathematically convenient to distinguish two chromosome sets, although we cannot speak of paternal and maternal chromosomes in this situation. Thus, there are now three types of loci containing deleterious alleles, (1, 0), (0, 1), and (1, 1), and genotypes are represented by three numbers (n, m, k), where n and m represent the numbers of deleterious mutations on the two chromosome sets at the heterozygous loci and k is the number of loci that are homozygous for deleterious mutations. Note that n and m have identical marginal distributions, but are not independent. Since there are only females we have to deal with only one PGF. The distribution of (n, m, k) just before selection in generation t is described by a three-dimensional PGF: Ft(z1,z2,z3)=Et(z1nz2mz3k),z1,z2,z3[0,1]. (11)

Similar to the previous case, selection changes this function to Ft(z1,z2,z3)=Ft((1hs)z1,(1hs)z2,(1s)z3)Ft((1hs),(1hs),(1s)), (12) and subsequent mutation gives Ft(z1,z2,z3)=e(λ2)(z1+z22)Ft(z1,z2,z3). (13)

The consequences of recombination and reproduction differ for terminal fusion, central fusion, and second meiotic spindles fusion.

If the number of homozygous loci in the parent is k, the different gamete combinations that can occur with terminal fusion are a with b, which (see Table 1) gives zygotes with genotype (y3 + y4, y2 + y5, k + y1), and c with d, leading to zygotes with genotype (y2 + y4, y3 + y5, k + y6). Each of these combinations occurs with chance 1/2. As a consequence Ft+1(z1,z2,z3)=12E[z1y3+y4z2y2+y5z3k+y1]+12E[z1y2+y4z2y3+y5z3k+y6], (14) which gives Ft+1(z1,z2,z3)=12Ft(2π2(r)(z1+z2)+π1(r)z3+π3(r),2π2(r)(z1+z2)+π3(r)z3+π1(r),z3)+12Ft(2π2(r)(z1+z2)+π3(r)z3+π1(r),2π2(r)(z1+z2)+π1(r)z3+π3(r)z3,z3) (15) (see appendix b).

For central fusion there are four possible combinations (see Table 1 and Figure 1): a with c gives zygotes with genotype (y1 + y3, y2 + y6, k + y4); a with d gives zygotes of type (y1 + y4, y5 + y6, k + y3); b with c gives zygotes of type (y1 + y5, y4 + y6, k + y2); b with d gives zygotes of type (y1 + y2, y3 + y6, k + y5). Since gametes are combined at random, these events all have probability 1/4. It can be shown that this leads to Ft+1(z1,z2,z3)=Ft((π1(r)+π2(r))z1+(π2(r)+π3(r))z2+π2(r)(1+z3),(π2(r)+π3(r))z1+(π1(r)+π2(r))z2+π2(r)(1+z3),z3). (16)

Finally, for second meiotic spindles fusion, there are three possible combinations: a with a gives zygotes with genotype (0, 0, k + y1 + y3 + y4); d with d gives zygotes of type (0, 0, k + y3 + y5 + y6); b with c gives zygotes with genotype (y1 + y5, y4 + y6, k + y2). The first and second combinations occur each with chance 1/4, and the third one with chance 1/2 (see Figure 1). This means that Ft+1(z1,z2,z3)=14Ft((1ρ(r))z3+ρ(r),ρ(r)z3+(1ρ(r)),z3)+14Ft(ρ(r)z3+(1ρ(r)),(1ρ(r))z3+ρ(r),z3)+12Ft((π1(r)+π2(r))z1+(π2(r)+π3(r))z2+π2(r)(1+z3),(π2(r)+π3(r))z1+(π1(r)+π2(r)+π2(r))z2+π2(r)(1+z3),z3), (17) with ρ(r) defined as in (8).

For each of the mechanisms, combination of the relation between Ft+1 and Ft with (12) and (13) leads to a recursion equation for the PGF. Their derivation is completely analogous to that of (10).


It is a well-known and very general result that in populations with clonal asexual reproduction and no epistasis the expected viability does not depend on the strength of the selection against deleterious mutations, but only on the mutation rate (Haldane 1937). This result, known as Haldane’s (or the Haldane-Muller) principle does not depend on the distribution of new mutations per genome. We prove that a similar principle holds for thelytokous reproduction: when hs is fixed, the value of the selection coefficient s does not affect expected population viability.

To illustrate the proof of Haldane’s principle, we give an outline for terminal fusion. The proof is generalized straightforwardly to the other mechanisms and to mutation distributions other than the Poisson (see appendix c). First, we introduce some new notation. Note that the shape of the PGF Ft(z1, z2, z3) depends on the values of the parameters λ, hs, s, and r. For clarity of the proof, we make the dependency on s explicit and write Ft(z1, z2, z3; s).

Filling in 0 for the values of z1, z2, and z3 in Equation 15 gives Ft+1(0,0,0;s)=12Ft(π3(r),π1(r),0;s)+12Ft(π1(r),π3(r),0;s). (18)

Substitution of (12) and (13) gives Ft+1(0,0,0;s)=CFt((1hs)π3(r),(1hs)π1(r),0;s)+Ft((1hs)π1(r),(1hs)π3(r),0;s)Ft(1hs,1hs,1s;s), (19) where C = ½e(λ/2)(π3(r)+π1(r)-2), and so, if we denote the stable PGF by F, in the long run the expected fitness is F(1hs,1hs,1s;s)=CFt((1hs)π3(r),(1hs)π1(r),0;s)+Ft((1hs)π1(r),(1hs)π3(r),0;s)F(0,0,0;s). (20)

In appendix c we prove that the conditional probabilities Pr[n = x, m = y|k = 0] (x = 0, 1,...; y = 0, 1,...) do not depend on s. This implies that we can write F(x,y,0;s)=i,jxiyjPr[n=i,m=jk=0]Pr[k=0;s], (21) so Pr[k = 0; s] cancels in the fraction on the right-hand side of (20), which, as a consequence, does not depend on s.


Iterations of PGFs were performed with Mathematica 4.0. To this end the PGFs were discretized on a grid of z-values with intervals of 0.05 (and, where necessary, to improve the accuracy of the result, 0.005). In all cases, the functions converged. Iterations continued until the sum of the absolute differences between successive values of the PGFs was <10-7. This usually occurred within 2000 iterations. Figure 2 shows the stable results with λ= 1 and hs = 0.01, which are the values that are most commonly used in the literature.

Figure 2.

—Expected viabilities of females for all reproductive modes. Parameter values are λ= 1, hs = 0.01, and for haplodiploidy σs = 0.05. Haplodiploidy here implies arrhenotoky or PGE with somatically haploid males. Mechanisms equivalent to cloning or amphimixis are PGE with somatically diploid males, gonoid thelytoky, and premeiotic doubling.

Populations with haploid males have the highest expected viabilities. Viabilities increase with the selection pressure on the males (when σs = 1, the expected viability is 0.980 for all r, results not shown). In thelytokous populations, the expected viability is the lowest with apomixis (regardless of the ploidy level) or automixis with gonoid thelytoky or premeiotic doubling. Then central fusion follows, which does increasingly better with higher recombination chances. Terminal fusion attains the same viability as gamete duplication at r = 0, with only a slight decrease as r increases. The expected viability for second meiotic spindles fusion lies in between central and terminal fusion for small values of r, but at high recombination rates it nearly reaches the same value as gamete duplication. The expected viabilities attained with amphimixis or PGE with somatically diploid males are just as low as with apomixis.

Recombination affects only the expected viabilities of populations with terminal, central, or second meiotic spindles fusion. It can be seen from Figure 1 that when there is no recombination (r = 0) terminal fusion is equivalent to gamete duplication, whereas central fusion corresponds to clonal reproduction. Recombination is disadvantageous with terminal fusion, since it reduces the chances on homozygous loci. This decreases the effectiveness of selection against deleterious mutations. With central fusion, however, recombination enhances the chances of homozygosity and therefore improves selection. It can be seen from Figure 2 that, whereas the disadvantage for terminal fusion is only slight, there is a huge advantage of recombination when central fusion is used. As r tends to infinity, the right-hand sides of Equations 15 and 16 both converge to the same expression, which implies that for completely free recombination the two mechanisms are equivalent. Their expected viabilities converge to 0.601. Second meiotic spindles fusion can be considered as a mixture of gamete duplication and central fusion. Here, recombination provides only a very slight advantage. The limit value of the expected viability is in this case 0.605.

Figure 3.

—Illustration of the effect of recombination on the distribution of male genotypes. Dashed lines, distributions of numbers of deleterious mutations on gametes formed from maternal and paternal chromosomes; solid lines, distributions of number of deleterious mutations in males. (A) Without recombination. (B) With recombination.

Recombination rate does not affect expected viability in organisms with amphimixis or PGE with male diploidy. Contrary to the standard results for sexual reproduction, however, Figure 2 illustrates that it does have an effect on the expected viability of females when males are haploid. Further numerical study indicates that the magnitude of this effect decreases with σs (results not shown). The reason for the adverse effect of recombination is illustrated in Figure 3: because of the different selection pressures on males and females, the distributions of deleterious mutations on paternal and maternal chromosomes in females are different. Genotypes of sons are formed by random selection from the gametes formed from these two sets of chromosomes. The effect of recombination on these distributions is to make their expectations more similar and their variances broader, with a net effect of making the variance of the mixture distribution smaller. Thus, the variance in numbers of deleterious mutations in males is decreased by recombination and this makes selection less effective. This ultimately leads to a higher mutation load in females.


Our main conclusion is that, when reproductive modes other than the “standard” sexual and asexual ones are considered, both recombination and segregation may be found to affect mutational load in large populations, even in the absence of epistasis.

Whereas recombination is neutral in sexual systems when males as well as females are diploid, it is disadvantageous as soon as males are somatically haploid. Our explanation for this, illustrated in Figure 3, may also account for other phenomena, such as the negative effect of recombination on expected viability in systems with strong synergistic or positive epistasis in diplodiploid systems (Otto and Feldman 1997; Philipset al. 2000). In that case the marginal distributions of numbers of deleterious mutations on paternal and maternal genomes are identical, but it appears that under specific conditions recombination may still lead to a lower variance of their mixture. This is a subject of further study.

Recombination affects mutational load in asexuals in three of the six known cytological mechanisms of thelytoky. In one of these cases increased recombination strongly decreases mutational load and is therefore very advantageous. This mechanism, central fusion, occurs in, e.g., thelytokous strains of the parasitic wasp Venturia canescens (Beukeboom and Pijnacker 2000; Schneideret al. 2002) and in the cape honeybee Apis mellifera capensis (Tucker 1958; Verma and Ruttner 1983). It would be interesting to examine whether there are indeed high recombination rates in such species. A technical difficulty is, however, that this cytological mechanism in the long run leads to homozygosity at loci far removed from the centromere, which hampers the estimation of recombination rates. In the other cases, terminal fusion and second meiotic splindles fusion, recombination has only a slight effect.

Expected viabilities are highest in haplodiploids, provided that males are somatically haploid and under a stronger selection pressure than the heterozygous females. This shows that segregation can provide an advantage even without recombination. It also suggests that deleterious mutations may have played an important role in the origin of arrhenotokous systems (see also Goldstein 1994).

Our findings indicate that in general haplodiploidy and automixis should be more successful than amphimixis or apomixis. This, however, is strongly at odds with the empirical evidence. Automixis and apomixis both appear to be associated with recent lineages, indicating that in the long run they are relatively unsuccessful against sexual reproduction. Furthermore, haplodiploidy is less abundant than diplodiploidy. Thus, the model appears to be unable to explain the general pattern concerning reproductive modes, and alternative explanations may have to be considered (see Hurst and Peck 1996; Pecket al. 1998; Westet al. 1999).

Our predictions can also be compared to empirical results on frequencies of evolutionary transitions between different reproductive modes. Up to now, the only study that is complete enough to do so is Normark’s (2003) review on insect reproduction. He found a relatively large number of evolutionary transitions from amphimixis to thelytoky, which is in agreement with our predictions. However, most of these transitions appear to be to apomixis rather than to automixis. Figure 2 shows that it is much harder for thelytoky to invade arrhenotokous systems. This may play a role in the maintenance of arrhenotoky and is consistent with the findings of Normark (2003), that transitions from arrhenotoky to thelytoky are much less frequent than those from amphimixis. (We disregard transitions caused by Wolbachia infection, where horizontal transmission is involved.)

A full study of the evolutionary consequences of our findings, however, should take differences in fecundities of sexuals and asexuals into account. The “twofold cost of sex” argument is based on the assumption that these are equal. Empirical evidence suggests that this is not always true: thelytokous females may initially have a much lower fecundity (Suomalainenet al. 1987). Further, initially individuals with an alternative reproductive mode will have the same mutational load as the resident population. We are presently studying the evolutionary consequences of such factors. The findings presented in this article, however, already indicate that all the automictic systems can sustain a considerable reduction in fecundity and still outcompete amphimixis, since the expected viabilities that are shown in Figure 2 correspond to the relative growth rates if fecundities of thelytokous females are reduced by half compared to sexual ones.

We assumed that offspring that are homozygous for deleterious alleles, although they may suffer a larger selection pressure, are nevertheless viable. Alternatively, homozygosity for deleterious mutations may lead to purging at the zygotic stage, so that no or few homozygous offspring are produced. This may provide an explanation for the initial reduction in fecundity of asexual females mentioned above. Purging may provide a considerable advantage for thelytokous reproduction. For instance, since gamete duplication leads to immediate homozygosity (see Figure 1), complete purging implies that there are no deleterious mutations in the offspring and the expected viability becomes one. The evolutionary consequences of the trade-off between reduced fecundity and increased viability in connection with purging are the subject of further study.

We proved that as long as the strength of selection on heterozygous loci hs remains constant, the selection coefficient on homozygous loci s does not affect the expected viability in thelytokous automictic systems. This generalized version of Haldane’s principle is valid regardless of the mutation distribution. It implies, e.g., that a population with a selection coefficient s = 0.1 and a coefficient of dominance h = 0.1 will in the long run have the same expected viability as a population with s = 0.01 and h = 1. Note, however, that contrary to the situation for apomictic systems hs does affect expected viabilities in automictic systems with terminal fusion, central fusion, or second meiotic spindles fusion. For instance, when hs is increased from 0.01 to 0.02, the limit value of the expected viability for r ∞ drops with → and 0.006 for terminal fusion and central fusion with 0.001 for second meiotic spindles fusion.

Iteration of recursion equations such as given in Equation 10 provides estimates of the stable PGFs of genotypes. Therefore, our method provides complete information of their stable distributions, without the need for simulations. As demonstrated, it can be used to derive analytical as well as numerical results. The calculation of higher moments of the distributions of genotypes, however, involves the derivatives of the PGFs. To estimate these accurately, it is best to derive recursion equations for the derivatives themselves. This can be done straightforwardly. Our method can easily be extended to examine other multilocus models with multiplicative allelic effects. When epistasis occurs, however, it is no longer possible to derive linear recursion equations for the PGFs.

To study the effects of recombination we used a model with unlinked loci, where chiasmas were assumed to occur according to a Markov process (see appendix a). These assumptions are not so realistic for small positive values of the recombination rate r. However, since any realistic model must give the same results for r = 0 (no recombination) as well as in the limit as r tends to infinity (completely free recombination), we expect that qualitatively the conclusions remain the same when linkage is taken into account.

We assumed that population sizes are large. It has been argued repeatedly that small asexual populations suffer from Muller’s ratchet, due to random loss of the individuals with the lowest numbers of deleterious mutations (Muller 1964; Maynard Smith 1978). As long as deleterious alleles occur in only heterozygous form, however, the class with zero mutations can be restored in thelytokous populations with terminal fusion, central fusion, and second meiotic spindles fusion (see Figure 1 and Table 1). Muller’s ratchet starts to operate only once there are individuals in the population that are homozygous for those alleles. Thus, the strength of the ratchet is reduced considerably, and our results are expected to remain valid for initially small populations. (This is a subject of further study). This finding reinforces our argument that it is important to include alternative modes of reproduction in studies of the evolution of sex.


We assume that the two chromatids that participate in a chiasma are chosen independently from chiasma to chiasma; i.e., there is no chromatid interference (see, e.g., Lange 2002). Consider the state at one specific locus. At each chiasma only two alleles can switch places. For instance, if the current state is 1 and chromatids a′ and c′ participate, the state becomes 2 (numbering of states as in Table 1). This happens with probability 1/4. Continuing in this way it can be derived that the matrix of transition probabilities between the states given that a chiasma occurs is M=01414141401412000141401200141400120141400012140141414140. (A1)

If, furthermore, we assume that the distance between chiasmas is exponentially distributed, the chromatid state at a locus will change according to a Markov process. The parameter of the exponential distribution is the probability per length unit that a chiasma occurs. Let r denote this parameter times the distance between the centromere and the considered locus; then it follows from the theory of Markov chains (see, e.g., Feller 1968) that the distribution of the gamete states is p(r)=p(0)exp[r(MI)], (A2) where M is the transition matrix given in (A1) and I the identity matrix. If the initial chromatid state is 1, p(0) = (1, 0,..., 0), which leads to p(r) = (π1(r), π2(r),..., π2(r), π3(r)), where the values of πi(r) are as given in (1). If the initial chromatid state is 6, p(0) = (0,..., 0, 1), which leads to p(r) = (π3(r), π2(r),..., π2(r), π1(r)).


Haplodiploid reproduction: Selection: Let gt(n, m) represent the probability of a female of type (n, m) at the start of the tth reproduction period (before recombination) and ht(k) the probability of a male of type (k) at this time. After selection these distributions become gt(n,m)=(1hs)n+mgt(n,m)n,m(1hs)n+mgt(n,m),ht(k)=(1σs)kht(k)k(1σs)kht(k). (B1)

Using the definitions of the PGFs given in (3) gives (4).

Mutation: The assumption that new mutations occur independently of already existing ones implies that the generating function of the total number of mutations is the product of the generating functions of the numbers of new and the old mutations. (This can be derived from the fact that the expectation of the product of independent variables is the product of their expectations.) We assume that numbers of new mutations on the two sets of chromosomes in females are independent, so the PGF of the additional numbers of mutations in females is E[z1u1z2u2]=i=0j=0z1iz2j1i!j!(λ2)i+jeλ=exp[λ2(z1+z2)λ]. (B2)

Multiplication with Gt (n, m) gives the PGF for the females. An analogous derivation gives the PGF for the males.

Recombination and reproduction: From assumption vi it follows that conditionally on the total number of deleterious alleles on the maternal chromosomes n, the ni are multinomially distributed with the probabilities determined by πi(r), i.e., Pr(n1,,n6n)=n!n1!n6!(π1(r))n1×(π2(r))n2+n3+n4+n5(π3(r))n6. (B3)

Furthermore, the mi are independent of the ni, with the same type of distribution [substitute mi for ni and m for n in the equation above, and exchange π1(r) and π3(r)]. It can be shown straightforwardly that if yi = ni + mi, for arbitrary nonzero constants c1,..., c6, E[c1y1c2y2c3y3c4y4c5y5c6y6n,m]=(π1(r)c1+π2(r)(c2+c3+c4+c5)+π3(r)c6)n×(π3(r)c1+π2(r)(c2+c3+c4+c5)+π1(r)c6)m. (B4)

This follows from the independence of the ni and mi and their multinomial distribution. Therefore, conditionally on the mother’s genotype (n, m) the first expectation in (6) equals E[zy1+y3+y4n,m]=(π1(r)z+2π2(r)z+2π2(r)+π3(r))n×(π3(r)z+2π2(r)z+2π2(r)+π1(r))m=((1ρ(r))z+ρ(r))n(ρ(r)z+(1ρ(r)))m, (B5) with ρ(r) as defined in (8). Proceeding in the same way for the other expectations we find that, if n′ denotes the number of deleterious mutations in sons, E[znn,m]=12((1ρ(r))z+ρ(r))n((1ρ(r))+ρ(r)z)m+12((1ρ(r))+ρ(r)z)n((1ρ(r))z+ρ(r))m. (B6)

The PGF for the male offspring is found by taking expectations over n, and m, which gives (7). The PGF for the female offspring corresponds to E[z1nz2k] . Since n′ and k are independent, this PGF is calculated by substituting z1 for z in (B6), multiplying with z2k , and subsequently taking expectations over n, m, and k, which is the same as multiplication by Ht(z2) .

Thelytokous reproduction: We denote the genotype of the daughters by (n′, m′, k′). For terminal fusion we find, using (14), that, taking expectations over the multinomial distributions conditionally on the mother’s genotype (as was done in the calculations for haplodiploidy), E[z1nz2mz3kn,m,k]=12{(π1(r)z3+2π2(r)(z1+z2)+π3(r))n×(π3(r)z3+2π2(r)(z1+z2)+π1(r))m+(π3(r)z3+2π2(r)(z1+z2)+π1(r))n×(π1(r)z3+2π2(r)(z1+z2)+π3(r))m}×z3k. (B7)

For central fusion: E[z1nz2mz3kn,m,k]=(π1(r)z1+π2(r)(z1+z2+z3)+π3(r)z2)n×(π3(r)z1+π2(r)(z1+z2+z3)+π1(r)z2)mz3k. (B8)

For second meiotic spindles fusion: E[z1nz2mz3kn,m,k]={14(π1(r)z3+π2(r)(2z3+2)+π3(r))n×(π1(r)+π2(r)(2z3+2)+π3(r)z3)m+14(π1(r)+π2(r)(2z3+2)+π3(r)z3)n×(π1(r)z3+π2(r)(2z3+2)+π3(r))m+12(π1(r)z1+π2(r)(z1+z2+z3+1)+π3(r)z2)n×(π1(r)z2+π2(r)(z1+z2+z3+1)+π3(r)z1)m}×z3k. (B9)

Taking expectations over n, m, and k, in (B7), (B8), and (B9), and using the definition of the PGF in (11), gives (15-17).


We first prove that the conditional probabilities Pr[n = x, m = y|k = 0] (x = 0, 1,...; y = 0, 1,...) do not depend on s. Since a female with more than zero homozygous loci can never get offspring with zero homozygous loci, the following recursion equations hold in all instances, Prt+1[n,m,0]=i,jPrt[i,j,0](1hs)i+jEt[(1hs)n+m(1s)k]Pr[(i,j,0)(n,m,0)], (C1) where Pr[(i, j, 0) → (n, m, 0)] denotes the probability that a female of type (i, j, 0) produces offspring of type (n, m, 0). If hs is fixed, these chances depend on mutation, recombination, and the cytological mechanism of thelytoky, but not on s. Therefore we can write Prt+1[n,m,0]=i,jPrt[i,j,0]α(i,j,n,m,λ,ρ,hs)Et[(1hs)n+m(1s)k], (C2) which leads to Prt+1[n,m,0]n,mPrt+1[n,m,0]=i,jPrt[i,j,0]α(i,j,n,m,λ,ρ,hs)n,mi,jPrt[i,j,0]α(i,j,n,m,λ,ρ,hs)=i,j(Prt[i,j,0]k,lPrt[k,l,0])α(i,j,n,m,λ,ρ,hs)n,mi,j(Prt[i,j,0]k,lPrt[k,l,0])α(i,j,n,m,λ,ρ,hs). (C3)

Thus, since Prt[i,j,0]j,lPrt[k,l,0]=Prt[n=i,m=jk=0], (C4) we find the recursion equation Prt+1[n=i,m=jk=0]=i,jPrt[n=i,m=jk=0]α(i,j,n,m;λ,ρ,hs)n,mi,jPrt[n=i,m=jk=0]α(i,j,n,m;λ,ρ,hs), (C5) and since α(i, j, n, m; λ, ρ, hs) does not depend on s, neither do the conditional probabilities. (Note that the unconditional chances Prt[i, j, 0] do depend on s.)

In the main text we gave only the proof for a specific case. To show that it holds for arbitrary mutation distributions, we rewrite (13) as Ft(z1,z2,z3;s)=ϕ(z1,z2)Ft(z1,z2,z3;s), (C6) where ϕ(z1, z2) is the PGF of the mutation distribution. Further, a closer look at (15-17) reveals that these are all equations of the form Ft+1(z1,z2,z3;s)=iaiFt(jbijzj+βij,jcijzj+γij,z3,s), (C7) where the ai, bij, βij, γij are constants. Combination of (C6) and (12) with (C7) and filling in the values 0 for z1, z2, and z3 gives the relation Ft(0,0,0;s)=iaiϕ(jβij,jγij)×F((1hs)jβij,(1hs)jγij,0;s)Ft(1hs,1hs,1s;s), (C8) and so, for the stable PGF we find F(1hs,1hs,1s;s)=iaiϕ(jβij,jγij)×F((1hs)jβij,(1hs)jγij,0;s)F(0,0,0;s). (C9)

From (21) and the fact that the conditional probability in that equation does not depend on s we can conclude that the right-hand side of (C9) does not depend on s and so neither does the expected viability.


We thank Leo Beukeboom, Tom van Dooren, Hans Metz, Bart Pannebakker, Laas Pijnacker, Jacques van Alphen, Stuart West, Bas Zwaan, and an anonymous referee for discussion and comments on previous versions of this work. The research of M.V.S. was supported by the Netherlands Science Foundation (NWO), grant no. 809.34.007.


  • Note added in proof: An article by Hopf et al. (F. A. Hopf, R. E. Michod and M. J. Sanderson, 1988, The effect of the reproductive system on mutation load. Theor. Popul. Biol. 33: 243-265) that contains mathematical results, derived by different means, which predate those in this article, recently came to our attention. Where we focus on the selective advantages or disadvantages of different reproductive modes for partially recessive mutations they focus on the generality of Haldane’s principle (in its original form), for completely recessive deleterious mutations, and its approximate validity in the case of partial recessivity. By using recurrence relations for the probability distributions, they show that in several cases the numbers of homozygous and heterozygous loci have independent Poisson distributions. Similar results follow from our recurrence relations for the generating functions, although the parameters of their and our stable distributions differ, due to differences in the detailed assumptions about the recombination process and the attendant consequences of recombination. If we take the perspective of their article, our most intriguing result is the generalized form of Haldane’s principle for thelytokous automictic systems. Furthermore, whereas they consider only Poisson mutation distributions, our method can be generalized to other distributions and our generalization of Haldane’s principle continues to hold.

  • Communicating editor: M. W. Feldman

  • Received June 6, 2003.
  • Accepted October 20, 2003.


View Abstract