The Anomalous Effects of Biased Mutation

Corresponding author: D. Waxman, University of Sussex, Brighton BN1 9QG, Sussex, United Kingdom., d.waxman{at}sussex.ac.uk (E-mail)

Communicating editor: M. W. FELDMAN

	ABSTRACT

TOP ABSTRACT MODEL RESULTS DISCUSSION APPENDIX LITERATURE CITED

A model is presented in which alleles at a number of loci combine to influence the value of a quantitative trait that is subject to stabilizing selection. Mutations can occur to alleles at the loci under consideration. Some of these mutations will tend to increase the value of the trait, while others will tend to decrease it. In contrast to most previous models, we allow the mean effect of mutations to be nonzero. This means that, on average, mutations can have a bias, such that they tend to either increase or decrease the value of the trait. We find, unsurprisingly, that biased mutation moves the equilibrium mean value of the quantitative trait in the direction of the bias. What is more surprising is the behavior of the deviation of the equilibrium mean value of the trait from its optimal value. This has a nonmonotonic dependence on the degree of bias, so that increasing the degree of bias can actually bring the mean phenotype closer to the optimal phenotype. Furthermore, there is a definite maximum to the extent to which biased mutation can cause a difference between the mean phenotype and the optimum. For plausible parameter values, this maximum-possible difference is small. Typically, quantitative-genetics models assume an unconstrained model of mutation, where the expected difference in effect between a parental allele and a mutant allele is independent of the current state of the parental allele. Our results show that models of this sort can easily lead to biologically implausible consequences when mutations are biased. In particular, unconstrained mutation typically leads to a continual increase or decrease in the mean allelic effects at all trait-controlling loci. Thus at each of these loci, the mean allelic effect eventually becomes extreme. This suggests that some of the models of mutation most commonly used in quantitative genetics should be modified so as to introduce genetic constraints.

MANY mutations affect continuously distributed traits such as height and weight (LYNCH and WALSH 1998 ) that are under stabilizing selection. Theoreticians studying the evolution of continuously distributed traits have generally assumed that, while some mutations tend to increase the value of a trait, others decrease it, and as a consequence, the average mutational effect is zero (BULMER 1980 , BULMER 1989 ; TURELLI 1984 ). Thus most theoretical treatments to date have assumed that mutations do not tend to cause any directional change in the mean value of a trait. In this sense, mutation has been assumed to be unbiased. There is, however, no a priori reason to assume unbiased mutation (except, perhaps, for mathematical convenience). Furthermore, experimental data suggest that, contrary to the hypothesis of unbiased mutation, mutations do often affect the mean value of phenotypic traits (SANTIAGO et al. 1992 ; LYMAN et al. 1996 ; MACKAY 1996 ; KEIGHTLEY and OHNISHI 1998 ). There is thus a compelling rationale to investigate the impact of biased mutation.

In this study, we have adopted a standard model of stabilizing selection, where an optimal value of the trait exists. We find that the dependence of the population's mean phenotypic value on the degree of mutational bias is nonmonotonic. As such, under some conditions, increasing the extent of mutational bias can actually lead to a reduction in the deviation of the population's mean phenotypic value from its optimal value.

We use a modified version of the model of mutation that was originally introduced by CROW and KIMURA 1964 and employed in a large number of important articles (e.g., LANDE 1975 ; TURELLI 1984 ). Biased mutation of arbitrary degree was included in the model of CROW and KIMURA 1964 and investigated by KIMURA 1965 . Here we show that the analysis presented by KIMURA 1965 is mathematically inconsistent. In particular, we show that if this model incorporates mutational bias then the mean allelic effects do not equilibrate (as Kimura assumed). Instead, they tend to increase indefinitely in absolute value (the mean phenotypic value does, however, equilibrate). Thus over sufficient time, allelic effects can acquire extremely large magnitudes. Furthermore, this tends to happen even if only one locus, out of all of the loci affecting the trait, experiences biased mutation, and the degree of bias of this locus is very slight. It is thus an implication of the results presented here that one of the most commonly used mathematical models of mutation is biologically implausible and requires some modification. To overcome this biological implausibility, we introduce a modification, based on previous work by ZENG and COCKERHAM 1993 , and show that this does yield biologically plausible results, with the population approaching equilibrium at long times.

Consideration of biased mutation is common when the trait under consideration is fitness, as most researchers believe that the vast majority of fitness-altering mutations cause a decline in fitness (CROW 1979 ; KONDRASHOV 1988 ; KEIGHTLEY 1996 ). In addition, a few previous studies have also considered biased mutation in the context of a continuously varying phenotypic trait (other than fitness) that is under directional selection (KIMURA 1965 ; IWASA and POMIANKOWSKI 1991 ; POMIANKOWSKI et al. 1991 ). In such models, mutational bias is expected under directional selection, when favorable alleles achieve high frequencies. Here, however, we consider the effects of bias under symmetrical stabilizing selection, and the study that is closest to the work presented here was published by KIMURA 1965 . Kimura's analysis was restricted to a range of parameter values that are now considered to be biologically implausible, at least for outcrossing species. The analysis presented here focuses upon the parameter range that is currently believed to be closest to biological reality for outcrossing species (LYNCH and WALSH 1998 ).

	MODEL

TOP ABSTRACT MODEL RESULTS DISCUSSION APPENDIX LITERATURE CITED

Consider a randomly mating population of dioecious sexual organisms, with no sexual dimorphism. The population size is assumed to be sufficiently large such that stochastic effects (genetic drift) can be ignored. Individuals are subject to selection on the value of a single phenotypic trait. The phenotype of a particular individual is assumed to depend on the individual's "genotypic value," G, plus a normally distributed environmental noise component, . Using z to represent an individual's phenotypic value, we have z = G + . The distribution of is assumed to be independent of G and has a mean of zero and a standard deviation of V_e. Following convention and without loss of generality, we scale all variables so that V_e is set to unity.

Individuals are diploid with n freely recombining loci that additively affect the genotypic value. These loci are labeled 1, 2, ... , n. The DNA sequence of an allele determines its effect on genotypic value, and the effects of the maternally and paternally inherited alleles at locus i are denoted by x_i and y_i, respectively.

Following CROW and KIMURA 1964 and many subsequent authors, we take the allelic effects to be continuous and to have an infinite range: > x_i, y_i > -. Additivity in the determination of the genotypic value leads to . Apart, possibly, from the initial generation, maternally and paternally inherited alleles have identical distributions. Because of this, we need to refer only to the distribution of alleles of maternal origin.

Generations are discrete, with all parents dying soon after the birth of offspring. Some offspring die before reaching reproductive maturity due to stabilizing viability selection. We confine ourselves to parameter ranges for which selection is weak at the level of the trait (see below). We can therefore employ a quadratic function to describe stabilizing selection, following the example of many authors including KIMURA 1965 and BULMER 1989 . Thus an individual with phenotypic value z has a probability of surviving viability selection that is proportional to 1 - s*(z - z_opt)² (where s* 0). The value of s* is a measure of the strength of stabilizing selection on phenotypes and individuals of optimal phenotype have z = z_opt. Note that the probability of surviving should be set to zero for values of z yielding 1 - s*(z - z_opt)² < 0; however, for the parameter values considered in this work, the probability with which this occurs is of order 10^-9 and thus negligible for practical purposes. Thus, the simple quadratic viability function is taken to hold without restriction on z.

We focus, in this study, on the impact of relatively weak selection: s* << 1. Under weak selection, a quadratic selection function gives results that are very close to those produced by a Gaussian selection function; however, a quadratic function is, mathematically, more amenable to analysis.

The preceding assumptions allow us to derive the effect of selection on the distribution of genotypic values. By averaging over environmental effects, it can be shown that the probability of survival for an individual with genotypic value G is proportional to

(1)

Here s = s*/(1 - s*) is a measure of the strength of selection on genotypic values.

Gamete formation involves standard Mendelian segregation and free recombination. The population members that have survived viability selection—termed adults—undergo random mating and proceed to produce new offspring.

Each of an individual's 2n alleles is a copy of an allele present in one or the other of the individual's parents. The effect of an allele in an offspring is identical to that of the parental allele, of which it is a copy, unless a mutation occurred during its production. The per-allele rate of mutation at locus i is denoted by µ_i, where 1 µ_i 0. The expected number of new mutations that affect the trait, per individual, per generation, U, is given by .

Let us now specify the mutation function. We have chosen a relatively general formulation that encompasses a number of previous approaches. In particular, we make the usual assumptions that mutations to different alleles occur independently and that mutant allelic effects are continuously distributed (CROW and KIMURA 1964 ; KIMURA 1965 ) and have a Gaussian form (LANDE 1975 ; BULMER 1980 ; TURELLI 1984 ).

Let x represent the effect of a particular allele at locus i in a particular offspring. Let x* represent the effect of the parental allele from which the offspring allele was copied. If no mutation of the allele occurred in the production of the offspring then x = x*. If a mutation did occur, then the value of x is chosen from a normal distribution with variance m_i² and mean x* + b_i, where 1 0. In other words, at locus i, the probability density function for the allelic effects of new mutations is given by

(2)

Let us consider the implications of this formula. If = 1 and b_i = 0 for all values of i, then we have the model of mutation most commonly used in quantitative genetics (e.g., LANDE 1975 ; TURELLI 1984 ). In this model the mean mutant allelic effect is simply the effect of the parental allele, x*.

Next, consider the case where = 1 and b_i 0. In this case the distribution of mutants has a mean of x* + b_i and we can interpret b_i as the bias introduced by mutation: on average, mutations at locus i change allelic effects by an amount b_i. We note that while CROW and KIMURA 1964 and KIMURA 1965 did not directly concern themselves with moments of the mutation distribution higher than two, they did allow the distribution to have a nonzero mean and thus implicitly incorporated mutational bias into their calculations.

Consider now the case < 1 and b_i = 0. This model of mutation was initially formulated by ZENG and COCKERHAM 1993 and the distribution of mutants has a mean of x*. The model incorporates the idea of genetic constraints, so that very extreme allelic effects are unlikely to arise as a cumulative consequence of mutation. Thus, even when selection is absent (s = 0), alleles with very extreme effects will not become common in the population as a result of mutation. Instead, allelic effects will remain clustered around zero. It does seem reasonable to incorporate some sort of genetic constraint. Otherwise one gets biologically implausible implications such as large populations yielding extremely large amounts of phenotypic variation on traits that are not under selection. Note that a special example of < 1 is the case = 0 and this corresponds to the house-of-cards model of mutation (KINGMAN 1978 ).

Finally, let us consider the case where < 1 and b_i 0. This model is a combination of a Gaussian mutation model and ZENG and COCKERHAM's (1993) regression model. In this case the distribution of mutants has a mean of x* + b_i. Thus b_i can be interpreted as the mean deviation of mutant alleles from x*. The model allows us to consider mutational bias in situations where very extreme allelic effects are unlikely to arise. Thus, of the models discussed here, this combined mutation model is the most realistic.

In what follows, all summary statistics that describe the population (phenotypic values, genetic variance, etc.) are measured immediately after the birth of the offspring and before any selection has taken place.

	RESULTS

TOP ABSTRACT MODEL RESULTS DISCUSSION APPENDIX LITERATURE CITED

It has been possible to produce analytical approximations of the model in a number of relevant cases and these have been supplemented with numerical studies. The appropriate analytical approximations depend on the combination of parameter values. The standard reference on this subject (LYNCH and WALSH 1998 ) suggests plausible ranges of the relevant parameters. In particular, the following is thought to be likely for many (or most) quantitative traits:

1. The expected number of new mutations per generation per individual that affect the trait, , satisfies U << 0.05.

2. The variance of mutant effects at any locus, m_i², satisfies m_i² << 1 (LANDE 1975 ; TURELLI 1984 ).

3. The strength of selection on genotypic values, s, satisfies s << 1 (GARCIA-DORADO and MARIN 1998 ).

4. The strength of selection acting on allelic effects at any locus, i, is much stronger than the effects of mutation at the locus, such that µ_i/(sm_i²) << 1 (TURELLI 1984 ).

We assume that conditions i–iv hold. In addition, we assume that bias is not large compared with the strength of selection, in the sense sb² < 0.05. As we shall see, all behavior of interest occurs when b is substantially smaller than the requirements of this inequality. Therefore this assumption does not place any important limitation on the scope of this work.

The results presented below apply when the preceding assumptions are met. See Table 1 for notation.

Results for equivalent loci with = 1:
One case where considerable analytic progress is possible is where = 1 and all loci have identical parameter values governing mutation. This is the case where the parameters µ_i, m_i, and b_i do not have any variation in value across loci. We refer to these universal values as µ, m, and b, respectively. It would be surprising to find a case of exactly equivalent loci in nature; however, as we shall see, results for the case of equivalent loci are helpful in predicting the outcome of other, more realistic, cases. In addition, in one case (albeit a degenerate one) a lack of variation in parameter values among loci automatically arises. This is where only a single locus is under selection (n = 1). In the next section we consider equivalent loci with < 1 and again the value of the analysis is for the insight gained for more realistic cases.

For equivalent loci, with = 1, we can produce estimates of summary statistics using the analysis presented in Equivalent loci with = 1 in the Appendix In particular, we estimate the equilibrium variance in genotypic values, G, among offspring, V_G, and also estimate the equilibrium genetic load, , where w(G) is given in Equation 1, is its mean equilibrium value, and L is proportional to the fraction of the population that fails to survive viability selection. The estimates given in the Appendix for V_G and L are, to first order in µ, unaffected by the degree of bias, b, and thus, to this order, identical to well-known approximations that have appeared in the literature (LYNCH and WALSH 1998 ): V_G 2nµ/s and .

We turn now to the behavior of the mean phenotypic value. Our estimate of the equilibrium mean phenotypic value among offspring, denoted , is given by

(3)

where

(4)

and corrections to in Equation 3 are O(µ²). Note that D(ß) can be written in terms of special functions; however, we have found the form given in Equation 4 to be most useful since all of its most important properties are readily derivable from this expression.

To obtain an estimate of the error of the expressions given above for V_G, L, and , we have compared the analytical results, given above, with highly accurate numerical results. The latter followed from numerical solution of the equation of the Appendix governing the distribution of allelic effects, Equation A1, using the method of WAXMAN 2003 . We note that under the approximations made in the Appendix, V_G and L simply follow from a sum over one-locus quantities, while is independent of the number of loci, n. As a consequence the estimates of errors are effectively on one-locus quantities and hence independent of n. With the plausible parameter values µ = 10^-5, m = 0.2, and s = 0.025 (TURELLI 1984 ) and biases over the range 4 b/m 0, we have found that for V_G, L, and the difference between numerical and analytical results is <3%.

Note that an implicit assumption underlies the results for the summary statistics presented above, for equivalent loci. This is that the distribution of allelic effects in gametes (a function of the allelic effects at the n loci that affect the trait) eventually comes to equilibrium. We have used numerical methods to test this assumption and the results are in accord with those from previous studies (PHILLIPS 1996 ). In particular, for all the cases of equivalent loci that we have numerically examined, the distribution of allelic effects in gametes does indeed come to equilibrium. However, the "position" of this equilibrium—which is specified by the mean effects of alleles at all n loci—depends on the initial distribution of allelic effects.

Let _i be the mean equilibrium effect of alleles of maternal origin, at locus i, in newborn offspring. It is identical to the corresponding quantity of paternal origin. With the mean equilibrium phenotypic value, as given by Equation 3, it follows that any set of n mean allelic effects satisfying is a possible end point of the dynamics (an equilibrium). In this sense, the equilibrium is "neutral." However, analysis in the Appendix shows that given that an equilibrium is achieved, and given mutationally equivalent loci, the shapes of the distributions of allelic effects at all loci (i.e., the marginal distributions of allelic effects) are identical. Thus, while the second and all higher central moments of the allelic-effect distributions are the same for all loci, the means of these distributions are generally different.

Let us now consider the effects of bias, ignoring the case b = 0, where mutation is unbiased and . From Equation 3 it follows that a finite positive b yields a mean phenotypic value that is larger than the optimal phenotypic value, i.e., > z_opt. When b is increased from zero, initially - z_opt increases approximately linearly with b. However, as b becomes larger the rate of increase in - z_opt declines until b reaches a critical value, and any further increase in b produces a decrease in - z_opt. For sufficiently large values of b, the value of - z_opt is, to a good approximation, proportional to b^-1. An example of this sort of nonmonotonic behavior is given in Fig 1.

View larger version (16K): In this window In a new window Download PPT slide   Figure 1. A plot, for mutationally equivalent loci, of the deviation of the equilibrium mean phenotypic value, , from its optimal value, zopt, as a function of mutational bias, b. The latter is measured in units of the standard deviation of mutant effects, m. The solid curve describes the case with regression parameter = 1 and is based on Equation 3. The dashed curve is for = 0.2, i.e., covered by the case < 1, and is based on Equation 6. The parameter values common to both curves are µ = 10-5, m = 0.2, s = 0.025, n = 10, and zopt = 3.

Results for b < 0 closely parallel those for b > 0 and can be determined from the latter from the relation , which follows directly from the property D(-ß) = -D(ß). This is apparent in Fig 1.

Note that the maximum value of the function D(ß), of Equation 4, is given by D_max 0.77, a value that is independent of all parameters. The maximum of D(ß) occurs when ßD(ß) = 1, i.e., when ß 1.31. The existence of a maximum of D(ß) implies a definite limit to the degree of deviation that biased mutation can cause in the equilibrium mean phenotype among offspring, , from the optimal phenotype, z_opt. For the parameter ranges specified previously, the absolute value of the maximum deviation, in terms of phenotypic standard deviations, is given by

(5)

The above equation indicates that the maximum possible absolute deviation of from the optimum, not surprisingly, becomes larger when µ is increased or s is decreased, since both of these changes enhance the role of mutation relative to selection. The maximum deviation decreases as n increases because the equilibrium phenotypic variation increases with the number of loci controlling the trait, while | - z_opt|_max is found to be independent of n. Thus increasing n decreases the effect of bias, when measured in phenotypic standard deviations. More intriguing is the effect of the standard deviation of mutant effects, m. The maximum possible deviation of from z_opt is highest when m is small. Furthermore, since D(b/m) reaches a maximum when b 1.3m, we find that the maximum deviation from the optimum occurs when both the degree of bias and the standard deviation of mutant effects are small. This intuitively makes sense since when m and b are small, many mutant offspring born to parents with nearly optimal genotypes will also have nearly optimal genotypes. This allows for the survival of a large proportion of the mutants and for the accumulation of biased mutations over the course of many generations.

Results for equivalent loci when < 1:
When all the selected loci are equivalent, but < 1, the preceding results are modified and analysis covering this case is contained in Equivalent loci with < 1 in the Appendix In particular, numerical investigation indicates that the distribution of genotypes no longer simply comes to a neutral equilibrium. Instead, the distribution comes, after some time, to a unique and stable equilibrium where mean allelic effects, _i, at all loci are uniquely determined. The equilibrium distribution is thus independent of the initial distribution and equivalence of loci results in the _i having identical values for all loci: . A unique equilibrium arises since a regression parameter, < 1, corresponds to an additional evolutionary force in the system that directly couples to the allelic effects. It destroys the property of the case with = 1 that, at equilibrium, any sets of mean allelic effects, _i, that lead to the equilibrium value of are equally good candidates for an equilibrium.

As far as summary statistics are concerned, we find that under the same analytical approximations used for the case = 1, and hence the same accuracy, the genetic variance, V_G, and genetic load, L, are, to leading order in the allelic mutation rate µ, unaffected by the degree of bias, b, and again given by V_G 2nµ/s and . Furthermore, the mean phenotypic value, , is now given by

(6)

where D(ß) is given in Equation 4. Thus, there is still a nonmonotonic dependence of upon the value of bias, b. However, when < 1, it is no longer the case that . This is illustrated by the dashed curve in Fig 1.

Results for nonequivalent loci for = 1:
We have, so far, confined ourselves to consideration of cases where all the mutational parameters (µ_i, m_i, and b_i) are the same at all loci. Let us now relax this assumption by allowing variation among loci in the mutational parameters for the case = 1. Analysis covering this case is contained in Nonequivalent loci with = 1 in the Appendix In particular, it is shown in the Appendix that in this more general case, we have an intriguing result: the assumption that allelic distributions at all loci come to an equilibrium generally leads to a mathematically inconsistent set of equations. Thus, in general, the population cannot equilibrate when loci are nonequivalent and = 1.

Given this lack of attainment of equilibrium, what is the long-time behavior? Numerical investigation for n > 1 mutationally nonequivalent loci shows the reason for the inconsistency mentioned above. Nonequivalent loci typically lead to a situation where the allelic distributions at every locus continue to change indefinitely. More specifically, the typical long-term behavior is that at every locus there is a roughly linear change in the value of the mean allelic value, _i, with time. This appears to generally occur at a rate smaller than the mutation rate. This is caused by a continual turnover of alleles at every locus, such that common alleles become rare and new mutations multiply and become more common (see Fig 2).

View larger version (17K): In this window In a new window Download PPT slide   Figure 2. A plot of the mean allelic effects, i (solid curves), as a function of time, in generations. The case illustrated is for four nonequivalent loci with = 1. The dashed curve shows the behavior of , where is the mean phenotypic value. While the mean allelic effects exhibit continuous change over the timescale of 106 generations, the mean phenotypic value settles down to its asymptotic value on a much shorter timescale (on the order of 103 generations). This figure was produced using a continuous-time approximation to the dynamical equations, which is highly accurate for the parameter values given below. For all loci, allelic effects were discretized, with a splitting (separation of adjacent allelic effects) of 0.04, and other parameter values adopted were s = 0.025, zopt = 0.267, [u1, u2, u3, u4] = [3.9, 2.4, 0.9, 1.1] x 10-5, [m1, m2, m3, m4] = [0.15, 0.24, 0.11, 0.13], and [b1, b2, b3, b4] = [0.13, 0.02, 0.15, 0.10]. The initial distributions of allelic effect were taken to be independent Gaussians, with different means and different variances.

Despite the turnover in alleles, when the inequalities given above in the first paragraph of RESULTS apply, the numerical studies indicate that the genetic variance and genetic load are reasonably close to the results that apply in the absence of bias: and . As an example, for the parameter values used in Fig 2 (see Fig 2 legend), the differences between the analytical predictions for V_G and L and the numerical results are <4%.

The result given in Equation 3 for the equilibrium mean phenotypic effect, , for equivalent loci with = 1 may be used to estimate the corresponding quantity when loci are not equivalent. The most straightforward procedure is to evaluate Equation 3 at the mean values of the mutational parameters µ_i, m_i, and b_i. As an illustration of this, we have compared the long-time value of of Fig 2 (which is the outcome of numerical solution) with Equation 3, evaluated at the mean values of the mutational parameters used in Fig 2 (see Fig 2 legend). The difference between the two values of is found to be <1%.

Results for nonequivalent loci when < 1:
What are the consequences of genetic constraint ( < 1) when loci are nonequivalent? This question is very difficult to fully address although some progress is made in Nonequivalent loci with < 1 in the Appendix In particular, extensive numerical study strongly suggests that when < 1, the population comes to a unique and stable equilibrium. Thus when < 1 we do not see the perpetual turnover of alleles that occurs when = 1 since mutational regression is evidently sufficient to stop turnover. We note that under the same approximations used in previous cases possessing an equilibrium, it is possible to analytically determine that the genetic variance and genetic load are, to O(µ), unaffected by mutational bias, so and . Furthermore, at each locus a distribution of allelic effects, characteristic of that locus, always becomes established, and this occurs regardless of the initial genotypic distribution. For the mean phenotypic value, we are able to determine the approximate bound

(7)

which is independent of the value of and again indicates the highly limited extent to which mutational bias can affect the mean phenotypic value.

Let denote the average of b_i over all loci, weighted by the mutation rate

(8)

Then the dependence of on is qualitatively similar to the nonmonotonic dependence of on b that was found for equivalent loci with < 1; see Equation 6. As an example, we have considered n = 4 loci, with a substantial level of constraint, namely = 0.2 and with mutation rates, µ_i, mutational standard deviations, m_i, and an optimal phenotypic value, z_opt, that are identical to those used in Fig 2 (see Fig 2 legend). We have produced a number of sets of randomly generated mutational biases, [b₁, b₂, b₃, b₄] and plotted the long-term values of , in Fig 3, against the weighted average of bias, . For comparison, we have also plotted an estimate of motivated by Equation 6: namely z_opt + (µ/(sm))D(/m - (1 - )z_opt/(2nm)), where the values of µ and m are taken as arithmetic mean values of the corresponding mutational parameters across loci. As is evident from Fig 3, the form for given for equivalent loci with < 1 (Equation 6) provides a useful estimate of the results for loci that are genetically constrained and mutationally nonequivalent.

View larger version (21K): In this window In a new window Download PPT slide   Figure 3. A plot, for mutationally nonequivalent loci, of the deviation of the equilibrium mean phenotypic value, , from its optimal value, zopt, as a function of weighted mutational bias, (Equation 8). The latter is measured in units of the mean standard deviation of mutant effects, m, across loci. The case illustrated is for four nonequivalent loci with = 0.2 with mutation rates, µi, mutational standard deviations, mi, and an optimal phenotypic value, zopt, which are identical to those used in Fig 2 (see Fig 2 legend). The vertical bars give the values of - zopt for 31 different sets of randomly selected mutational biases of common sign (each set contains the values of b1, b2, b3, and b4). For each set of mutational biases we calculate the value of /m, and this determines the height of each vertical bar. For comparison, the dashed curve is an estimate of - zopt motivated by the result for equivalent loci, Equation 6. In this equation, the values of µ and m are taken as arithmetic means, across loci, of the corresponding mutational parameters.

	DISCUSSION

TOP ABSTRACT MODEL RESULTS DISCUSSION APPENDIX LITERATURE CITED

In this study we have considered the evolutionary consequences of a biased-mutation process that affects a trait that is undergoing stabilizing selection. Stabilizing selection tends to bring the mean phenotypic value closer to the optimal phenotypic value (z_opt). However, when mutation is biased, the new mutations arising during every generation tend to change the value of the trait. The overall degree of bias in mutations (taking all pertinent loci into account) can be characterized by , which is given in Equation 8 and is the weighted average of the degree of bias at each locus, with the weighting determined by the mutation rate at each locus. If ||, the absolute value of , is large, then, on average, mutations tend to cause substantial directional changes in trait values. If || is small, then the average effect of mutations on trait values is also small.

We have assumed that the values of the parameters that describe the mutation process are within bounds that are currently believed to be biologically realistic. We have also assumed that, at every locus, the degree of bias, |b_i|, is not very large compared with the standard deviation of mutational effects, m_i. Under these assumptions, phenotypic variance and mean fitness are almost unaffected by the value of . On the other hand, the deviation of the mean phenotype from the optimum, | - z_opt|, is sensitive to the value of . This dependency turns out, however, to be nonmonotonic. While a small amount of bias (e.g., slightly in excess of zero) tends to move the value of in the direction of the bias, a point is always reached where any further increase in the degree of bias will actually bring the value of closer to z_opt (see Fig 1 for the special case of equivalent loci). Furthermore, for plausible parameter-value choices, and for all models of mutation studied here, the maximum-possible deviation of from z_opt is quite small: <1% of a phenotypic standard deviation. The small effect of bias depends on the existence of just one optimal phenotype, as assumed throughout this article. If multiple optima exist, then a small amount of bias may have very large long-term evolutionary effects (YAMPOLSKY and STOLTZFUS 2000 ).

The reason for the nonmonotonic response of to mutational bias is easiest to understand if one relaxes one of the assumptions of the model and considers the fate of mutations when the degree of bias is large at every locus (|b_i| very large at all loci). The extreme degree of bias leads to an equilibrium mean phenotype, , that is very close to z_opt, at least among adults. In this case, the reason is very clear. Nearly every mutation causes such a large deviation from the optimum that it induces death before maturation. Thus, only nonmutant offspring tend to survive, and so the adults have a mean phenotype that is very close to the optimum. When the values of the bias parameters are smaller, not every mutation induces fatality, and so mutant effects can accumulate each generation. Thus, the effect of mutational bias upon the equilibrium mean phenotype is largest when the degree of bias takes on an intermediate value.

Our claim that mutational bias cannot cause a large deviation of from z_opt depends on our assumption that the standard deviation of the mutant allelic effects at locus i, namely m_i, is not very small. This is embodied in the assumption m²_i >> µ_i/s. While the assumption of a substantial value of m_i is consistent with much of the data on mutation, it should be recognized that very-small-effect mutations are hard to identify, and so current estimates of m_i may be much too large (TURELLI 1984 ; LYNCH and WALSH 1998 ; DAVIS et al. 1999 ). If many trait-affecting mutations have very small effects, then the m_i may be much smaller than current estimates allow. The implication of small m_i is that mutational bias may have a substantial effect on (see Equation 3, Equation 6, and Equation 7, but recall that these results have been derived assuming the m_i are not small).

Another point to keep in mind is that our analysis applies to very large (effectively infinite) populations. It is possible that the long-term impact of mutational bias upon phenotype is greatly enhanced when population size is small. This is because of the action of genetic drift. In a finite population mutations have dynamics that are similar to those of strictly neutral alleles if their effect on fitness is small in comparison to the reciprocal of the effective population size (CROW and KIMURA 1970 ). Thus, if the population size is sufficiently small and selection is sufficiently weak, then biased mutation may be able to move the mean phenotype a considerable distance from z_opt before selection becomes sufficiently strong to stop it. Furthermore, it seems possible that this enhanced effect of biased mutation will be seen even if the population is very large, but is broken into small subpopulations, which are weakly connected by mutation.

In the past, models of quantitative genetics have typically assumed that the mutation process is not affected by any kind of genetic constraint. This means, for example, that the probability that a mutation will increase the effect of an allele upon a trait is independent of the premutation effect of the allele (this is the = 1 case of our model). Our results show that, when mutation is biased, a lack of genetic constraint typically (i.e., with nonequivalent loci) leads to the evolution of ever-more-extreme allelic values at every locus that affects the trait. This is so even if the degree of bias is very small and even if mutation is biased at only one locus, with mutations at all other loci being unbiased. The reasons for this continuous change in allelic effects lie within the nonlinear mathematics describing the problem. However, we note that the neutral equilibria, described in Equivalent loci with = 1 in RESULTS, lie at the heart of the phenomenon. It is evidently the case that for nonequivalent loci with = 1, the neutral equilibria acquire dynamical behavior at the level of the alleles, which themselves are not directly under selection but which underlie the trait. However, there is no manifestation of this dynamical behavior at the level of the trait—which is directly under selection.

Further insight into the processes responsible for the various results can be obtained by considering the details of the analysis, as presented in the Appendix. We note, however, that biased mutation generally induces skew into the distributions of allelic effects, with each distribution generally asymmetric about its mean. One exercise that is particularly instructive is to try to simplify the analysis, for the case = 1, by ignoring skew in the distributions of allelic effects that become established at long times. With the neglect of skew, the equilibrium value of - z_opt - z_opt is the outcome of a dynamical balance between a term proportional to the additional selection coefficient that is induced because - z_opt 0, namely 2s x V_G x ( - z_opt), and the mutational input into the mean phenotypic effect per generation, . Thus, under the neglect of skew, - z_opt is given by . Combining this with our (numerically verified) finding that moderate mutational bias causes negligible change in the genetic variance indicates that the mean phenotypic value at equilibrium, , no longer behaves nonmonotonically as the degree of bias is increased. Thus, ignoring skew leads to incorrect results. The behavior of the skew of the distributions of allelic effects is also strongly implicated in the continuous change in allelic effects exhibited when there are nonequivalent loci with = 1. Thus, some of the more intriguing behavior produced by the model depends on the skew in allelic distributions that is induced by mutational bias.

Of course, a model that leads to ever-increasing (or ever-decreasing) allelic effects is obviously not biologically reasonable. There must be some constraint on the effect that alleles, at any given locus, can have on a trait. Furthermore, biased mutation is known to occur in a variety of cases (SANTIAGO et al. 1992 ; LYMAN et al. 1996 ; MACKAY 1996 ; KEIGHTLEY and OHNISHI 1998 ), and it seems very unlikely that perfectly unbiased mutation is a common phenomenon. These considerations lead to the conclusion that reasonable models of mutation for quantitative trait loci must include some form of genetic constraint. The only reason that unconstrained mutation models, such as CROW and KIMURA's (1964) classic model, have been used in previous studies, without much difficulty, is that these studies either have made the biologically unlikely assumption of strictly unbiased mutation or have not addressed the inconsistency of assuming an equilibrium (this inconsistency of equilibrium is addressed in the Appendix).

In this study, we incorporated a simple model of genetic constraint that was suggested, in the absence of bias, by ZENG and COCKERHAM 1993 . The degree of constraint is characterized by the parameter . If = 1, then there is no constraint, while slightly <1 means that genetic effects can become quite extreme before genetic constraint has much effect on the distribution of mutants. If = 0, then the degree of constraint is maximal, and parental allelic values have no effect on mutant allelic values (this is the house-of-cards model of KINGMAN 1978 ).

We found that, whenever < 1, allelic effects tend to come to an equilibrium. The equilibrium distribution of allelic effects at each locus appears to be independent of initial conditions. However, if the degree of constraint is not very large ( slightly <1), then, at equilibrium, when there are more than one nonequivalent loci, the allelic effects at each locus tend, generally, to be rather extreme. This is because extreme effects at one locus tend to be compensated by extreme, but opposing, effects at other loci. Despite this, the equilibrium values of the mean phenotype and the genetic variance are, typically, not much affected by mutational bias.

The results raise some intriguing possibilities in the realm of molecular evolution. When mutations are constrained ( < 1) the equilibrium distributions of allelic effects, at the various loci, depend on the parameters of the model. These parameters can, quite plausibly, change over time. For example, a change in temperature might affect mutation rate (GROGAN et al. 2001 ) and/or the optimal value of the phenotype. We have undertaken very preliminary studies, and these have shown that with < 1, small changes in the parameters that govern selection or mutation may cause very large changes in the distribution of alleles at two or more of the n loci that control the trait. This effect is most dramatic when is close to unity.

Because changes at the two loci compensate each other, even very large changes in mean allelic effects are virtually invisible at the phenotypic level. However, the large changes in mean allelic effects imply that a great deal of molecular evolution has occurred at all loci involved. This molecular evolution can take a long time, as the time-to-equilibration seems to roughly scale as µ^-1, which is as one would expect. Thus, if the per-allele mutation rate is on the order of 10^-5, then, after a small environmental change, it may take many thousands of generations for loci to undergo large mutually compensatory changes before approaching equilibrium. We intend to explore these phenomena further in future studies. We also intend to explore the possible implications of this sort of "quasi-neutral" genetic change for speciation in a model of the sort studied by DOBZHANSKY 1936 and MULLER 1939 , MULLER 1940 .

To sum up, it seems that biased mutation may hold the key to understanding some of the phenomena that fascinate evolutionary biologists today. In a single large population, biased mutation seems to have little effect on phenotypes, but it might have a substantial long-term impact on molecular evolution. In a subdivided population biased mutation has the potential to magnify small environmental differences between the habitats of different subpopulations, and so it might lead to speciation and to related behavioral mechanisms of reproductive isolation. In the light of this, it would not be surprising if further investigations implicate biased mutation as a prime mover of evolution in additional areas that have not been anticipated here.

	ACKNOWLEDGMENTS

It is a pleasure to thank Nick Barton and John Welch for some very helpful discussions. We also thank two anonymous reviewers for suggestions that have significantly improved the manuscript.

Manuscript received March 6, 2003; Accepted for publication April 5, 2003.

	APPENDIX

TOP ABSTRACT MODEL RESULTS DISCUSSION APPENDIX LITERATURE CITED

Here we provide the theoretical background to the results presented in the main body of the article. We adopt the convention that unless specified to the contrary, all integrals range from - to .

In the model (of which a full description is given in the main text), viability of individuals of genotypic value G, namely w(G), has been taken as a quadratic function of G rather than as alternatives such as a Gaussian function. Nevertheless, if the strength of selection is small enough to satisfy sm²_i < 0.05, then there is not a substantial difference between results calculated from a Gaussian w(G) and those calculated from its quadratic approximation. In particular, we have considered the typical magnitude of the quartic term that is omitted in a quadratic approximation to a Gaussian w(G). We estimate that summary statistics, such as the mean phenotypic value, the genetic variance, and the genetic load, differ by <10% between the results calculated from a Gaussian w(G) and its quadratic approximation.

In all of the analysis of this work, we have followed the classic treatments by making the approximation of global linkage equilibrium (TURELLI and BARTON 1990 ). This holds to good accuracy when selection is weak in the sense sV_G << 1 where V_G is the genetic variance (BULMER 1989 ). We have numerically investigated the validity of the approximation of linkage equilibrium by numerical iteration of the full dynamical equations of a two-locus model with 17 discrete-effect alleles at each locus and thus in excess of 4 x 10⁴ different possible genotypes. This model is a direct discretization of the biased continuum-of-alleles model of mutation at the center of this work. We adopted typical parameter values for the simulations (µ_1,2 10^-5, m_1,2 0.2, s = 0.025, b_1,2 m_1,2, 0.9) and iterated the full dynamical equation for >10⁶ generations. The levels of linkage disequilibria observed, as measured by the correlation - , were <10^-3 of the allelic variances, ², at either locus in either the presence or the absence of bias. This is in accordance with what we should theoretically expect: the level of linkage disequilibrium generated by selection between any two unlinked loci depends on the product of the variances of allelic effects of the two loci (see, e.g., Equation 32 of WAXMAN 2000 ). Under the house-of-cards approximation (see below), allelic variance changes negligibly due to mutational bias, and hence the level of pairwise linkage equilibrium changes negligibly due to mutational bias. Thus the estimates in the literature for the very low levels of linkage disequilibria remain equally valid in the presence of bias.

Let us proceed with the analysis, under the assumption that the distribution of allelic effects in gametes approaches an equilibrium solution, which we write as (x). Numerical results are in accordance with this, except in the case of nonequivalent loci with = 1 and for this case we provide, below, a nonequilibrium analysis.

The approximation of global linkage equilibrium (BULMER 1989 ; TURELLI and BARTON 1990 ) leads to (x), being given by . Here _j(x_j) is the equilibrium distribution of allelic effects of maternal origin at locus j in offspring (it is identical to the distribution of allelic effects of paternal origin at the same locus) and is a nonnegative and normalized function: _j(x_j) 0, _j(x)dx = 1. The distribution _j(x_j) obeys an effective one-locus haploid equation that arises by integrating the equilibrium equation for (x) over allelic effects of all loci, with the exception locus i. The effective haploid equation reads

(A1)

Here an overbar denotes an average of the respective quantity over the distribution in zygotes, for example, . The quantity - _i appearing in Equation A1 represents a genetic background contribution that arises from all alleles except one of the alleles at locus i. The quantity µ_i is the allelic mutation rate of locus i and f_i(x_i - u - b_i) is the distribution of mutant allelic effects given in Equation 2 of the main text.

Equation A1 coincides, in form, with the equilibrium equations of KIMURA 1965 . It has the same general form as Equation 20 of BULMER 1989 and Equation 2.11 of TURELLI 1984 , when the latter is specialized to a quadratic fitness function, and terms of order µ_i x (selection coefficient) are neglected because allelic mutation rates are small and selection is weak.

We can write Equation A1 in the useful form

(A2)

where . Noting that f_i(x) f_i(0) it follows, from Equation A2 and normalization of _i(x_i), that _i x µ_i/(m_is). It is the smallness of _i, compared with m_i, that lies at the heart of the house-of-cards approximation (TURELLI 1984 ). The house-of-cards approximation applies when the strength of selection acting on allelic effects at any locus is much stronger than the effects of mutation at the locus, i.e., when µ_i/(sm²_i) << 1 (which is equivalent to _i << m_i). When this condition applies, the variance of mutant allelic effects is large compared with the equilibrium variance of allelic effects. Thus the allelic effect of a mutant is virtually unrelated to the parental allelic effect, and this is very similar to the exact behavior of the house-of-cards mutational model of KINGMAN 1978 , hence the name of the approximation. In the problem under consideration, we implement the house-of-cards approximation by replacing f(x_i - u - b_i)_i(u)du in Equation A2 by f(x_i - _i - b_i), leading to

(A3)

and the requirement of normalization determines the value of _i.

We now investigate Equation A3 for some particular cases.

Equivalent loci with = 1:
Consider the case where = 1 and all loci are mutationally equivalent; i.e., at loci affecting the trait, allelic mutation rates, mutational variances, and mutational biases are all given by µ, m, and b, respectively. Then Equation A3 takes the form

(A4)

(we omit the subscript i, on , in the case of equivalent loci). We proceed by multiplying Equation A4 by (x_i + - z_opt - _i) and integrating over all x_i. Using the substitution . Combining the integral with the integral with y -y yields

(A5)

In this form, we make two additional, but well-controlled approximations that lead to errors in - z_opt of O(µ²). The first approximation is to neglect ² within the integral. The rationale is that when |y| >> , the presence of is irrelevant. Furthermore, when y we can estimate, by Taylor expanding f(y - + z_opt - b) - f(y + - z_opt + b) to linear order in y, that neglecting ² results in an error of O(µ) = O(µ²). Thus neglect of ² is well justified in Equation A5. The second approximation is to note that - z_opt is, by Equation A5, O(µ) and hence neglecting - z_opt where it appears on the right-hand side of Equation A5 again leads to errors in - z_opt, on the left-hand side, of O(µ²). As a consequence, the mean phenotypic value, , which coincides with the mean genotypic value, , is given by - z_opt µ(2s)^-1 y^-1 [f(y - b) - f(y + b)]dy. Manipulations of this integral show

(A6)

where D(ß) is given by Equation 4 of the main text. In this way we arrive at Equation 3 of the main text.

An additional result follows by multiplying Equation A4 by (x_i + - z_opt - _i)², integrating over all x_i, and neglecting ², as justified above. This shows that, to leading order in µ, the allelic variance at any locus, , is given by µ/s, i.e., independent of locus label, i, and unchanged from its unbiased (b = 0) value. As a consequence, the genetic variance and genetic load are, to leading order in µ, unchanged from their unbiased values: V_G 2nµ/s and L = 1 - 2nµ.

Note that in the above calculations, the mean allelic effects, _i, do not appear in the final results and as a consequence are not determined by the equilibrium calculations. This is an exact property of Equation A1 for equivalent loci and = 1. It may be seen to follow from the change of variable , which eliminates _i from the equation. This property is a manifestation of the fact that in a dynamical calculation, the constant values the _i achieve at long times are dependent on initial data and numerical solution of the dynamical equations exhibits this feature.

Equivalent loci with < 1:
When = 1, and all loci are mutationally equivalent, the approximate distribution _i(x_i), following from (A3), takes the form _i(x_i) (µ/s)f(x_i - _i - b)[(x_i + - z_opt - _i)² + ²]^-1. Following closely the analysis for the case of equivalent loci with = 1, we find

(A7)

and a variance in allelic effects of . Thus genetic variance and genetic load are, to leading order in µ, unchanged from their unbiased values: V_G 2nµ/s and .

Note that for equivalent loci but with < 1, the _i cannot be eliminated from the equations by a coordinate transformation (unlike the case where = 1). In particular, _i appears explicitly in Equation A7 and numerical results verify that all _i are uniquely specified at equilibrium.

Equation A7 holds for i = 1, 2, ... , n and it is plausible that the equilibrium mean allelic effects, _i, are equal at all loci and given by /(2n) and this is numerically confirmed in all cases considered.

Since Equation A7 yields we can replace on the right-hand side of Equation A7) by . This introduces only errors of O(µ²) in - z_opt and using Equation 4 we obtain Equation 6 of the main text—which is simply Equation 3 of the main text with the "regression correction" b b - (1 - )z_opt/(2n).

Nonequivalent loci with < 1:
In the case where < 1, and loci are not mutationally equivalent, the approximate distribution _i(x_i), following from (A3), takes the form _i(x_i) (µ_i/s)f_i(x_i - _i - b_i)[(x_i + - z_opt - _i)² + ²_i]^-1. Proceeding as previously, we find numerically that the _i are not indeterminate and making the same approximations as previously, we find so , and

(A8)

where D(ß) is given in Equation 4. While similar in form to the result for equivalent loci with < 1 (Equation 6), we note that Equation A8 is substantially more complicated since there is no simple approximation for the _i. To determine the _i, it is necessary to simultaneously solve the set of equations, Equation A8 for i = 1, 2, ... , n, supplemented with . This is generally nontrivial because D(ß) and hence the set of equations are nonlinear. Despite this, it is possible to draw a general conclusion from Equation A8. Noting