Genetic Variation Maintained in Multilocus Models of Additive Quantitative Traits Under Stabilizing Selection
Reinhard Bürger, Alexander Gimelfarb


Stabilizing selection for an intermediate optimum is generally considered to deplete genetic variation in quantitative traits. However, conflicting results from various types of models have been obtained. While classical analyses assuming a large number of independent additive loci with individually small effects indicated that no genetic variation is preserved under stabilizing selection, several analyses of two-locus models showed the contrary. We perform a complete analysis of a generalization of Wright’s two-locus quadratic-optimum model and investigate numerically the ability of quadratic stabilizing selection to maintain genetic variation in additive quantitative traits controlled by up to five loci. A statistical approach is employed by choosing randomly 4000 parameter sets (allelic effects, recombination rates, and strength of selection) for a given number of loci. For each parameter set we iterate the recursion equations that describe the dynamics of gamete frequencies starting from 20 randomly chosen initial conditions until an equilibrium is reached, record the quantities of interest, and calculate their corresponding mean values. As the number of loci increases from two to five, the fraction of the genome expected to be polymorphic declines surprisingly rapidly, and the loci that are polymorphic increasingly are those with small effects on the trait. As a result, the genetic variance expected to be maintained under stabilizing selection decreases very rapidly with increased number of loci. The equilibrium structure expected under stabilizing selection on an additive trait differs markedly from that expected under selection with no constraints on genotypic fitness values. The expected genetic variance, the expected polymorphic fraction of the genome, as well as other quantities of interest, are only weakly dependent on the selection intensity and the level of recombination.

MANY quantitative characters in natural populations are apparently subject to stabilizing selection toward an intermediate optimum (e.g., Endler 1986). This means that extreme phenotypes have lower fitness than those near the population mean. Therefore, stabilizing selection is expected to exhaust genetic variation. This view has been substantiated by classical analyses based on the assumption that loci are independent and each allele contributes only an extremely small amount to the total genetic variance (Fisher 1930; Haldane 1932; Robertson 1956). Further support came from Wright’s (1935) study of the so-called quadratic optimum model in which two diallelic loci contribute additively to the character whose fitness deviates in a quadratic way from its maximum value. By contrast, most quantitative traits exhibit relatively high levels of genetic variability in nature.

This apparent contradiction has been a fundamental problem in evolutionary genetics, and several mechanisms that can potentially contribute genetic variability under stabilizing selection have been proposed and investigated. In principle, variation can be maintained either by mechanisms acting directly on the considered trait or as a side effect of polymorphisms that are independent of the observed character. Among the direct mechanisms are migration, mutation, frequency-dependent selection, genotype-environment interaction, and epistasis. Considerable progress has been made in elucidating the potential and the limitations of these mechanisms to maintain genetic variation (cf., Lewontin 1974; Lande 1975; Barton and Turelli 1989; Gimelfarb 1989; Falconer and Mackay 1995; Maynard Smith 1998).

However, relatively little is known about the ability of stabilizing selection per se to maintain polymorphisms and genetic variability in quantitative traits controlled by more than one locus acting additively. Conflicting results have been obtained, depending on model assumptions about the number of involved loci, magnitude of allelic effects, and linkage equilibrium.

All the classical analyses of additive quantitative traits assumed many loci in linkage equilibrium, with each locus having a very small effect on the trait. We review these first following Bulmer’s (1971, 1980) generalization and formalization of the approaches of Fisher (1930), Haldane (1932), and Robertson (1956). The model assumes that at each locus contributing to the trait, there are two alleles with sufficiently small effects on the phenotype, such that the density function of the trait in the subpopulation with genotype AiAj is simply shifted by the amount gij, where gij is the deviation of the genotypic contribution of AiAj from the mean (cf. Nagylaki 1984 for proof that, in the absence of linkage disequilibrium, this approximation is correct to first order). Assuming additionally that the phenotypic distribution is normal, Bulmer (1971) showed that the polymorphic equilibria are unstable under stabilizing selection with allelic frequencies converging to either zero or one. Therefore, no genetic variability can be maintained under these assumptions.

Further support for the view that stabilizing selection alone cannot maintain genetic variation in additive traits comes from analyses of mutation-selection-balance models that assume n loci in linkage equilibrium and with mutational effects drawn from a continuous probability distribution (cf. Crow and Kimura 1964; Lande 1975; Fleming 1979; Turelli 1984; Bürger 1986; Turelli and Barton 1990). Bürger and Hofbauer (1994) proved for such a model under quite general assumptions that the genetic variance at any equilibrium decreases to zero as the per-locus mutation rates decrease to zero.

Wright (1935, 1952) studied a model in which a quantitative trait is determined additively by two diallelic loci with equal and symmetric effects with respect to the optimum. The fitness of the trait was assumed to decrease quadratically as the phenotype deviates from its optimum value. He showed that no stable polymorphic equilibrium could exist under these assumptions, thus supporting the view that stabilizing selection depletes genetic variation. Barton (1986) generalized Wright’s quadratic optimum model to many unlinked equivalent loci and omitted the assumption that allelic effects are symmetric with respect to the optimum. He showed that for half of the possible positions of the optimum (within the range of genotypic values) selection maintains no genetic variation at all, while for the other half of the positions variation is maintained at exactly one locus.

In general, it might be expected that in a genetic system with very many loci contributing additively to a trait and the effects of alleles varying among loci, the optimum phenotype can be matched closely by the genotypic value of some homozygote and, hence, either no genetic variation will be maintained by such a system or at most a polymorphism will be maintained at just one locus. In contrast to the classical view and to the results discussed above, several analyses of two-locus systems that allowed for linkage or nonequivalent loci arrived at very different conclusions. Gale and Kearsey (1968) and Kearsey and Gale (1968) observed that in a two-locus diallelic additive model with a triangular fitness function, both loci can be stably polymorphic if their effects are sufficiently diverse. They also noted that the amount of diversity required to maintain the polymorphism decreases as linkage becomes tighter. Nagylaki (1989) investigated two-locus diallelic models of stabilizing selection under the assumption of linkage equilibrium, but with rather general fitness functions that have their optimum at the value of the double heterozygote and decrease monotonically and symmetrically from the optimum. For a large class of such fitness functions, both loci may be stably polymorphic if the ratio of the effect of the “major” locus to that of the “minor” locus exceeds a critical value. Gavrilets and Hastings (1993, 1994) investigated an extended version of Wright’s (1935) quadratic optimum model that includes linkage and disposes of the assumption of equivalent loci. They showed that for sufficiently tight linkage and unequal allelic effects, stable two-locus polymorphisms may exist. Thus, several analyses of two-locus systems suggest that substantial amounts of genetic variation can be maintained under stabilizing selection.

A number of interesting and unresolved questions emerge from the above results: How many loci have to contribute to an additive quantitative trait in order for no (or almost no) genetic variance to be maintained at equilibrium by stabilizing selection per se? How does this depend on the distribution of allelic effects and the strength of stabilizing selection? What is the role of linkage? More generally, what can be said about the equilibrium structure of multilocus systems under stabilizing selection on an additive quantitative trait?

In this article, we address these questions by studying diallelic two-, three-, four-, and five-locus models with arbitrary recombination and allelic effects varying between loci. We assume the quadratic optimum model in which the optimum coincides with the totally heterozygous genotype. For two loci, the equilibrium structure can be determined explicitly, while for three, four, and five loci, we follow the approach of Gimelfarb (1998) and perform numerical iterations of the recursion equations for large sets of randomly chosen recombination rates and allelic effects. In this way, we obtain information on how the expected equilibrium properties under stabilizing selection depend on the number of loci and can compare the properties of multilocus systems under stabilizing selection with the properties of systems having the same number of loci but randomly assigned genotypic fitness, i.e., under selection without constraints on the genotypic fitnesses.


In an infinite, randomly mating diploid population, a quantitative character that is controlled additively by n diallelic loci is considered. The contribution of one allele at each locus is zero, whereas the contribution, β , of the other allele is a random number between zero and one. It is assumed that the minimum and maximum genotypic values are always zero and one. Therefore, the actual contribution by the second allele at locus is scaled to be α=12βΣk=1nβk . This implies that the genotypic value of the total heterozygote is always ½, and the average allelic effect among the n loci controlling the trait is α=1(2n) . Environmental variance is ignored, so that genotypic values and phenotypic values are identical. We assume the quadratic optimum model, i.e., the fitness of an individual with genotypic value is G is W(G)=1s(G12)2. (1) The strongest possible selection occurs for s = 4, when the minimum and maximum phenotypes are rendered lethal. This normalization has the advantage that the strength of selection on genotypes can be compared for different numbers of contributing loci.

Let gametes be designated by i, their frequencies among zygotes in consecutive generations by pi and pi′, and the fitness of a zygote consisting of gametes j and k by Wjk. Further, let R(i|j,k) denote the probability that a randomly chosen gamete produced by a (j,k) individual is i. The function R is determined by the pattern of recombination between loci. With these ingredients, the well-known system of recurrence equations describing the dynamics of the distribution of gametes under viability selection followed by recombination is given by pi=1WjkpjpkWjkR(ij,k), (2) where W = Rj,kpjpkWjk denotes mean fitness.


For two loci, fairly complete analytical results of the above general model can be obtained. If the alleles at the first and second locus are denoted by A1, A2, and B1, B2, respectively, the genotypic values of the four possible gametes, A1B1, A1B2, A2B1, A2B2 are 0, α2, α1, and α1 + α2 = ½, respectively. For definiteness, we assume α1 ≥ α2 and refer to these loci as major and minor, respectively. With the fitness assignments as above, the fitness values of the nine possible genotypes are given by A1A1A1A2A2A2(1d1b1a1c11c1a1b1d),B1B1B1B2B2B2 (3) where a = s1 - α2)2, b=sα12 , c=sα22 , d = ¼s. This is the symmetric viability model that has been intensively studied (see Karlin and Feldman 1970 for an extensive equilibrium analysis and further references). The parameters satisfy the additional relation a + d = 2(b + c), which is of technical importance. Wright (1935, 1952) was the first to investigate this model of stabilizing selection. More recently, Gavrilets and Hastings (1993, 1994) calculated all equilibria and explored the transient dynamics of the mean genotypic value and the genetic variance.

Our aim here is to study not only the equilibrium structure of this model, but in particular the amount of genetic variance that can be maintained through selection, as well as the magnitude of further quantities, like mean fitness, deviation of the mean phenotype from the optimum, and linkage disequilibrium. Proofs of the results presented below are given in the appendix, and a graphical representation is shown in Figure 1.

There are four types of equilibria: (i) a fully polymorphic equilibrium that is symmetric with all allele frequencies equal to 1/2 (in the appendix, this is denoted by F1); (ii) a pair of fully polymorphic equilibria that, in the language of the symmetric viability model, are called unsymmetric because their coordinates satisfy no simple symmetry relations (they are denoted by F2 and F3); (iii) a pair of equilibria (denoted F4 and F5) with the major locus polymorphic and the minor locus fixed; (iv) two monomorphic equilibria corresponding to fixation of A1B2 or A2B1 (these are denoted by F6 and F7).

The existence and stability of these equilibria depend on the relation between the parameters α21 and r/s, where r is the recombination rate. For consistency with the multilocus results, we use the standard deviation of allelic effects, σα = ½(α1 - α2) = ¼ - α2, as a measure for the disparity of effects. An important value is σα=112 , which is equivalent to α1 = 2α2. The ranges of stability are displayed in the bottom of Figure 1 for r/s ≤ ⅛, which covers the full range of parameters if selection is as strong as possible (s = 4). For weaker selection, larger values of r/s can occur. Fully polymorphic equilibria can be stable only if rs112 . Only one class of equilibria can be stable for a given parameter combination. The unsymmetric equilibria are stable whenever they exist, while the other equilibria always exist but may be unstable.

For each equilibrium, mean fitness, absolute deviation of the mean phenotype from the optimum, genetic variance, and linkage disequilibrium can be calculated (see appendix). For a range of recombination rates, the top and middle parts of Figure 1 display the equilibrium genetic variance and the deviation of the mean phenotype from the optimum as function of σα.

The following is a summary of the main findings:

  1. Substantial genetic variance is maintained for any recombination rate if the effects of the alleles at the two loci are sufficiently different (σα > 1/12), in which case one or both loci, depending on the recombination rate, are stably polymorphic. The reason is that no homozygote genotype is close to the optimum and, from double heterozygotes, a high recombination rate would generate many complete homozygotes with extreme phenotype and low fitness. Therefore, if recombination is strong relative to selection, the minor locus is fixed and the major locus is keptpolymorphic by overdominance. With tight linkage, both gametes A1B2 and A2B1 are maintained at high frequency.

  2. If the effects of the two loci are not very different (σα < 1/12), very tight linkage (r < 4sσ2α) is necessary to maintain genetic variation. In this case, a homozygote (A1A1B2B2 or A2A2B1B1) is close to the optimum and, therefore, fixed unless linkage is extremely tight.

  3. For any value of r, the equilibrium genetic variance decreases from its maximum value 12(α12+α22)=18 to 0 as σα decreases from its maximum value ¼ to 0.

  4. For fixed σα, the genetic variance has a maximum value for intermediate recombination rates. If the effects are not very different (σα < 1/12), then the variance is lower under loose linkage than under very tight linkage, while the reverse is true if the effects are sufficiently different. Thus, tighter linkage may increase the degree of polymorphism at equilibrium but at the same time reduce the genetic variance.

  5. The mean phenotype evolves rapidly to its equilibrium value (Hastings 1987). However, for unequal effects the equilibrium does not coincide with the optimum unless linkage is extremely tight. Convergence of the genetic variance and of the genotype frequencies may be very slow (Gavrilets and Hastings 1993, 1994). Therefore, if there is substantial initial genetic variability and the selection regime is such that it eventually depletes genetic variability, this process may be very slow even for relatively strong selection.

Figure 1.

—Equilibrium properties of two-locus two-allele systems under stabilizing selection. (Top) The ratio of the equilibrium genetic variance, Formula to the maximum possible variance, Formula , as a function of the standard deviation of allelic effects, σα = ¼ - α2, for five different rates of recombination and s = 4. (Middle) The corresponding absolute deviation of the mean from the optimum. (Bottom) The regions of stability of the four types of equilibria. The following abbreviations are used for the regions of stability. 0, two monomorphic equilibria; 1, two equilibria with one (the major) locus polymorphic; 2a, the symmetric equilibrium with both loci polymorphic; 2b, two asymmetric equilibria with both loci polymorphic. The upper boundary of region 2a is given by r1/s (A9) and the upper boundary of 2b by r2/s (A10). At σα, we have Formula . The boundary between the regions 0 and 1 is at Formula and intersects with region 2b at r/s = 1/36. The horizontal lines correspond to the curves in the top and middle parts. It should be noted that the areas indicating the different types of equilibria are not proportional to their respective probabilities, because the horizontal scale is σα. It may be noted that 0 ≤ r/s1/8 covers the whole range only if s = 4. For smaller s (weaker selection), the maximum value of r/s is larger than 1/8; hence the proportion of parameter values yielding stable two-locus polymorphisms becomes correspondingly smaller.


Usually, parameters of genetic systems controlling quantitative traits are unknown or can be inferred only indirectly. In addition, because the dimensionality of the parameter space and of the space of gamete frequencies increases rapidly as the number of loci increases, an explicit and analytical characterization of the equilibrium properties of multilocus models in terms of all parameters and initial conditions would be of limited value, even if the necessary analytical methods were available. Therefore, we used a different approach to evaluate the quantities of interest for randomly chosen parameter sets and initial conditions and, consequently, to obtain statistical results.

For a genetic system with a given number of loci (n = 2, 3, 4, 5), we constructed 4000 parameter sets (allelic effects of loci, coefficients of recombination between adjacent loci, the strength of stabilizing selection). For each parameter set, allelic effects were obtained by generating values β (= 1, 2,..., n) as independent random variables uniformly distributed between 0 and 1 and transforming them into the actual allelic effects α=12βΣkβk . The strength of stabilizing selection, s, was obtained as a random variable uniformly distributed between 1 and 4, with s = 4 corresponding to the strongest possible quadratic selection (the fitness of extreme phenotypes is zero), while s = 1 represents “weak” selection (the fitness of extreme phenotypes is 75% of the fitness of the best fit phenotype). On the basis of such obtained values of allelic effects and the strength of selection, genotypic fitnesses, Wij, were calculated and substituted into the recursion equations (2). Recombination coefficients between adjacent loci, r,+1(=1,,n1) , were obtained as independent random variables (no interference) uniformly distributed between 0 and 0.5. For four-locus genetic systems, we also constructed 4000 parameter sets with random allelic effects and random recombination, but with “strong” (s = 4) and weak (s = 1) selection, as well as with random allelic effects and random selection but with free recombination, r,+1=12 . In addition, we constructed 4000 parameter sets for two-, three-, four-, and five-locus systems with random recombination between adjacent loci and with genotypic fitness values, Wij, chosen as independent random variables uniformly distributed between 0 and 1, as described by Gimelfarb (1998).

For each of the 4000 parameter sets, the recursion equations (2) were numerically iterated starting from 20 random initial distributions of gametes. To make the initial distributions more evenly distributed in the gametic space, they were chosen such that the (Euclidean) distance between any two of them was no less than a predetermined value (0.25, 0.30, 0.35, and 0.30 for 2, 3, 4, and 5 loci, respectively). Starting from an initial distribution, the recursion equations (2) were iterated until either equilibrium was reached or the number of iterations exceeded 300,000. In the latter case, the parameter set was excluded from the analysis. For each parameter set, the number of different equilibria, the gametic frequencies at each equilibrium, and the number of trajectories (initial distributions) converging to each equilibrium were recorded. Using this database, the equilibrium properties of multilocus genetic systems were analyzed.

For two-locus genetic systems, results can be obtained not only by the statistical approach, but also by numeric integration of equations in the appendix, as well as by placing a fine grid over the parameter space and averaging over the appropriate quantities whose equations are given in the appendix. All three methods yielded, within statistical and numeric accuracy, identical results.


Results obtained by the statistical approach are summarized in Tables 1, 2, 3 and in Figures 2 and 3.

Table 1 shows the proportion of parameter sets for which there was at least one equilibrium with a given number of polymorphic loci (e.g., among 4000 of four-locus systems under stabilizing selection, 68% yielded at least one monomorphic equilibrium, while only 1% of such sets yielded at least one equilibrium with two polymorphic loci). Also shown in Table 1 are the proportions of trajectories that, starting from a random initial distribution, converge to an equilibrium with a given number of polymorphisms (e.g., 58% of trajectories for four-locus genetic systems under stabilizing selection converged to a monomorphic equilibrium).

Table 2 presents some of the parameters characterizing stable equilibria that are expected for genetic systems with a given number of loci: the number of different stable equilibria per parameter set; the polymorphic fraction of the genome (the probability for a locus to be polymorphic); and the mean fitness. For stabilizing selection, Table 2 also presents the expected deviation of the mean from the optimum; the genetic variance at equilibrium; the ratio of the genetic variance at equilibrium to the maximum genetic variance that can be maintained by the given genetic system; and the ratio of the allelic effect at a polymorphic locus to the average allelic effect among all loci in the system. The maximum genetic variance that can be maintained in linkage equilibrium by an additive trait controlled by n loci with allelic effects {α} is Vmax=12Σα2 , while the average allelic effect among loci is α=1(2n) .

Table 3 compares parameters of stable equilibria for four-locus genetic systems with different recombination rates and under stabilizing selection of different strength.

Equilibrium structure: The data in Table 1 demonstrate that, while the probability for a genetic system under stabilizing selection to maintain a monomorphism or a one-locus polymorphism is higher for systems with more loci, the probability of maintaining more than one polymorphic locus drops rapidly for systems with more loci. In fact, for a five-locus genetic system, the probability of maintaining a polymorphism in two or more loci is negligible. The same is not true under selection with randomly assigned fitnesses, in which case the probability of maintaining more than one locus polymorphic is higher for genetic systems with more loci. In addition, Table 2 shows that with increasing number of loci, the polymorphic fraction of the genome decreases at a much higher rate under stabilizing selection than under selection with random fitness. Indeed, as a function of n, the polymorphic fraction of the genome under stabilizing selection is very closely approximated by 0.99n-1.58, while for random fitnesses it decreases proportionally to n-0.45.

For two loci, we proved in the appendix that different types of equilibria cannot stably coexist (except for r = 0), and the maximum number of coexisting stable equilibria is two. It can be noted that the sums of the entries in the column “sets” for stabilizing selection in Table 1 are 1.00, 1.02, 1.33, and 1.59 for 2, 3, 4, and 5 loci, respectively. This indicates that the probability of several simultaneously stable equilibria with different degrees of polymorphism, while zero for 2 loci and very low for 3 loci, becomes substantial if the number of loci increases. Also, Table 2 shows that with increasing number of loci in the genetic system, the expected number of simultaneously stable equilibria increases faster under stabilizing selection than with randomly assigned fitnesses, and for 4 and 5 loci, it is actually higher under stabilizing selection. The maximum numbers of stably coexisting equilibria detected by numerical iteration starting with 20 initial distributions were 6, 12, and 15 for n = 3, 4, and 5, respectively, and there was exactly 1 polymorphic locus at all these equilibria. These results indicate that historic effects may be of paramount importance in the evolution of populations under stabilizing selection.

View this table:

Equilibrium structure

View this table:

Stable equilibria and their genetic variability

View this table:

The role of selection intensity and recombination

Because recombination rates between adjacent loci in our parameter sets were generated as independently and uniformly distributed random variables, the probability of obtaining an n-locus system with all loci tightly linked decreases exponentially as n increases. Therefore, we also performed some iterations for four-locus genetic systems with completely linked loci. All trajectories converged to fully polymorphic equilibria. However, only two complementary types of gametes, e.g., 0110 and 1001, were present at these equilibria and at equal frequency. Such equilibria correspond to the symmetric equilibrium F1 in the two-locus case. Sometimes convergence to such equilibria was extremely slow. By perturbation arguments (Karlin and McGregor 1972), these findings should extend to sufficiently small recombination rates.

Linkage disequilibrium: The average amount of linkage disequilibrium was calculated for equilibria with at least two polymorphic loci and was found to be very small (data not shown). Given that the probability of a stable equilibrium with two or more polymorphic loci is very small if more than two loci control the trait, linkage disequilibrium can be neglected in multilocus systems under stabilizing selection.

Strength of selection and recombination: As is well known and shown by Figure 1, for two loci the equilibrium structure depends on the relative magnitude, r/s, of recombination rate and selection intensity. For a given standard deviation of allelic effects, σα, the equilibrium genetic variance is not very sensitive to changes in the recombination rate (relative to the selection intensity) except for small σα and very small r. Furthermore, for any particular value of s, the equilibrium genetic variance, VG, averaged over all trajectories and all parameter sets {r, α1, α2} is almost independent of s. It decreases from VG = 0.0311 for s = 4 to VG = ln32 - ⅜ ≈ 0.0305 as s → 0. However, as shown by Figure 1 and discussed in analytic theory for two loci, the number of polymorphic loci at a stable equilibrium is strongly effected by r/s.

Employing the statistical approach, we investigated the effect of the strength of stabilizing selection and of recombination in the case of four loci. The results are summarized in Table 3 and demonstrate that, indeed, many expected equilibrium properties, notably the average number of equilibria, the average degree of polymorphism, and the average genetic variance at equilibrium, are almost independent of selection intensity and recombination rate. This is further supported by additional iterations with s = 4 and r = 0.5 (not shown). Presumably, this is so because most genetic variance comes from single-locus polymorphisms that are stable because of overdominance. In this case, allele frequencies are independent of s but depend only on the relative disadvantage of the two homozygotes. Different results are observed only if all loci are completely (or very tightly) linked, in which case, as mentioned above, all trajectories converge to fully polymorphic equilibria with gametes of only two complementary types present.

The explored values of selection intensities (1 ≤ s ≤ 4) cover the range from weak to the strongest possible selection but not extremely weak or no selection. Inclusion of values s < 1 would change our results very little because it would not much increase our parameter range and because all our results indicate that the equilibrium structure is nearly unaffected by changes in s unless selection is very strong relative to recombination. In addition, Tables 2 and 3 show that even strong selection imposes only a small genetic load on the population. Thus, purely on the basis of observation, one might conclude that these populations are under weak selection even if, in fact, selection is strong.

Figure 2.

—Ratio of the expected equilibrium genetic variance, VG, to the maximum genetic variance under linkage equilibrium, Vmax, as a function of the standard deviation of allelic effects, σα. Each data point represents the average over all trajectories for the given parameter set. Allelic effects, recombination rates between adjacent loci, and strength of stabilizing selection are independently drawn from uniform distributions as described in the statistical approach.

Deviation from the optimum: In general, the mean genotypic value deviates from the optimum. This deviation is displayed in Figure 3. It may be substantial for certain sets of parameters (see also below). However, as shown in Table 2, the average deviation from the optimum decreases as the number of loci increases. This is due to the fact that the optimum can be matched more closely in systems with more loci of varying effects. The present data quantify this statement and show that the deviation from the optimum declines slightly faster than the average allelic effect α . The ratio of the expected deviation from the optimum to the average allelic effect is 0.29, 0.31, 0.26, and 0.22 for n = 2, 3, 4, and 5, respectively.

Figure 3.

—Deviation of the expected mean genotypic value from the optimum as a function of the standard deviation of allelic effects, σα.

Genetic variance: Given that our scaling is such that the average allelic effect, α=1(2n) , is smaller for systems with more loci controlling a trait, it is expected that the genetic variance at equilibrium, VG(n), averaged over all trajectories and parameter sets (s,{r},{α}) , must decline with increasing number of loci, n. Table 2 shows that this indeed is the case. However, the decline occurs at a rate much faster than expected. To a close approximation, we have VG(n) = 0.52n-4.0. If the genetic variance is scaled relative to the maximum genetic variance that can be maintained for a given genetic system in linkage equilibrium, Vmax=12Σα2 , the resulting overall average value of VG(n)/Vmax [with VG(n) being the variance for a given parameter set, averaged over all trajectories] still decreases in proportion to n-2.4.

What is the reason for such a fast decline? First, we have already seen that under stabilizing selection the polymorphic fraction of the genome decreases with increasing number of loci at a rate of approximately n-1.58, which is much higher than under selection with random fitnesses. Second, if a locus maintains a polymorphism, i.e., contributes to the genetic variance, the allelic effect of such a locus (as compared to α , the average effect among all loci) is smaller for systems with more loci. If a trait is controlled by two loci then, as we have seen above, it is either the major locus or both loci that are segregating at a polymorphic equilibrium. If n = 3 or 4, Table 2 shows that the effect of a polymorphic locus is expected to be approximately equal to the average allelic effect, α , whereas if n = 5, it is expected to be smaller than the average (0.85α ). For comparison, if n = 5, the expected minimum allelic effect among polymorphic loci is 0.4α , while the expected maximum effect is 1.6α . It is also interesting to note that, in contrast to the two-locus case, for n ≥ 3 the expected genetic variance at an equilibrium is higher if this equilibrium maintains several loci polymorphic. For example, with n = 4, the expected genetic variance at equilibria with two polymorphic loci is almost four times as large as at equilibria with one polymorphic locus (data not shown).

Figure 2 displays the ratio VG/Vmax averaged over all trajectories for a given parameter set as a function of the standard deviation, σα, of allelic effects in the given set. The maximum genetic variance is attained if σα is maximal, which is the case if one locus has maximum effect (= ½), while the others have zero effect. The maximum value of σα is σα,max=(14n)(11n) , which is 0.25, 0.236, 0.217, 0.2 for n = 2, 3, 4, 5, respectively. If σα is maximum, then both alleles at the major locus are segregating at equal frequency and the mean phenotypic value coincides with the optimum.

Let us next discuss the properties of stable equilibria, in particular the equilibrium genetic variance and the deviation of the mean from the optimum as affected by the diversity of allelic effects measured by σα. First, we concentrate on two-locus systems with loci not very tightly linked (r/s > 1/36, which is the value at which the regions 0, 1, and 2b in the bottom of Figure 1 intersect). In such a system, the major locus or both loci are segregating if 112σα14 . The genetic variance decreases from its maximum value to zero, as σα decreases from 1/4 to 1/12, while the deviation from the optimum increases from zero and to its maximum value at σα=112 (Figure 1). If σα<112 , then only monomorphic equilibria are stable unless linkage is very tight. The smaller σα is, i.e., the more similar effects of loci are, the more closely the optimum can be matched by a completely homozygous genotype (the one that gives rise to the stable equilibrium). Therefore, the deviation from the optimum converges to 0 as σα → 0.

On the basis of this explanation, the two-locus panel of Figure 2 is easily interpreted. The thick line of densely packed points represents parameter sets with loose linkage because, for given σα, all two-locus systems with sufficiently large r/s have the same stable equilibria. Tightly linked loci are represented by the scattered dots.

With more than two loci, the complexity increases. Figures 2 and 3 show that there is no monotone relation between the standard deviation of allelic effects and the genetic variance. However, they also indicate that for sufficiently large σα, the equilibrium behavior is similar to the two-locus case. As σα decreases from σα,max, the effect of the major locus decreases from its maximum value of 1/2, and alleles with smaller effects are present at other loci. Therefore, the genetic variance declines. In analogy with the two-locus case, we suspect that if the loci are loosely linked, then it is the major locus that is polymorphic while the other loci are monomorphic. In the parameter range below σα,max, the mean increasingly deviates from the optimum as is illustrated by the almost straight lines ascending from σα,max to the left in the two-, three-, and four-locus panels of Figure 3. The deviation is caused by the asymmetries introduced by the monomorphic loci, because for such values of σα, no homozygote genotype will closely match the optimum. The scattered points below these lines represent parameter sets with linked loci. For sufficiently tight linkage, two or more loci may be polymorphic and, as in the two-locus case, at such equilibria the mean genotypic value may be close to or coincide with the optimum. This, however, is rarely the case if the number of loci is larger than three (see Table 1 and Figure 3).

If the number of loci is greater than two, a problem arises in that the same value of σα can be caused by different sets of allelic effects, and the smaller σα is, i.e., the more similar the effects of loci are, the more important this problem becomes. As σα decreases further, the genetic variance declines and eventually becomes 0. As shown in Figure 2, there is an interval of values of σα where this happens (for n = 3, 4, 5, this interval is contained in the range 0.11 < σα < 0.14). It is the same interval where the maximum deviation of the mean from the optimum occurs (cf. Figure 3). For a range of smaller σα values, the figures indicate that most stable equilibria are monomorphic and that the deviation of the mean now decreases. The few data points in this range of σα values that have a positive genetic variance come from tightly linked loci; their number decreases as n increases. As σα decreases further (below ≈ 0.1), the situation becomes more complicated as can be seen from the figures: obviously, for given σα, there are now many different types of stable equilibria; some maintain high levels of genetic variance (because they are polymorphic in more than one locus or because an allele of large effect is segregating), while most maintain relatively little variance. There also emerges another peak of the mean’s deviation from the optimum.

Finally, if all loci have very similar effects, i.e., if σα is close to zero, the equilibrium structure is simple again. If the number of loci is even, then all loci are fixed because there exists a homozygote genotype very close to the optimum, and VG/Vmax converges to 0. If, on the other hand, the number of loci is odd, then one of the (almost equivalent) loci is segregating, and the genetic variance converges to 1/(8n)2, i.e., VG/Vmax converges to 1/n, as observed in Figure 2. In both cases, the mean coincides with the optimum because of the symmetries of the model.


While previous analyses, assuming a quantitative trait controlled additively by a large number of loci of individually very small effects, indicated that stabilizing selection on such a character depletes genetic variation, several two-locus analyses showed the contrary (see Introduction). Because analyses of the first kind are based on assumptions such as linkage equilibrium and equivalent loci, it has not been obvious if their conclusions would still be valid without these assumptions. We demonstrated that under our model assumptions, which include linkage and loci with different effects, the expected genetic variance maintained under stabilizing selection declines very rapidly from a high value in two-locus systems to an extremely low value in five-locus systems, thus providing quantitative support for the classical results. In particular, already four- and five-locus systems exhibit equilibrium properties as expected under the infinitesimal model. In addition, with more than three loci, the probability of equilibria involving at least two polymorphic loci is almost negligible, implying that in such systems no linkage disequilibrium is to be expected. A further interesting finding is that many quantities, e.g., the average number of stable equilibria, the average polymorphic fraction of the genome, the average deviation of the mean phenotype from the optimum, and the average genetic variance are virtually independent of the strength of stabilizing selection and the level of recombination.

What are the implications of these findings for the classical problem of how genetic variation is maintained in traits that apparently experience stabilizing selection? Obviously, our results suggest that it is very unlikely, although not impossible, that appreciable levels of genetic variation are maintained at equilibrium by the selective forces resulting from direct stabilizing selection if the loci controlling the trait act additively. From this, however, no support can be deduced for any of the theories about the maintenance of genetic variation, be it by mutation-selection balance or by pleiotropic effects of polymorphisms maintained independently of the observed character.

First, our model provides a kind of null hypothesis in assuming that allelic effects at loci as well as recombination rates between “adjacent” loci are drawn from a uniform distribution. There is empirical evidence that allelic effects have a highly leptokurtic distribution (Mackayet al. 1992; López and López-Fanjul 1993), and it has been suggested that a (reflected) gamma distribution may be adequate to approximate distributions of mutational effects (Hill and Rasbash 1986). How would this change our results? Probably just a little. With such a distribution, the expected standard deviation of allelic effects is higher than under the present uniform distribution, while the mean, α , remains unchanged. [For instance, with n = 5, the expected σα for uniformly distributed β is 0.057, while for exponentially distributed β , it is 0.082. For n → ∞, it can be shown that the expected σα converges to σβ(2nβ) , where σβ is the standard deviation of the distribution of the β , and β is its mean.] Therefore, there would be more data points in Figures 2 and 3 with larger σα; in particular, it would be more likely to obtain a system with one major and many minor loci. As indicated by Figure 2, this might increase the expected genetic variance at equilibrium by increasing the average effect of a polymorphic locus (Table 2), but it could not offset the decline caused by the rapidly decreasing polymorphic fraction of the genome.

Little is known about linkage relations between loci affecting quantitative traits. If quantitative trait loci tended to be tightly linked, polymorphisms at several loci could be maintained. Except for the two-locus case, this would lead to higher levels of genetic variance.

Our analysis, as well as most previous analyses, is based on the assumption of additive loci. However, recent experimental results indicate that typical quantitative traits may be under control of loci with substantial epistatic effects (Mackay and Fry 1996; Routman and Cheverud 1997). Theoretical analyses of multilocus models under stabilizing selection have shown that with epistasis, high degrees of genetic variation can be maintained (Gimelfarb 1989). Future theoretical investigations of multilocus systems underlying quantitative traits should definitely pay more attention to epistatic effects. Because to date, only little is known about the nature of epistatic effects, such investigations had to be based on ad hoc assumptions.

Following the tradition in population genetics, we have explored the equilibrium properties of our model. But are natural populations ever in or close to equilibrium? The more loci are affecting a trait, the more likely it is that transient polymorphisms are contributing to genetic variation. Such transient polymorphisms can be present for a variety of reasons, for instance, new mutations sweeping through the population or balancing selection caused by forces that may be independent of the character under consideration. We have no quantitative information about the transient properties of our model populations, but we have observed that in some cases, particularly for tightly linked loci, convergence to equilibrium may be exceedingly slow, on the order of tens of millions of generations. Therefore, it might prove useful to devote more consideration to evolving populations and to study nonequilibrium properties.

Keeping all these reservations in mind, let us compare the amount of genetic variation that can be maintained under stabilizing selection alone with that maintained under mutation-selection balance. It is well known that the equilibrium genetic variance that can be maintained under mutation-selection balance in a large population is approximately 2nμ/s, provided the loci are not very tightly linked and the expected squared effect, α2=(1n)Σα2 , is larger than 10μ/s (cf. Latter 1960, Bulmer 1972; Turelli 1984; Bürger and Lande 1994, cf. Bürger 1998, for a review). In our simulations, we have α2=0.0132 if n = 5, and α2 is larger if n is smaller. If there are only two loci and if the character is not nearly neutral, then for all realistic mutation rates, the equilibrium genetic variance maintained by mutation-selection balance is much lower than the average value, VG = 0.031, that is maintained under pure stabilizing selection. If there are five loci, however, then the average genetic variance maintained in our model is only approximately VG ≈ 0.001. Such an amount of variance can be easily maintained by mutation-selection balance, namely whenever μ/s > 0.0001. This argument is not intended to imply that mutation-selection balance is the primary mechanism for maintaining genetic variation in quantitative traits (see Barton and Turelli 1989 and Bürger and Lande 1994 for discussions). It shows, however, that pure stabilizing selection can be excluded as an important mechanism, unless genetic systems are very different from those assumed in this article.

In accordance with previous theoretical investigations (cf. Turelli and Barton 1990; Bürger and Lande 1994), our results show that the average level of linkage disequilibrium is exceedingly low. The main reason is that with more than three loci affecting the trait and independent, uniformly distributed recombination rates, polymorphisms involving two or more loci are so rare.

Finally, we draw attention to the multitude of stable equilibria that may exist under stabilizing selection. We not only found that the expected number of coexisting stable equilibria is higher under stabilizing selection than under selection with randomly assigned fitnesses (see Table 2), but we also found that in a five-locus system, up to 15 of 20 randomly chosen initial conditions converged to different equilibria. This implies that historic effects may play an important role for populations under stabilizing selection. Therefore, it cannot be expected that closely related (but isolated) populations that have experienced similar selective regimes are genetically identical at all or most loci contributing to traits under stabilizing selection, even if the populations are large so that effects of random drift are negligible.


Let x1, x2, x3, x4 denote the frequencies of the four gametes, A1B1, A1B2, A2B1, A2B2, respectively. Collectively, all quadruples (x1, x2, x3, x4) satisfying 0 ≤ xi ≤ 1 and Rixi = 1 form the 4-dimensional simplex S4. This is the state space for the two-locus two-allele dynamics. Let D = x1x4 - x2x3 be the linkage disequilibrium measure and let Wi denote the marginal fitness of gamete i. Then, the well-known recursion equations can be written as Wxi=xiWiηirD,i=1,,4, (A1) where η1 = η4 = -η2 = -η3. The equations in this appendix are written without the constraint α1 + α2 = 1/2. However, it is assumed that the fitness optimum coincides with the value of the double heterozygote, i.e., W(G) = 1 - s[G - (α1 + α2)]2.

Equilibria: We assume r > 0. Then there exist up to nine equilibria, seven of which may be stable (but not simultaneously). We denote the possibly stable ones by F1-F7.

There always exists a symmetric equilibrium, F1, which is calculated to be given by F1:x^1=x^4=14+D^1,x^2=x^3=14D^1, (A2a) where D^1=14sα1α2[rs2α12α22+r2]. (A2b)

Two further interior equilibria, F2 and F3, may exist that are unsymmetric. They can be calculated explicitly using the results of Karlin and Feldman (1970). Let us introduce the coordinates x=x1x4,y=x2x3,z=x1+x4x2x3. (A3) Then we have x1=1+x+2x4,x2=1z+2y4,x3=1z2y4,x4=1+z2x4. (A4) Because of the specific form of the selection parameters a, b, c, and d in (3), the quadratic polynomial in Equation A7 of Karlin and Feldman simplifies dramatically (with a little help from Mathematica; Wolfram 1991) and becomes a linear expression in z. Its root, z^ , is given in (A6c) below. Implementing this value into Equation A2 of Karlin and Feldman yields two possible values for the equilibrium mean fitness at the prospective unsymmetric equilibria. It turns out that only the value We=198r+14rs(α12+α22)18s2(α12α22)2 (A5) leads to equilibria in S4. Inserting these expressions in (A3) and (A4) of Karlin and Feldman (1970) yields the desired solutions. These are F2:x^=(α1+α2)[rs(α1α2)2]R8s12r32α1α2 (A6a) y^=(α1+α2)[r+s(α1+α2)2]R8s12r32α1α2 (A6b) z^=(α1α2+α2α1)+3r4sα1α2+sα1α24r(α1α2α2α1)2, (A6c) where R=3r2+2rs(α12+α22)s2(α12α22)2. (A7) Substituting x^ , y^ , and z^ into (A4) gives the coordinates of the unsymmetric equilibria in terms of the chromosome frequencies x1, x2, x3, and x4. The equilibrium F3 is symmetric to F2 upon interchanging x^1 with x^4 and x^2 with x^3 (or, equivalently, taking - x^ and - y^ ). It is not difficult to show that F2 and F3 exist if and only if r1<r<r2, (A8) where r1=13s(α12+α22)+23sα14α12α22+α24 (A9) is the positive root of the equation R = 0 and r2=min{s(α1α2)2,13s(α12α22)}. (A10) The reader may observe that s(α1α2)213s(α12α22) if and only if α12α2 and 13s(α1α2)2<r1<13s(α12α22). (A11) We also note that the linkage disequilibrium at the unsymmetric equilibria is negative and given by = 14(z^+y^2x^2) ). This expression can be factorized and is negative between r1 and r2.

Next, there may exist two edge equilibria, F4 and F5, with the major locus polymorphic: F4:x^1=x^3=0,x^2=12+α2α1,x^4=12α2α1, (A12) F5:x^2=x^4=0,x^1=12α2α1,x^3=12+α2α1, (A13) These exist (i.e., are in S4) if and only if α1 > 2α2.

Finally, there are four corner (vertex) equilibria corresponding to fixation of one of the chromosomes. These always exist but only the vertices corresponding to fixation of A1B2 or A2B1 can be stable. These are denoted by F6 and F7, respectively, and their coordinates are F6:x^2=1,x^1=x^3=x^4=0, (A14) F7:x^3=1,x^1=x^2=x^4=0. (A15)

If r = r1 then F1 = F2 = F3 and if r = r2 then F2 and F3 coincide either with F4 and F5 (if α1 > 2α2) or with F6 and F7 (otherwise). All equilibria exhibit either negative linkage disequilibrium (the three interior equilibria) or are in linkage equilibrium (the boundary equilibria). Indeed, more can be proved, namely that all orbits eventually enter the region where D ≤ 0. We consider the function Z = x2x3/x1x4 and show that Z′ > Z if D > 0 [note that D = x1x4(1 - Z)]. This follows from Z=x2x3x1x4=(x2W2+rD)(x3W3+rD)(x1W1rD)(x4W4rD)>x2x3x1x4W2W3W1W4=ZW2W3W1W4 and a straightforward calculation that reveals that W2W3W1W4.

Next, we investigate the stability properties of these equilibria.

Stability: The precise condition for local stability of the symmetric equilibrium in the case a + d = 2(b + c) was derived by Karlin and Feldman (1970). Their condition (4.6) implies that F1 is locally stable if and only if rr1. (A16)

The eigenvalues of F6 and F7 are readily calculated and are 1sα121s(α1α2)2,1sα221s(α1α2)2,1r1s(α1α2)2. Local asymptotic stability requires that these eigenvalues are less in modulus than one. Therefore, the equilibria F6 and F7 are locally stable if and only if α12α2andr>s(α1α2)2. (A17)

Also, the eigenvalues of the boundary equilibria F4 and F5 can be calculated explicitly and are given by 1s(α123α22)112s(α122α22) and 1[112s(α122α22)][112r12sα12±r28α22rs+4α12α22s2]. Again, all these eigenvalues must be less in modulus than one. The first one satisfies this condition if and only if α1 > 2α2, i.e., if and only if the equilibria exist (in S4). The expression under the square root of the second and third eigenvalue is always ≥(r - 4sα22)2 if α1 ≥ 2α2. Therefore, these are real and it follows easily that F4 and F5 are locally stable if and only if they exist (α1 > 2α2) and r13s(α12α22). (A18)

The stability analysis of the unsymmetric equilibria is rather complicated because, in general, the eigenvalues cannot be determined explicitly. However, the determinant of the Jacobian matrix (which is the constant term of the characteristic polynomial) is of relatively simple structure. It is a polynomial in r of degree six and decomposes into six linear factors. Among the zeroes of this polynomial are r1, s1 - α2)2, and 13s(α12α22) . Hence, if r = r1 or r = r2, the eigenvalues can be determined explicitly. The largest eigenvalue is one and the others are less in absolute value than one (and agree either with those of F1 or those of F4, F5 or F6, F7). Taking the derivative of the determinant of the Jacobian with respect to r and evaluating at r1 and r2 shows that all three eigenvalues are less in modulus than unity if r is slightly larger than r1 or slightly smaller than r2. Numerical calculations suggest that the unsymmetric equilibria are stable whenever they exist, i.e., if r1 < r < r2.

Thus, if 0 < rr1 then F1 is locally stable (and apparently globally stable). When r = r1, F1 loses its stability and the unsymmetric equilibria F2 and F3 arise by a pitchfork bifurcation. Then these are stable and converge to F6 and F7 as rs1 - α2)2 if α1 ≤ 2α2, while they converge to F4 and F5 as r13s(α12α22) if α1 > 2α2. In the first case the equilibria F4 and F5 do not exist. When the pair of unsymmetric equilibria hits a pair of boundary equilibria (either F4 and F5 or F6 and F7), a transcritical bifurcation occurs and the boundary equilibria become stable and remain stable for all larger values of r, while the unsymmetric equilibria leave the simplex.

Therefore, in this model no more than two equilibria can be simultaneously stable. As long as linkage is sufficiently tight, i.e., if r<r2, (A19) there exist one or two stable polymorphisms in the interior and all boundary equilibria are unstable. If r > r2 then no two-locus polymorphism can be maintained. Because r2 < ⅓ always holds (sα121 ), no stable two-locus polymorphism exists if r > 1/3. If r > r2 and the disparity between the allelic effects of the two loci is sufficiently large (α1 > 2α2), then the major locus is kept stably polymorphic by marginal overdominance while the other is monomorphic. If the disparity of effects is not large enough (α1 ≤ 2α2), no polymorphism can be maintained and either chromosome A1B2 or A2B1 becomes fixed (cf. Gavrilets and Hastings 1993).

Genetic variance: The genetic variance and all other quantities are measured before selection, so that the population is in Hardy-Weinberg equilibrium. Therefore, it is sufficient to calculate the genetic variance among gametes and multiply it by two. The mean haploid genotypic value, h, is simply Gh=α2x2+α1x3+(α1+α2)x4. (A20) From this, the deviation of the mean genotypic value at equilibrium from the optimum is calculated to be 0, R(2rs)α2 , and α1 - α2 for F1, F2 and F3, F4 and F5, and F6 and F7, respectively.

The (additive) genetic variance is given by σG2=2(α22x2+α12x3+(α1+α2)2x4Gh2). (A21) This yields the equilibrium genetic variances σ^G2(F1)=12(α12+α22)+4α1σ2D^1, (A22) σ^G2(F2,F3)=38rs[r22rs(α12+α22)+s2(a12a22)2], (A23) σ^G2(F4,F5)=12(α124α22), (A24) and σ^G2(F6,F7)=0. (A25) Although two equilibria may be simultaneously stable, it is important to note that alternative stable equilibria have the same mean phenotype and the same genetic variance.


We thank Ellen Baake and an anonymous reviewer for their comments. R.B. was supported by the Austrian Science Foundation (FWF), Project P12865-MAT, and by the Adaptive Dynamics Network of the International Institute of Applied Systems Analysis (IIASA) in Laxenburg, Austria. A.G. was supported by National Science Foundation Research Training Grant DBI9602266.


  • Communicating editor: W. Stephan

  • Received October 29, 1998.
  • Accepted March 15, 1999.


View Abstract