Abstract
An approximate solution for the mean fitness in mutationselection balance with arbitrary order of epistatic interaction is derived. The solution is based on the assumptions of coupling equilibrium and that the interaction effects are multilinear. We find that the effect of morder epistatic interactions (i.e., interactions among groups of m loci) on the load is dependent on the total genomic mutation rate, U, to the mth power. Thus, higherorder gene interactions are potentially important if U is large and the interaction density among loci is not too low. The solution suggests that synergistic epistasis will decrease the mutation load and that variation in epistatic effects will elevate the load. Both of these results, however, are strictly true only if they refer to epistatic interaction strengths measured in the optimal genotype. If gene interactions are measured at mutationselection equilibrium, only synergistic interactions among even numbers of genes will reduce the load. Oddordered synergistic interactions will then elevate the load. There is no systematic relationship between variation in epistasis and load at equilibrium. We argue that empirical estimates of gene interaction must pay attention to the genetic background in which the effects are measured and that it may be advantageous to refer to average interaction intensities as measured in mutationselection equilibrium. We derive a simple criterion for the strength of epistasis that is necessary to overcome the twofold disadvantage of sex.
THEORETICAL population genetics is in many ways the biological equivalent of theoretical physics. However, despite its mathematical and statistical sophistication, population genetics differs from the physical sciences in its lack of a theory for measuring its fundamental parameters. Usually, model parameters such as gene effects, selection coefficients, and mutation rates are introduced without consideration of their domain of application, their scale, or how they are to be operationally measured (Wagner and Laubichler 2000). Any population genetical model must be understood as a representation of some simple subset of the genetic system that is embedded in an unspecified genetic and environmental background. The model parameters are then implicitly defined with respect to this background and can only be operationally measured with reference to it. We have previously argued that epistatic gene interactions need to be defined operationally in terms of their dynamical effects (Wagneret al. 1998; Hansen and Wagner 2001) and we have developed a model of functional epistasis that makes the role of the genetic background explicit (Hansen and Wagner 2001). In this article we use this model to study the effects of epistasis on the mutation load. We show that an explicit consideration of the “reference genotype” can be a powerful conceptual tool that allows us to describe the relationship between epistasis and the mutation load in a way that is not feasible in standard population genetics theory. We also demonstrate that the reference genotype in which gene effects are measured has important consequences for the dynamical and evolutionary interpretations of epistasis.
Although mutation is a fundamental prerequisite for evolvability, most new mutations are deleterious, and every organism carries a load of deleterious mutations (Haldane 1937; Muller 1950). Estimates of the total genomic mutation rate are still controversial but some estimates indicate that the number of new deleterious mutations may average more than one per individual per generation in animals (Crow and Simmons 1983; Crow 1993; EyreWalker and Keightley 1999; Lynchet al. 1999). Under mutationselection balance this translates into a large mutational load, which may have a variety of evolutionary implications. The mutation load may be involved in the evolution and maintenance of recombination and sexual reproduction (Kondrashov 1988), in the evolution of senescence (Rose 1991), in the evolution of reproductive effort late in life (Charlesworth 1990a), in inbreeding depression (Charlesworthet al. 1990), and in the evolution of mate choice (Hansen and Price 1999).
Epistasis is fundamentally implicated in many of the evolutionary consequences of mutation load. If deleterious mutations interact synergistically they may be more efficiently removed by selection and the load is decreased (Kimura and Maruyama 1966; Charlesworth 1990b; Gavrilets and deJong 1993). Kondrashov (1982, 1984, 1988) argued that this effect might be sufficiently strong to outweigh the twofold disadvantage of sex, and help maintain sexual reproduction. Similarly, synergistic interactions among segregating mutations may enhance the differences in the load between individuals. They may, for example, elevate age and sex differences in genetic fitness with obvious consequences for the evolution of mate choice and life histories (Hansen and Price 1999).
THE MULTILINEAR GENOTYPEPHENOTYPE MAP
The multilinear model (Hansen and Wagner 2001) was introduced with two aims. One aim was to provide a dynamically tractable representation of functional epistatic interactions, and the other was to ensure that parameters and variables were operationally defined. Toward the second aim we introduced the concept of a reference genotype. This is a real or abstract genotype that serves as a yardstick in which model parameters can be measured. The variables of the model are “reference effects” of genes. The reference effect of a singlelocus genotype is defined as the phenotypic effect of substituting this genotype into the reference genotype. Let the variable ^{i}y represent the reference effects of genotypes at locus i. These variables can assume any real value. We can describe a given genotype, g, in terms of the reference effects of all the loci at which it is (potentially) different from the reference genotype. If there are n such loci, the genotype is described as the set g = {^{1}y,..., ^{n}y}. A genetic substitution at a locus from ^{i}y to ^{i}y′ can also be assigned a reference effect as ^{i}δ = ^{i}y′  ^{i}y.
To achieve a tractable representation of functional epistasis we made the assumption that a change in the genetic background can at most lead to a linear transformation of the effects of a gene substitution at a locus. In other words, if a gene substitution has an effect ^{i}δ in the reference genotype, then the same gene substitution in another genetic background, g, will be ^{g}^{→}^{i}f ^{i}δ, where the epistasis factor ^{g}^{→}^{i}f, to be defined below, is a constant that depends only on the genetic background g. Thus, all gene substitutions at this locus will be modified by the same factor. Together with the weak, but essential, assumption that the effect of a genotype is independent of the order of substitutions that lead to the genotype, we showed that the map from a genotype g to a univariate trait x can be represented as
The reference effects are defined with respect to wholelocus genotypes only. Substitution of single alleles is a special case of this, and the theory makes no assumptions about dominance. The linearity assumption, however, implies that additive and dominance effects are combined in the same fashion. Hence, there is no distinction among A × A, A × D, and D × D epistasis.
By the “genetic background” of a single or multilocus genotype we mean the state of all other loci in g. The effect of a genetic background, g, on a singlelocus genotype, i, can be described in terms of an epistasis factor
If the reference genotype is changed, say from a genotype r to a genotype r′, the reference effects and epistasis coefficients measured in r′ are
Proofs and further details are given in Hansen and Wagner (2001). A further result that becomes important below, and is proven in appendix a, is that the multilinear model can also be written in terms of epistasis factors as
THE MUTATION LOAD
Let w be fitness with optimal (= maximal) value w_{0}. Further, let the genetic architecture of fitness be described as a multilinear form in the reference effects of n diallelic loci. We start by using the optimal genotype as reference.
Let ^{i}a be the deleterious reference effect of a heterozygote at locus i. Technically this reference effect is then defined on a negative fitness scale (i.e., the reference effect is ^{i}a on the fitness scale). This is convenient as it allows us to represent deleterious effects as positive numbers. Let ^{i}p be the frequency of the deleterious nonrecessive allele at locus i, and let its mutation rate per allele per generation be ^{i}u. The mean deleterious reference effect of the locus is then ^{i}y¯ ≅ 2^{i}a^{i}p. Observe that the actual fitness effect of the heterozygote in a given genetic background, g, is (^{g}^{→}^{i}f ^{i}a). Using this, and assuming all loci are in coupling equilibrium and ignoring back mutations, we have
Again assuming coupling equilibrium, the average fitness is
The genetic load is
Equation 12 does not actually predict the equilibrium mean fitness given a priori mutation rates and other genetics parameters. The reason is that the epistasis coefficients in (12) are defined and measured at mutationselection equilibrium and are thus strictly spoken functions of the model variables rather than constant parameters. However, there is an interpretation of this equation that is nevertheless useful. The first term,
A further result, which is established in appendix c, is that if w(^{1}y′,..., ^{n}y′) is the fitness of a genotype as measured in the equilibrium, then the fitness of the optimal genotype is
EFFECTS OF EPISTASIS ON THE MUTATION LOAD
The first term in (11) and (12) corresponds to the classical result (Haldane 1937) L = U for the mutation load of a trait with an additive genetic basis. If epistasis is present this can be modified in various directions. To see the effects more transparently, we can rewrite (12) as
A novel result of some interest is that higherorder epistasis potentially may have very strong effects on the load if U is large. This is because the effect of mthorder epistasis is proportional to U^{m}. Of course, this interpretation depends critically on the effective epistasis coefficients of higher order not being vanishingly small. There are several reasons why coefficients of higher order may be expected to be smaller. First, note that the effective epistasis coefficients are strict averages of the individual epistasis coefficients only if all loci have equal mutation rates and the number of loci is much larger than the order of the interaction. The second of these caveats is unlikely to be important as fitness components are usually influenced by many loci. However, unequal mutation rates may reduce higherorder effective coefficients well below the average interaction strength. Second, functional epistatic interactions among many loci may be very infrequent, making the effective epistasis coefficients near zero. Nevertheless, provided the interaction density is not too low, we suggest that higherorder interactions may strongly influence the load when U is large.
A surprising result is that the sign of the effects of synergistic epistasis seems to depend on whether the interaction is among an odd or even number of genes. The opposite effects of even and oddorder synergistic interactions on genetic load need careful interpretation. If we start with (11) describing the solution when effects are measured in the optimal genotype, we can observe that there are two effects of epistasis of a given order. One is given directly and the other is given as a modification of lowerorder effects through the epistasis factors. By differentiation of (11) it can be shown that increasing any epistasis coefficient will lead to a reduction in the load provided all the singlelocus epistasis factors are relatively close to one. Thus, at least for weak epistasis, synergistic interactions of any order will reduce the load.
However, as can be seen from (12) and (15), oddorder synergistic epistasis as measured with reference to the equilibrium population always elevates the load. To understand how this is possible we may consider what will happen if we take a population where both second and thirdorder epistasis are synergistic with reference to the optimal genotype and then measure the interaction strengths in the equilibrium population. What we will observe is that secondorder interactions will be much stronger (relative to a population without synergistic thirdorder interactions; synergistic pairwise epistasis will in fact be weaker at equilibrium if there are no higherorder interactions), while thirdorder interactions will be weaker. Thus, the two results are consistent. This underscores the importance of knowing exactly in what sort of genetic background epistatic effects are measured.
Effects of variation in epistasis: It has been argued that variation in epistatic interactions across loci will elevate the load and make the evolution of recombination less favorable (Otto and Feldman 1997; Phillipset al. 2000). As gene interaction is likely to be extremely variable, this is a potentially fatal argument against the hypothesis that recombination is an adaptation to reduce the mutation load. These results were based on the analysis of twolocus models, and it is of interest to see how they hold up in a multilocus analysis.
The effect may be seen most clearly with secondorder epistasis measured with reference to the optimum genotype [see (11)]. Observe that the epistasis coefficients are weighted with the inverse of singlelocus epistasis factors. If there is a lot of variation in epistatic interactions, the epistasis factors will be variable. Some will be larger than one and some will be smaller than one. If a locus has a lot of antagonistic epistatic interactions with other loci, the epistasis factor acting on this locus will tend to become less than unity. This will enhance the effect of all interactions involving this locus. Therefore, the loadelevating effects of antagonistic epistasis will tend to be enhanced by variation in epistatic interactions. Similarly, the loadreducing effects of synergistic epistasis will be diminished, and we predict that variation in pairwise epistasis will increase the load.
However, this argument pertains to epistatic interactions measured with reference to the optimum genotype. Variation in epistasis measured in the equilibrium genotype has no inherent tendency to elevate the load [see (12)]. For this reason, there are heuristic advantages to measuring directional epistasis in the average genotype, as this can be used to assess its effects on the load in an unbiased fashion.
THE AMOUNT OF GENE INTERACTION NECESSARY TO MAINTAIN SEX
Given (15) one may ask how strong epistasis needs to be to compensate for the twofold advantage of asexual reproduction. Let W(N) be the per capita growth rate of a mutationfree sexual population as a function of population density, N. All else being equal we expect the growth rate of a mutationfree asexual clone to be 2W(N). Assume that the deleterious mutations have the same effects on fitness at all densities. In mutationselection equilibrium the mean fitness of a sexual population will then be
This implies that sex is likely to be maintained by synergistic epistasis only if the total genomic mutation rate is large. Given estimates of U and of directional epistasis, (19) provides a test criterion that can be used to evaluate whether synergistic epistasis could be strong enough to maintain sex. To interpret this result it is useful to observe that an epistasis coefficient equal to 2, which would be necessary if U = 1, means that an average deleterious mutation with a 1% effect on fitness would increase the effects of other mutations at least 2% on average. Given that not all mutations affect each other, this must be considered strong epistasis. If U = 10, a mutation with a 1% effect needs to increase the effect of other mutations at least 0.11% on average.
We add the caveat that the above analysis presupposes that both sexuals and asexuals are in mutationselection equilibrium when they start to compete. If an asexual clone arises directly from the sexual population itself, it may be able to invade before it has acquired an elevated load. Thus, complicated dynamics may easily ensue, where asexuals first spread and then get outcompeted as their loads increase.
Equation 18 indicates that sex is most likely to be favored when the per capita growth rate is low. If competition between sexual and asexual individuals is limited to periods of rapid population growth, as may happen in an rselected species, sexual reproduction will be more difficult to maintain.
In conclusion, sex is likely to be favored if the genome is large (i.e., U is large) and the growth rate is low. Assuming that there is no intrinsic difference in the average degree of gene interaction between organisms with small genomes and those with large genomes, sex is expected to be more readily maintained in slowly reproducing organisms with many genes. This is consistent with a number of natural history facts. For instance, the frequency of asexually reproducing species is lower among higher animals with large genomes and slow reproduction. Another pattern is the prevalence of parthenogenetic species in pioneer species that are characterized by high reproductive rates. We note, however, that there are other hypotheses that can account for these differences (Maynard Smith 1976; Bell 1982).
THE FIXATION LOAD
Another consequence of synergistic epistasis is its potential ability to halt Muller’s ratchet in asexuals (Kondrashov 1994). However, this may be unlikely as it depends critically on the unrealistic assumption that the effects of all mutations are the same (Butcher 1995). Here, we suggest another potentially important role for epistasis in the fixation of deleterious mutations.
We suggest that synergistic epistatic interactions tend to become reduced and antagonistic interactions elevated as we move away from the optimum genotype. This comes about because the effect of a new deleterious mutation in a background g is ^{g}^{→}^{i}f ^{i}δ. Due to the fact that the epistasis coefficients that determine the epistasis factor are symmetrical across loci (^{ij}ε = ^{ji}ε, and so on), this means that a mutation on a locus with a lot of antagonistic interactions with other loci will tend to have reduced effect due to segregation or fixation of mutations on other loci. Therefore, its fixation probability will be increased. For the same reason, a mutation with a lot of synergistic interactions will have reduced fixation probability. Due to the fact that the fixation probability of deleterious mutations is an extremely nonlinear function of the effect of these mutations (Bürger and Ewens 1995), this may lead to extreme differences in the fixation probabilities of mutations with different sorts of epistatic interactions. The ratchet clicks faster with antagonistic epistasis.
The implication of this is that the process of fixing deleterious mutations may be associated with the evolution of antagonistic epistatic interactions, and this will be true for epistasis of all orders. If there are suites of mutations with mutually antagonistic interaction (e.g., compensatory mutations), then the fixation, or even segregation, of some of these will alter the genetic background in a way that may greatly elevate the fixation probability of the other mutations. Thus, in small populations, we may expect a genetic architecture where many alleles have strong deleterious effects with reference to an optimal or ancestral genotype, but where antagonistic epistasis keeps the fixation load bounded. Note also that, in such a population, the individual epistatic interactions as measured in the optimal genotype do not have to be strong. Thus, it is essential to measure epistatic interactions with reference to the population in which they occur. The effect of such a process on the evolution of recombination needs investigation.
DISCUSSION AND CONCLUSION
Epistasis is an evolving entity. The phenotypic effects of gene interactions depend on the genetic background in which the interactions take place and will change when the genetic background is changing. This has also been found in empirical studies (Moreno 1994; Polaczyket al. 1998). Therefore, both empirical and theoretical studies should pay close attention to the genetic background in which epistasis is measured. The concept of a reference genotype makes this explicit, and has led to several new insights about the effects of gene interaction on the mutation load.
What do the results imply for the maintenance of sex and recombination? The most important implication is that the pattern of higherorder epistasis may be important. The deterministicmutation hypothesis (Kondrashov 1988) depends on a large value of the total deleterious mutation rate, U, and this is precisely the situation in which higherorder interactions start to affect the load. How large U needs to be before higherorder interactions become important depends critically on the relative strength of higherorder epistasis. Although difficult, the estimation of higherorder epistasis coefficients may be helpful to understand the role of epistasis in the maintenance of sex and recombination.
The notion that synergistic epistasis is reducing the load in a sexually reproducing population has also been refined. We have confirmed the insight that the mutation load is reduced by synergistic epistasis, at least as long as the epistatic effects are not very strong. However, this is strictly true only if the statement refers to gene interactions measured in the optimal genotype. As measured in the equilibrium population, synergistic interactions among oddnumbered genes will in fact always elevate the load. Note, though, that this may then be compensated by a relative increase in evenordered interaction strengths.
The asymmetry between odd and evenordered epistasis is indeed a puzzling result. One interpretation (suggested by J. Hermisson, personal communication) is that it may be due to an asymmetric effect of even and oddordered epistasis on beneficial alleles. At equilibrium, a synergistic interaction among an even number of deleterious alleles implies an antagonistic interaction among alleles that are beneficial relative to the mean. However, a synergistic interaction among an odd number of deleterious alleles implies a synergistic interaction among beneficial alleles. It is possible that this asymmetry can explain the different effects on the load of even and oddorder interactions as measured in the “average genotype.”
Variation in pairwise epistatic interactions may indeed be said to elevate the load if the statement refers to pairwise epistasis measured with reference to the optimal genotype (as in Otto and Feldman 1997; Phillipset al. 2000). However, we found no obvious relationship between the size of the load and variation in epistasis as measured in the equilibrium population.
Epistatic interactions have been measured in a number of studies including quantitative trait loci data (e.g., Cheverud 2001) or in genotypes in which mutations have been induced or been allowed to accumulate Mukai 1969; Clark and Wang 1997; Whitlock and Bourguet 2000 in Drosophila; De Visser et al. 1996, 1997a in Chlamydomonas; De Visseret al. 1997b in Aspergillus; Elena and Lenski 1997 in Escherichia coli). Before such estimates can be related to the mutation load and used to test evolutionary hypotheses, the genetic background in which they are measured should be precisely characterized. If the genetic background can be assumed to have reached mutationselection equilibrium, then the effects of induced mutations can be interpreted with our equilibrium reference theory, and their impact on the load and the advantage of recombination can be evaluated by (15) and (19).
Unfortunately, with the possible exception of Elena and Lenski’s (1997) study of asexual prokaryotes, no existing studies can be said to measure gene interactions with reference to an equilibrium background or to an optimal background. There is a need for measures of epistasis in wildtype backgrounds, or at least in backgrounds where the fitness effects of the genetic differences to a wildtype or optimal genotype are known, so that the effects of a change in reference can be assessed.
The most useful measures of epistasis may very well be those that are made in natural equilibrium populations or at least in lab populations that can be assumed to have reached mutationselection balance in a stable environment. This is because the effects on the load are then not obscured by variation in epistatic effects, and because this is the situation in which evolution actually takes place. It is therefore encouraging that largescale studies of epistasis in natural populations are now emerging (e.g., Fenster and Galloway 2000).
APPENDIX A
Result: The genotypic value for a trait under multilinear epistasis in genotype g can be written as
Proof: We need to show that the formula is equivalent to the multilinear model,
APPENDIX B
If there are multiple alleles segregating at a locus we can derive (8) in the following way. Let ^{i}q be the frequency of the optimal allele at locus i, and let ^{i}ā be the average deleterious reference effect of the locus when a deleterious allele is present. The mean reference effect at the locus is then ^{i}y¯ = 2^{i}ā(1  ^{i}q). Let ^{i}x be an indicator variable that takes the value 1 if a random allele at locus i is the optimal allele and 0 otherwise. The expectation of ^{i}x is then ^{i}q. The Price equation then gives
The local stability of the equilibrium reference effects can be investigated through the Jacobian matrix for the system of equations
Stability requires that all eigenvalues of the Jacobian are inside the unit circle in the complex plane. Note first that, as long as the epistasis factors are positive, the diagonal elements are always less than one. Gerschgorin’s theorem says that all eigenvalues must be inside the union of the disks centered at the diagonal elements with radii equal to the sum of the absolute values of the offdiagonal elements in the corresponding row. Thus, a sufficient condition for local stability is that the following relations hold for all loci, i,
APPENDIX C
Result: If w(^{1}y′,..., ^{n}y′) is the fitness of a genotype as measured in mutationselection equilibrium, then the fitness of the optimal genotype is
Proof: When we use the mutationselection equilibrium genotype as reference genotype, the mean reference effects are 0 at all loci. Thus the mean fitness is w(0,..., 0). To find the fitness of the optimal genotype, all we need to do is to find the reference effects of the optimal genotypes at each locus and substitute into w(^{1}y′,..., ^{n}y′).
Obviously, the reference effects of the singlelocus optimal genotypes are 0 when measured in the optimal genotype. We use this as a starting point, and apply the change of reference formula (4) to compute the reference effects measured at equilibrium.
The variable ^{i}Δ is the effect of the new reference as measured in the old reference. Thus, in (C2), ^{i}Δ is the mean effect at the locus as measured in optimal genotype. This is given in (8) as
QED
Acknowledgments
The authors are grateful to Reinhard Bürger, David Houle, Alex Kondrashov, and an anonymous reviewer for helpful comments on previous versions of the manuscript; and to Joachim Hermisson for many stimulating discussions and helpful insights that improved the manuscript. A postdoctoral fellowship from the Norwegian science council supported T.F.H.
Footnotes

Communicating editor: M. K. Uyenoyama
 Received September 9, 2000.
 Accepted February 15, 2001.
 Copyright © 2001 by the Genetics Society of America