## Abstract

The accumulation of deleterious mutations plays a major role in evolution, and key to this are the interactions between their fitness effects, known as epistasis. Whether mutations tend to interact synergistically (with multiple mutations being more deleterious than would be expected from their individual fitness effects) or antagonistically is important for a variety of evolutionary questions, particularly the evolution of sex. Unfortunately, the experimental evidence on the prevalence and strength of epistasis is mixed and inconclusive. Here we study theoretically whether synergistic or antagonistic epistasis is likely to be favored by evolution and by how much. We find that in the presence of recombination, evolution favors less synergistic or more antagonistic epistasis whenever mutations that change the epistasis in this direction are possible. This is because evolution favors increased buffering against the effects of deleterious mutations. This suggests that we should not expect synergistic epistasis to be widespread in nature and hence that the mutational deterministic hypothesis for the advantage of sex may not apply widely.

THE interactions between deleterious mutations, known as epistasis, can have important effects on many evolutionary processes. Epistasis affects the mutation load of a sexual population (Kimura 1961; Kimura and Maruyama 1966), the evolution of ploidy (Kondrashov and Crow 1991; Jenkins and Kirkpatrick 1995), speciation (Wagner *et al*. 1994), and Muller's ratchet (Felsenstein 1974; Kondrashov 1994; Muller 1964). Epistasis also affects the linkage disequilibria between deleterious mutations (Barton 1995; Peters and Lively 2000). This is the basis of the mutational deterministic theory for the evolution of sex, which finds that sex can provide an advantage when interactions between deleterious mutations are synergistic (*i.e*., mutations have a stronger effect in combination than would be expected from their individual effects) (Kimura and Maruyama 1966; Kondrashov 1982, 1988, 1993; Charlesworth 1990).

This has led to a number of experimental studies to quantify the type and strength of epistasis between deleterious mutations in a variety of organisms (deVisser *et al*. 1996, 1997a,b; Elena and Lenski 1997; Elena 1999; Whitlock and Bourguet 2000; Wloch *et al*. 2001; Burch and Chao 2004). Unfortunately, epistasis is difficult to measure (West *et al*. 1998). The measurements that exist show no general patterns. Rather, some studies have shown mostly synergistic epistasis (Mukai 1969; deVisser *et al*. 1996, 1997a; Whitlock and Bourguet 2000), some have found no epistasis (Elena 1999; Wloch *et al*. 2001), some have found antagonistic epistasis (Bonhoeffer *et al*. 2004; Burch and Chao 2004; Jasnos and Korona 2007), and others have found that some mutations interact synergistically and others antagonistically, with no clear bias toward one or the other type of interaction (deVisser *et al*. 1997b; Elena and Lenski 1997). Since the experimental evidence is unclear, in this article we focus instead on a theoretical approach, asking under which circumstances synergistic or antagonistic epistasis would tend to evolve. By understanding whether evolution favors one type of epistasis, we can shed light on whether or not it is plausible that this type of epistasis is common.

Analyzing the evolution of epistasis can also help us understand the forces shaping patterns of genetic architecture. We can view epistasis as a property of the types of genetic networks in an organism: for example, more modular genetic pathways will tend to favor different types of epistasis than less modular ones (Wagner and Altenberg 1996; Liberman and Feldman 2005; Hansen *et al*. 2006). Alternatively, we can think of epistasis as a way of looking at general properties of fitness landscapes: for example, rougher and smoother landscapes correspond to different patterns of epistasis between mutation (Whitlock *et al*. 1995). Thus by studying the evolution of epistasis between deleterious mutations, we are studying whether certain types of fitness landscapes or types of genetic networks are likely to be favored.

Epistasis is caused by a variety of different mechanisms. This makes it unclear precisely what the evolution of epistasis means. Epistasis between many loci could be determined in a complex way by the states of these and possibly multiple other loci. In this article, we do not attempt to model all possibilities. Rather, we use a simple modifier locus approach, in which there are a number of “major” loci at which deleterious mutations can occur, and another “modifier” locus that controls the epistatic interactions between the major loci. This approach may be realistic in certain situations. For example, one gene may code for an enzyme that catalyzes the reaction between several others; mutations in this modifier gene make mutations in certain residues in the other interacting proteins interact more synergistically or more antagonistically. Alternatively, mutations that affect crosstalk between distinct functional modules, or that cause duplications or eliminate redundancy in a genetic network, could be modifier loci that affect the epistatic interactions between many other loci (Gu *et al*. 2003; Wagner 2005; Sanjuan and Elena 2006). However, our main goal with this particularly simple model is to focus on whether generally more synergistic or more antagonistic epistasis is favored by evolution, and by how much, rather than to accurately account for all the details of how epistasis is determined.

Our approach is related to the modifier approach of Liberman and Feldman (2005, 2006; Liberman *et al*. 2007). These authors considered two loci that control fitness in diploids, coupled to a third modifier locus that affects the epistatic interactions between them. They analyzed the fate of a rare allele at the modifier locus that changes the epistasis between the other two. This work considered more general types of epistasis than we study, but only within the context of a simple two-locus system. Others have used an alternative approach based on the theory of quantitative inheritance, which views epistasis as a contribution to variance in phenotype and asks about the evolution of this epistatic variance (Hansen and Wagner 2001a,b; Hermisson *et al*. 2003; Carter *et al*. 2005; Hansen *et al*. 2006).

After defining our model more precisely in the next section, we begin our analysis by studying the development of mutation–selection balance in a population with a given type and strength of epistasis between the major loci, initially in an asexual population and then with varying levels of recombination. We next study the evolution of epistasis by considering whether a mutation at the modifier locus that changes the epistasis between the major loci can invade a population that is initially homogeneous at the modifier locus.

## MODEL

We consider a haploid population evolving under the pressure of mutation to alleles at many loci which may interact epistatically. Suppose that there are *L* loci, each with two possible alleles. At each locus *i*, allele *A _{i}* is assumed to be more fit and mutates to the less fit allele

*a*at rate μ per individual per generation. In isolation, each allele

_{i}*a*carries a fitness cost

_{i}*s*, and we assume . That is, an individual with one deleterious mutation has fitness

*e*

^{−s}≈ 1 −

*s*.

We introduce epistasis by assuming that the fitness of an individual with *k* deleterious mutations is(1)The parameter ϵ describes the sign and strength of the epistasis. When ϵ < 0, the mutations in combination have less effect than the product of their separate effects; this corresponds to antagonistic (positive) epistasis. When ϵ > 0, the mutations in combination have a larger effect than the product of their separate effects; this is synergistic (negative) epistasis. This is illustrated in Figure 1. This model is identical to that of Lenski *et al*. (1999), with their β equal to our 1 + ϵ. It differs from earlier models where fitness is a quadratic function of *k* (Charlesworth 1990; Elena and Lenski 1997); these models present a less clear distinction between synergistic and antagonistic epistasis and can introduce artifacts for large *k* (Wilke and Adami 2001).

We assume that selection is strong compared to mutation, such that . In this regime, haplotypes with more deleterious mutations are always rare compared to haplotypes with less, so that we can neglect back mutations from *a* to *A*. In addition, because for essentially all individuals in the population, multiplicative and additive epistasis are equivalent:(2)Finally, assuming means that stochastic effects and Muller's ratchet can be ignored even in moderate-sized populations, and hence we can use deterministic rather than stochastic equations to describe the evolution of the population. This focus on moderate to larger populations makes sense, since this is where the mutational deterministic hypothesis for the advantage of sex and other important effects of epistasis apply.

## THE EPISTATIC MUTATION–SELECTION BALANCE

We begin by studying the evolution under mutation and selection of a population that is homogeneous at the modifier locus, with epistasis ϵ. We start with the asexual case and then consider the effects of recombination.

#### An asexual population:

Denote by *F _{k}* the frequency in the population of individuals with

*k*deleterious mutations. The frequency in the next generation can then be written in terms of the frequency in the previous generation, giving the recursions(3)where is the mean fitness of the population, measured after mutation.

Since we assume , we can neglect terms in the product of *s* and μ in Equation 3 and neglect *k* compared to *L*. This produces the simplified recursions(4)(5)

At steady state, this yields(6)(7)where . Note that for ϵ = 0, this reduces to the classical result in the absence of epistasis (Haigh 1978).

It is easy to calculate from Equation 7 to demonstrate that Equation 6 is satisfied, which confirms the consistency of our solution. Note that while and the ratio do not depend on ϵ, the full distribution does (Figure 2). More antagonistic epistasis leads to the maintenance of more multiple deleterious mutants in the population, and vice versa. However, these different distributions are such that the mean fitnesses of populations at different values of ϵ are identical.

#### Recombination:

To study the effects of recombination, we must assume an arrangement for the location of the various loci. Suppose that each individual has a single linear chromosome and that all of the sites *A _{i}*, which we refer to as the “major loci,” are arranged along this chromosome, with a recombination rate

*r*between each pair of adjacent loci. There is an additional locus

*M*, which we refer to as the “modifier locus,” at which different alleles confer different values of ϵ. This locus is assumed to be at the end of the chromosome, with a recombination rate

*R*between

*M*and the rest of the major loci (Figure 3).

Recombination between the modifier locus and the major loci clearly matters only when the modifier locus is polymorphic, but recombination among the major loci will have an effect on the distribution of deleterious alleles even when the population is homogeneous for a single value of ϵ. In fact, when *r* ≠ 0 (*i.e*., there is recombination between the major loci), it is no longer true that all major loci are equivalent, and the frequencies of deleterious alleles at different loci will not be equal. Instead, deleterious alleles will either tend to cluster together along the chromosome or tend to spread themselves apart. By doing so, they will either minimize or maximize the effective recombination rate among themselves.

To analyze this effect, we use an iterative approach, formally valid for small recombination rates *r*. We begin with the solution to the *r* = 0 case, where the frequencies of deleterious alleles and the linkage disequilibria are uniform across the chromosome. We write down the recursions describing the dynamics of the frequencies of these alleles for nonzero *r*. These are impossible to solve exactly, but we solve for the steady-state (now position-dependent) frequencies of the deleterious alleles by assuming the linkage disequilibria are still given by their *r* = 0 values. From these approximate frequencies, we can calculate a better approximation for the position-dependent linkage disequilibria. We can then plug this new formula for the disequilibria into the recursions for the frequencies of the alleles and solve these again to give a more accurate expression for the frequencies of the deleterious alleles. We can repeat this iterative process; for small *r* each successive iteration produces a smaller correction, so we can iterate as far as needed to get the desired accuracy. Fortunately, for small *L*μ/*s*, the recombination rate can be quite large before this approach fails; the precise conditions are described below.

To begin this analysis, we define *F _{k}*(

*i*

_{1}, …,

*i*) to be the frequency of individuals with a total of

_{k}*k*mutations at loci

*i*

_{1}, …,

*i*∈ {1, …,

_{k}*L*}. Note that , where the sum is over all possible configurations of

*i*

_{1}<

*i*

_{2}< … <

*i*. We define the linkage disequilibrium between two loci

_{k}*i*,

*j*as(8)Note that this definition of

*D*neglects terms involving mutations at other sites; these will all be higher order in

*L*μ/

*s*. Without recombination (

*i.e*., for

*r*= 0), we find from Equation 7 that to lowest order in μ/

*s*, the linkage disequilibrium is a constant,(9)Note that

*D̂*has the opposite sign to ϵ, in agreement with earlier work (Feldman

*et al*. 1980; Liberman and Feldman 2005). We expect that recombination will act to reduce the absolute value of the steady-state linkage disequilibrium for all pairs of loci, but will not change the sign (this is confirmed by simulations, not shown). Note that this means that in steady state, we expect

*D*to be at most of order (μ/

*s*)

^{2}, regardless of

*r*.

We assume a model of recombination where recombination acts separately from selection (*i.e*., selection acts only on nonrecombining sequences; our model also describes the case where recombination acts jointly with selection but *r* is small enough that we can neglect products of *s* and *r*). The recursions for the frequencies of the deleterious alleles thus become, to order (μ/*s*)^{2} (and ignoring products of *s* and μ),(10)(11)(12)These recursions cannot be solved exactly. However, as noted above, for small *r* we can find an approximate equilibrium by assuming that *D* takes its *r* = 0 value given in Equation 9. We detail this in appendix a and find(13)(14)

We can use these expressions for and (and the values of that they imply, ) to calculate a better approximation for , which will now be position dependent. We find(15)We see that |*D*| is equal to its *r* = 0 value minus a position-dependent correction due to recombination that is small when . As claimed, the effect of recombination is to reduce the absolute value of *D* without changing its sign. We can insert this corrected *D* back into our recursions to calculate a better approximation for and and repeat this iterative process until the desired accuracy is reached. The above formulas make clear that the condition required for this iterative approach to converge is ; this is what is meant by low recombination rates. In fact, whenever , the above formulas are already a good approximation to the disequilibria and the frequencies of the deleterious alleles . A comparison of these results with simulations is shown in Figure 4.

Some of the results above are valid even when is not small compared to 1. We see from Equation 13 that for , the correction to due to nonzero *r* and the position dependence of the frequencies of these single mutants are small compared to the *r* = 0 value . For small, any biological values of *r* (0 < *Lr* < 0.5) satisfy this small-*r* condition. Because recombination can only reduce |*D*|, this means that even this zeroth order iterative result gives accurate results for the average value of , and hence also for the mean fitness of the population, for any biological *r*. On the other hand, the position dependence in Equation 14 will be small compared to the average only when . If this condition is not satisfied, then the number of double mutants formed by recombination of single mutants will be comparable to the number of double mutants formed by mutation of single mutants. In this case, the average value of [and therefore ] may differ substantially from its *r* = 0 value, and there may be significant position dependence in both quantities.

We can also calculate how recombination affects the mean fitness of the population. From the recursion Equation 10, we find that with recombination , where the difference in equilibrium mean fitness from the *r* = 0 value is(16)This correction is at most order , so for any biological value of *r*. For we can approximate by its *r* = 0 value to obtain(17)This analytical result for agrees with simulations (data not shown).

Note that has the opposite sign to *D* and therefore the same sign as ϵ. Thus, recombination increases the equilibrium mean fitness of populations with synergistic epistasis and decreases the mean fitness of populations with antagonistic epistasis. In this sense, *synergistic epistasis favors recombination*. This agrees with earlier work on the mutational deterministic hypothesis for the advantage of sex (Kimura and Maruyama 1966; Kondrashov 1982, 1988, 1993; Charlesworth 1990). This result also indicates that if a population consists of separate subpopulations with different values of ϵ and there is recombination within but not between subpopulations, then the subpopulation with the highest value of ϵ will be favored over the others.

We can also see from Equations 13 and 14 that the position dependence in the distribution of deleterious alleles depends on the sign of ϵ. For synergistic epistasis, recombination will increase the relative frequency of single mutants near the center of the chromosome and decrease the relative frequency of double mutants with large separation between their two mutations, and vice versa for antagonistic epistasis. This is confirmed by simulations, as shown in Figure 4, and makes intuitive sense. When there is synergistic epistasis, bringing together multiple mutations into the same individual decreases mean fitness, while breaking up combinations of deleterious alleles (dispersing them into separate individuals) increases it. The population therefore decreases the effective recombination rate between single mutants by clustering their mutations together in the middle of the chromosome and increases the effective recombination rate for double mutants by clustering their mutations at opposite ends of the chromosome. For antagonistic epistasis, mean fitness decreases when combinations of deleterious alleles are broken up and increases when they are brought together, so we see the opposite effect.

## THE EVOLUTION OF EPISTASIS

When mutations that change the sign or strength of epistasis are possible, epistasis can evolve. A mutation at locus *M* may produce an allele whose effect is to change ϵ, and the subpopulation with this different ϵ may be favored or disfavored by evolution. In the long term, the population will tend to evolve toward those values of ϵ that are most favored. Our analysis is done within the context of this modifier locus approach, where ϵ is controlled by the locus *M*, which gives us a well-defined mutation and recombination structure in which to make concrete calculations. However, the basic idea is more general: we ask whether or not more synergistic or more antagonistic epistasis will tend to be favored, and by how much, regardless of the cause of the change in the epistasis. Naturally the details of the evolution of epistasis in any particular system will depend on whether or not mutations changing ϵ in a given way are possible. Our analysis demonstrates when such mutations will be favored by evolution, given that they have occurred.

We have seen above that, as measured by mean fitness, synergistic epistasis favors recombination. Surprisingly, we find below that the converse is not true. Rather, recombination favors the evolution of antagonistic (or less synergistic) epistasis. Thus, although synergistic epistasis will tend to favor recombination, this recombination will, if possible, eliminate the synergistic epistasis.

#### Asexual populations:

In the absence of recombination the mean fitness in steady state does not depend on ϵ. This means that if two subpopulations with different values of ϵ are placed in competition, neither will be favored over the other in the long term. Depending on the initial mean fitness of the two subpopulations (which may depend, for example, on the genotype of the individual in which the mutation changing ϵ occurs), one or the other may initially become more common. However, they will both eventually reach their mutation–selection steady state, in which they will both have the same mean fitness. This is confirmed with simulations (not shown). Once this happens, the fraction of the population with epistasis ϵ_{n} will drift neutrally.

#### Recombination:

This simple result no longer holds when recombination is possible. With more than one allele possible at the modifier locus, there are two types of recombination events: those between the modifier and the set of all major alleles and recombination between the major alleles. Consider a population in which there are two alleles *M*_{0} and *M*_{1} producing different values of ϵ segregating at the modifier locus. If these two subpopulations are initially at their individual mutation–selection balances, described by and , respectively, they have equal mean fitness and neither is favored. However, recombination brings the frequency distributions *F _{k}* and

*G*of these two subpopulations closer together, as it converts individuals of type

_{k}*F*to type

_{k}*G*and vice versa. This increases the mean fitness of the subpopulation with the smaller (or more negative) value of ϵ, because in this subpopulation the number of individuals with deleterious mutations is reduced. This is illustrated in Figure 2. In contrast, the mean fitness of the subpopulation with the larger value of ϵ is decreased. This means that the subpopulation with more negative ϵ is favored by selection. Thus, recombination favors the development of more antagonistic (or less synergistic) epistasis between deleterious mutations.

_{k}When recombination is rare, we can estimate the strength of the selection for lower ϵ. We imagine a population dominated by one value of epistasis, ϵ_{0}. We introduce, at low frequency, individuals with a different value of epistasis ϵ_{1}, where ϵ_{1} < ϵ_{0}, and calculate the strength of the selection for individuals with the more antagonistic epistasis, ϵ_{1}. We denote by *F _{k}*(

*i*

_{1},…,

*i*) the frequency of individuals with ϵ

_{k}_{0}and

*k*deleterious mutations at major loci

*i*

_{1},…,

*i*and by

_{k}*G*(

_{k}*i*

_{1},…,

*i*) the frequency of individuals with ϵ

_{k}_{1}and

*k*deleterious mutations at

*i*

_{1}, …,

*i*and

_{k}. F*G*denote the total fraction of the population with epistasis ϵ

_{0}and ϵ

_{1}, respectively. When recombination is rare , assuming as before that recombination acts separately from selection, the recursions for

*G*are given to leading order in , , and by(18)(19)(20)where we have defined

_{k}*D*(

_{G}*i*,

*j*) =

*G*

_{0}

*G*

_{2}(

*i*,

*j*) −

*G*

_{1}(

*i*)

*G*

_{1}(

*j*) and in the last recursion we have assumed

*i*<

*j*. The

*F*recursions are the same, with all

*G*'s replaced by

*F*'s and vice versa. Note that we can neglect

*G*and

_{k}*F*for

_{k}*k*≥ 3 because these terms are all higher order in , , and .

When *M*_{1} is rare, we can ignore the terms containing *G _{k}* in the

*F*recursions and the terms proportional to in the

_{k}*G*recursions. This is because when

_{k}*M*

_{1}is rare, recombination events of ϵ

_{1}haplotypes are almost always with ϵ

_{0}haplotypes, and recombination events of ϵ

_{0}haplotypes are also almost always with ϵ

_{0}haplotypes. With this simplification, the recursions for

*F*become simply those given in Equations 10–12 and can be solved to give the steady-state distributions, Equations 13 and 14. The ϵ

_{k}_{0}subpopulation remains in this steady-state distribution as long as

*M*

_{1}is rare. In this steady state, the mean fitness of the ϵ

_{0}individuals is, using Equation 17,(21)

We can now substitute this steady-state *F _{k}* distribution into our recursions for

*G*. This leads to a quasi-steady-state

_{k}*G*distribution, because, as we shall see, the distribution of

_{k}*G*(

_{k}*i.e*., the values of

*G*/

_{k}*G*) quickly reaches its steady state relative to the timescale on which

*G*or

*F*change. For simplicity, we describe this distribution by

*g*

_{1}(

*i*) ≡

*G*

_{1}(

*i*)/

*G*

_{0}and

*g*

_{2}(

*i*,

*j*) ≡

*G*

_{2}(

*i*,

*j*)/

*G*

_{0}. Upon substituting the steady-state

*F*distributions into the recursions for

_{k}*G*, and solving for the quasi-steady-state

_{k}*G*, we find that in this quasi-steady state,(22)where

_{k}*i*<

*j*in the second equation. Note that the effect of

*R*on

*g*

_{1}(

*i*) is higher order in than that of

*r*. We can find the strength of selection for the

*M*

_{1}allele by calculating the mean fitness of the

*G*individuals and comparing it to the mean fitness of the

*F*individuals found in Equation 21,(23)where Δϵ = ϵ

_{1}− ϵ

_{0}. This means that when

*M*

_{1}is rare (

*F*≈ 1), its frequency increases roughly exponentially, at a rate(24)which is the strength of selection for the more antagonistic epistasis. Note that this selection pressure is weak compared to

*s*, justifying our quasi-steady-state assumption.

As the frequency of ϵ_{1} changes, the distribution *g _{k}* also slowly changes. Eventually, the ϵ

_{1}individuals become common and the ϵ

_{0}individuals become rare. In this opposite limit, the distribution of the

*G*is given by the steady-state result, Equations 13 and 14, and we can find a quasi-steady-state distribution for the ϵ

_{k}_{0}individuals. This yields a roughly constant selection pressure for ϵ

_{1}individuals throughout most of the course of the evolution. In Figure 5, we compare these analytical predictions for the strength of selection to simulations.

## DISCUSSION

The patterns of epistasis in natural populations contain information about the way in which genetic and regulatory networks are organized. This epistasis, in turn, has important implications for the future evolution of the population, such as the forces acting on recombination rate and mutational load and driving speciation. The aspects of the biology that determine epistatic interactions are also themselves under evolutionary pressures, both for unrelated pleiotropic effects they have on fitness and because of their impact on the epistasis itself. We have studied this latter effect.

A variety of types of epistasis are important for these questions. Epistasis between beneficial mutations, or the possibility that some mutations may be neutral or deleterious by themselves but beneficial in the right combinations (“sign epistasis”), may be very important to adaptation and speciation (Iwasa *et al*. 2004; Weinreich and Chao 2005). However, combinations of mutations involved in sign epistasis are likely to be rare compared to deleterious mutations. We have focused exclusively on the latter case of the tendency of deleterious mutations to interact either synergistically or antagonistically.

The epistasis that exists among these deleterious mutations is determined by several factors. First, there are presumably only a limited number of ways in which biological networks can change that affect epistasis. It may be that there are fewer ways to change existing genetic networks to make epistasis more antagonistic than to make it more synergistic, or vice versa. Second, mutations that change epistasis may also tend to have direct fitness impacts, which could be stronger than the selection due to their effect on epistasis. These direct fitness effects might “average out,” in the sense that the direct fitness effects of mutations that make epistasis more synergistic are on average the same as those that make epistasis more antagonistic. But this does not have to be the case; it could be that the evolution of epistasis is simply a by-product of these other direct impacts on fitness. Third, mutations that change epistatic interactions are selected because of this change, independent of any direct effects on fitness. This final effect has been the focus of this article.

The patterns of epistasis in natural populations are ultimately an empirical question. Given particular observed patterns of epistasis, an understanding of the direction and strength of one of the evolutionary forces that may have acted to generate this epistasis can help us understand what the observed epistasis might mean. If, for example, we observe an organism where epistasis between deleterious mutations is predominantly synergistic, our analysis has shown that this cannot be because evolution favors such synergistic interactions. Rather, there must be some other reason that we can look for, such as direct advantages to the types of networks that happen to produce synergistic epistasis.

Since experimental evidence on the structure of epistasis in natural populations is at present unclear, our analysis is also useful in showing that there is no theoretical reason to believe that synergistic epistasis is likely to be widespread in nature. Rather, we expect evolution to have favored the development of more antagonistic epistasis—“buffering” between the effects of deleterious mutations. While synergistic epistasis may indeed be present when it happens to be a by-product of other selectively advantageous traits, we can expect that it has been eliminated by evolution where possible.

We have found that the strength of selection against synergistic epistasis and in favor of antagonistic epistasis increases with the recombination rate. Thus in mostly asexual populations, synergistic epistasis is more likely to persist. This poses a problem for the mutational deterministic hypothesis for the advantage of recombination, which posits that sex be advantageous because it increases fitness in the presence of synergistic epistasis between deleterious mutations. While this remains true, recombination will create selection pressure against the synergistic epistasis and thus tend to eliminate the source of its own evolutionary advantage.

The analytical work in this article is strictly valid for , and , that is, for a large number of loci with potential epistatic interactions and subject to rare slightly deleterious mutations. It is only in this parameter range that Equation 24 accurately describes the strength of selection for antagonistic epistasis. However, the more general result that recombination gives a selective advantage to more antagonistic epistasis by breaking up associations between antagonistic epistasis and deleterious mutations should not depend strongly on the values of the parameters in the model. In particular, we expect it to hold for , and . In fact, the effect should be stronger for larger values of *R*, *r*, and *s*.

In a computational study of the evolution of artificial gene networks, Azevedo *et al*. (2006) recently came to the opposite conclusion that evolution favors synergistic epistasis. They found that selection acts to increase the mutational robustness in sexual populations, which tends to be correlated with synergistic epistasis. Their work is not necessarily inconsistent with ours, because in their simulations synergistic epistasis evolves only due to other fitness effects of the changes that create this type of epistasis.

Mutational robustness is closely related to genetic canalization, which refers to the extent to which a genotype is buffered against mutations (Wagner *et al*. 1997). Clearly canalization is also related to epistasis. Synergistic epistasis can be considered to represent canalization, as it means that one or a few mutations have a relatively small effect on fitness, compared to the larger effects of multiple mutations that take a genotype away from the buffered state (Burch and Chao 2004). Since the only selection on canalization is through its effects on fitness through the epistasis that it creates, our analysis thus shows that selection should tend to act against canalization.

However, this conclusion is sensitive to definitions—specifically, to what extent we distinguish the effects of single mutations from their epistatic interactions. In our model, the effects of single mutations are fixed and the epistasis between them is allowed to evolve. If we instead were to allow the effects of single mutations to change, while fixing the fitness of double and triple mutants, we would find that reducing the fitness cost of single mutations would be favored, consistent with the work of Gardner and Kalinka (2006). This would then have the side effect of producing synergistic epistasis and hence canalization. In our view, however, it would be misleading to classify this as the evolution of synergistic epistasis. Rather, it is selection in favor of direct fitness impacts that happen to also be correlated with epistasis. When attention is focused specifically on the selection on epistasis, as in our model, more antagonistic epistasis is instead favored.

Regardless of definitions, we see that buffering of genotypes is favored. When the fitness effects of single mutations are allowed to evolve, evolution favors the buffering of genotypes against their effects—in other words, the reduction of their fitness cost. If we assume that the fitness of double mutants is not correspondingly reduced, this leads to synergistic epistasis and hence canalization or “mutational robustness.” On the other hand, if the fitness effects of single mutations are fixed and their epistatic interactions are allowed to evolve, selection favors an increase in the buffering between the effects of the mutations—that is, more antagonistic epistasis.

We can view our results as a specific example of a more general process by which recombination favors genetic robustness. In a population with multiple segregating modifier alleles affecting the degree of robustness to deleterious mutations, mutations will accumulate in individuals carrying the modifier that confers the most robustness. Recombination will break up the associations between this modifier and deleterious mutations, transferring them to less robust individuals for whom they will cause a greater loss of fitness. Thus, recombination favors the evolution of increased robustness to mutation.

## APPENDIX A: STEADY-STATE FREQUENCIES OF DELETERIOUS ALLELES WITH RECOMBINATION

In this section we outline the iterative approach that leads to an approximate steady-state solution, valid for small *r*, to the recursion equations for the position-dependent frequencies of the deleterious alleles in the presence of recombination, Equations 10–12. Since *D*(*i*, *j*) appears only in the recursions multiplied by *r*, to leading order in *r* we can use the *r* = 0 value, . The equation for then becomes(A1)where . Note that(A2)which gives(A3)Since (see Equation 16 and the surrounding discussion), the in the denominator does not affect the leading term in , and we can drop it, giving Equation 13.

For *F*_{2}(*i*, *j*), we have(A4)From (13), we can see that the contribution of the position dependence of to will be higher order in than the contribution of the *r*| *j* − *i*|*D* term, so we can just use the *r* = 0 value . Plugging this in gives(A5)As before, is higher order in and can be ignored, giving Equation 14.

If desired, we could now iterate this process to find the correction to second order in *r*. We recalculate *D* using Equations 13 and 14, substitute this into the recursion equations, Equations 10–12, and repeat the process outlined above.

## APPENDIX B: SIMULATION METHODS

We tested our analytical results with both deterministic (infinite population) and stochastic (finite population) simulations. In the deterministic simulations, individuals with zero, one, or two mutations were assumed to evolve according to the recursion equations defined in the text, while individuals with more than two mutations were assumed to be nonviable (including higher-order mutants that had a negligible effect on the results.) The simulation results shown in Figures 4 and 5 are from these deterministic simulations. We also checked the results for the average frequencies of deleterious alleles and the strength of selection for antagonistic epistasis using fully stochastic finite-population simulations. For these, we assumed that the population evolved according to a Wright–Fisher model with recombination, mutation, and selection (occurring in that order) and again assumed that individuals with more than two mutations were nonviable. To find the strength of selection, we first ran the simulations until they approached quasi-steady-state distributions and then recorded the difference in mean fitness of the populations with different values of ϵ. The results of these finite-population simulations (not shown) are consistent with our analytical predictions.

## Acknowledgments

We thank Daniel Fisher for many useful ideas and discussions. This work was supported by National Institutes of Health grant GM 28016. D.W. is supported by a National Science Foundation graduate research fellowship.

## Footnotes

Communicating editor: N. Takahata

- Received May 10, 2007.
- Accepted August 8, 2007.

- Copyright © 2007 by the Genetics Society of America