Balancing Selection at Closely Linked, Overdominant Loci in a Finite Population
 Montgomery Slatkin⇓
 Address for correspondence: Montgomery Slatkin, Department of Integrative Biology, University of California, 3060 Valley Life Sciences Bldg., Berkeley, CA 947203140. Email: slatkin{at}socrates.berkeley.edu
Abstract
High levels of allelic diversity and strong linkage disequilibrium are found in the major histocompatibility (MHC) system in humans and other vertebrates. This article proposes several descriptive statistics that quantify the extent and pattern of strong linkage disequilibrium between pairs of highly polymorphic loci. It also develops an approximate analytic theory incorporating the effects of balancing selection, mutation, recombination, and genetic drift at two closely linked loci and compares the theoretical predictions with published surveys of the MHC class II loci, DQA1 and DQB1, in humans and nonhuman primates. The descriptive statistics proposed include the fraction of complementary haplotypes (haplotypes with D″ = 1), the fraction of excess haplotypes, and the numbers of alleles at each locus in complementary haplotypes with one or more alleles at the other locus. The model assumes the infinite alleles model of mutation and the symmetric overdominance model of selection. Analytic approximations in some cases are obtained in the strong selection, weak mutation (SSWM) limit introduced by J. Gillespie. The predictions of the approximate analysis are confirmed by simulation. Both the analytic theory and simulations show that relatively few haplotypes will be found when selection is strong and recombination is weak relative to genetic drift. The model can reproduce many of the observed patterns at DQA1 and DQB1 provided that the recombination rate is assumed to be very small.
BALANCING selection resulting from either overdominance in fitness or genetic selfincompatibility maintains many more alleles at a locus than would be expected at a neutral locus with the same mutation rate. Several theoretical studies provide approximations that allow prediction of the number of alleles maintained by strong balancing selection in finite populations (Kimura and Crow 1964; Takahata 1990; Sasaki 1992; Vekemans and Slatkin 1994; Slatkin and Muirhead 1999). Balancing selection affecting two closely linked loci can maintain substantial allelic diversity at each locus and also maintain strong linkage disequilibrium between them, but the problem is much harder to analyze. There are many theoretical studies of linkage and selection affecting diallelic loci in infinite populations (Lewontin and Kojima 1960; Karlin and Feldman 1970) and a few studies of linkage and selection when one of the loci has more than two alleles (Christiansen and Feldman 1975; Feldmanet al. 1975). There are also several studies of linkage and mutation in finite populations (Hill and Robertson 1968; Griffiths 1981). But there are no analytic studies that model the combined effects of selection, mutation, and linkage in finite populations. In this article, I develop an approximate analytic theory for two closely linked overdominant loci, each of which has a potentially large number of alleles. The theory predicts the numbers of alleles maintained at each locus and the pattern of allelic association in haplotypes. The theory is based on Gillespie's (1984) strongselection, weakmutation (SSWM) approximation, which has been shown to be useful for modeling balancing selection, mutation, and genetic drift at a single locus (Sasaki 1989, 1992; Slatkin and Muirhead 1999).
Linkage disequilibrium in MHC: The present theoretical study was motivated in part by observations of linkage disequilibrium in the major histocompatibility (MHC) region in humans and other vertebrates. The MHC region contains numerous loci that are closely linked and highly polymorphic. Particularly good evidence of strong linkage disequilibrium is found in the class II loci, DQA1 and DQB1, in humans and some other primates. In humans, these two loci are ~20 kb apart on the short arm of chromosome 6. Although no recombinants between them have been reported, the recombination rate can be estimated roughly to be 0.0002 if 1 Mb is assumed to correspond to a rate of 0.01 (1 cM). Recombination rate is known to vary with location in the MHC region in humans (Cullenet al. 1997), so the recombination rate between DQA1 and DQB1 might well be <0.0002.
At DQA1, 16 alleles are known, 8 of which are in appreciable frequency worldwide; at DQB1, 25 alleles are known, 15 of which are in appreciable frequency worldwide (Marsh 1998). The two loci code for polypeptides, DQα and DQβ, that form a heterodimer. Kwok et al. (1993) showed that alleles at each locus could be grouped and that haplotypes formed from alleles in incompatible groups did not produce stable cell surface heterodimers. These “prohibited” haplotypes had not been present in population surveys (FernandezViñaet al. 1991), but recently Grahovac et al. (1998) found several of the supposedly prohibited haplotypes in their study of Siberian populations.
In humans, relatively few haplotypes are present and there is significant excess of almost all of them. Linkage disequilibrium is not complete and some alleles at DQA1 are in strong disequilibrium with two or more alleles at DQB1 (Begovichet al. 1992). Some typical data are presented in Table 1 and will be discussed in more detail later. DQA1 and DQB1 in chimpanzees and gorillas show no significant linkage disequilibrium (Gyllensten and Erlich 1993). In the rhesus monkey (Macaca mulata) DQA1 and DQB1 are in very strong linkage disequilibrium, and the pattern is different from that found in humans (Sauermann 1998). In a sample of 134 individuals, Sauermann (1998, Table 4) found 14 different haplotypes present in more than one copy. Each of the 14 haplotypes contained alleles at both loci that were in no other haplotypes in the sample. Sauermann (1998) attributed the difference between the patterns in rhesus monkeys and in humans to stronger selection acting on these two loci in the rhesus monkey.
All of these studies are based on samples of individuals collected for other purposes, so the data cannot be regarded as representing random samples from welldefined populations. The studies of Gyllensten and Erlich (1993) and Sauerman (1998) were of individuals bred in captivity, so some of the chromosomes sampled were possibly identical by descent because of past inbreeding. The nonindependence of chromosomes would affect the formal statistical analysis of the data. In this article, the data are compared only qualitatively to theoretical predictions. They are used to illustrate the application of the theory rather than to draw strong conclusions about the operation of selection in the different species discussed. Although larger samples would be desirable, they are not currently available.
LINKAGE DISEQUILIBRIUM AT MULTIALLELIC LOCI
If there are only two alleles at each of two loci, the coefficient of linkage disequilibrium, D, defined by Lewontin and Kojima (1960), is sufficient to describe the degree of nonrandom association between those loci. All other measures of nonrandom association can be expressed as functions of D and the two allele frequencies. With more than two alleles per locus, no single measure of linkage disequilibrium provides a complete description; what measure or measures are best depends on the pattern in the data and the purpose of the analysis. If the goal is to detect nonrandom association and to test its significance, then a measure based on the χ^{2} statistic, as proposed by Feldman et al. (1975) and others, or the probability value from Fisher's exact test, used by Peterson et al. (1995), will serve. If, instead, the goal is to focus on a particular allele at one locus and characterize the extent of disequilibrium with markers at another locus, then the δ statistic, suggested by Bengtsson and Thomson (1981), is more appropriate.
In other situations, there is obvious strong linkage disequilibrium involving most or all alleles at both loci and there is no reason to focus on any particular allele at either locus. As in the data presented in Table 1, very few of the many possible haplotypes are found and almost every haplotype shows significant excess or deficiency. Furthermore, many or most of the values of Lewontin's D′ (1964) statistic (defined to be the ratio of the coefficient of linkage disequilibrium to its maximum possible value) for haplotypes present in more than one copy are 1 or nearly 1. Very strong linkage disequilibrium of this kind can be generated by strong overdominant selection affecting both loci, as shown in the next several sections. To compare theory and data, some descriptive statistics are necessary to quantify this pattern.
Fraction of complementary haplotypes: I adapt the terminology of Feldman et al. (1975) and call haplotypes, which are present in more than a specified number of copies and for which D″ = 1, complementary haplotypes. Allowing for a threshold is necessary because haplotypes found in only one or a few copies may be artifacts or aberrant for some reason. For the data in Table 1, I use a threshold value of two copies, which means that two haplotypes, 05010302 and 03010303, are ignored. One descriptive statistic is the fraction of haplotypes that are complementary. For the data in Table 1, there are 14 haplotypes found in more than one copy, and 8 of them have D″ = 1, so that fraction, which is denoted by CF (for complementary fraction), is 0.571.
Fraction of extra haplotypes: Another descriptive statistic is the number of haplotypes relative to the total possible number. The minimum possible number of haplotypes is the larger of k_{A} and k_{B}, and is 11 for the data in Table 1. The maximum number is k_{A}k_{B}, which is 66 for the data in Table 1. An index describing the proportion of haplotypes greater than the minimum, which I call the fraction extra (FE) haplotypes, is FE = (k_{H} − k_{M})/(k_{A}k_{B} − k_{M}), where k_{H} is the number of haplotypes and k_{M} is the larger of k_{A} and k_{B}. This fraction is 0 if the minimum number is found and 1 if all possible haplotypes are found. For the data in Table 1, FE = 0.05.
MultiA, multiB, and unique haplotypes: It is useful to characterize the complementary haplotypes in which each allele is found. If an allele at A is in a complementary haplotype with only one allele at B, that haplotype will be called a unique haplotype, and the number of unique haplotypes will be denoted by K. If an allele at A is in complementary haplotypes with i alleles at B, that will be defined as one class i multiA haplotype, indicating that one allele at A is associated with multiple alleles at B. The number of multiA haplotypes will be denoted by I, and the fraction of multiA haplotypes that are class i will be denoted by ϕ_{i}. If an allele at B is in complementary haplotypes with j alleles at A, that will be defined to be a class j multiB haplotype, indicating that one allele at B is associated with multiple alleles at A. The number of multiB haplotypes will be denoted by J, and the fraction that are class j will be denoted by ρ_{j}. An example of the use of this notation is given in Table 2.
In Table 1, J = K = 0, I = 2, and ϕ_{2} = ϕ_{3} = 0.5. The allele 0102 at DQA1 is found in complementary haplotypes with three different alleles at DQB1, and 0103 at DQA1 is found in complementary haplotypes with two different alleles at DQB1. Note 0101, 0201, and 0301 at DQA1 do not contribute to I because the 01010503, 02010201, and 03010301 haplotypes are not complementary. In Sauermann's (1998) study of rhesus monkeys, I = J = 0, and K = 14.
DETERMINISTIC THEORY OF SELECTION
Before considering the effects of mutation and genetic drift, I review and extend slightly results from the deterministic theory of overdominant selection and recombination. The selection model in this article is the symmetric overdominance model (Christiansen and Feldman 1975; Feldmanet al. 1975), in which the relative fitness of doubly heterozygous individuals is 1, of doubly homozygous individuals is 1 − α, of individuals heterozygous only at the A locus is 1 − β, and of individuals heterozygous only at the B locus is 1 − γ. The recombination rate between the loci is c.
Two alleles per locus: This model with only two alleles at each locus has been analyzed by several authors beginning with Lewontin and Kojima (1960). We restrict our attention to the range of parameter values that seems most appropriate for the MHC region and assume α, β, and γ are all between 0 and 1 and further that α > β, γ.
The quantity ε = β + γ − α plays an important role in the deterministic theory. If ε > 0, which is sometimes called positive epistasis, the fitnesses of the single homozygotes are closer on average to the fitnesses of the double homozygotes than they are to the fitness of the double heterozygote. If ε < 0, negative epistasis, the opposite is the case. At present, there is no basis in saying whether loci in the MHC or other systems are better modeled by assuming positive or negative epistasis, but both existing deterministic theory and the new theory developed later for a finite population show that the sign of ε is important.
Lewontin and Kojima (1960) showed that, if ε > 4c, there is an equilibrium at which alleles at both loci are in frequency 1/2 and at which D is nonzero. Karlin and Feldman (1970) proved that this equilibrium is locally stable under the restricted ranges of values of α, β, and γ considered here. If ε ≤ 0, no equilibrium with D ≠ 0 exists for any c > 0.
More than two alleles per locus: Feldman et al. (1975) and Christiansen and Feldman (1975) studied this selection model when one locus has more than two alleles. These articles show that positive epistasis is necessary for the existence of stable equilibria with nonzero values of D. Christiansen and Feldman (1975, p. 195) found that, in a model with two alleles at one locus and m alleles at the other, stability of the equilibrium with only the haplotypes A_{1}B_{1},…,A_{1}B_{m/2}. A_{2}B_{m/2+1},…,A_{2}B_{m} (m even) is possible only if c<ϵ/(2m). The upper bound on c decreases inversely with the number of alleles at the B locus.
We can obtain a similar result for the generalization of the Lewontin and Kojima model to the case with K alleles at each locus. The symmetry of the problem suggests that the frequencies of all coupling haplotypes, A_{1}B_{1},…,A_{K}B_{K} are the same, z, and the frequencies of all recombinant haplotypes (A_{l}B_{m} for l ≠ m) are the same. There is always an equilibrium at which there is complete linkage equilibrium, z = 1/K^{2}, and it is easy to show that if c < ϵ/(4(K − 1)), there is a second equilibrium at which there is permanent linkage disequilibrium:
Small recombination rates: Karlin and McGregor (1972) showed that the equilibrium solutions to a model of overdominant selection for small c are very close to those for the same model of selection with c = 0, and, if an equilibrium is locally stable for c = 0, then it will also be locally stable for sufficiently small values of c. The KarlinMcGregor theory can be applied to the equilibrium given by (1) to prove that it is locally stable for small values of c.
The KarlinMcGregor theory tells us that a thorough analysis of a deterministic model with c = 0 provides a complete understanding of the behavior of a deterministic model of selection with very small values of c. We will see that the KarlinMcGregor theory is a useful starting point for analyzing overdominant selection in finite populations, but that new results are found when genetic drift and mutation are included.
SSWM APPROXIMATION
We begin with the analysis of the case in which c = 0, for which an analytic approximation can be obtained, and then use simulations to determine the effects of small amounts of recombination. With c = 0, the model reduces to a model of selection at one locus, for which many general results are known (Nagylaki 1992).
Deterministic equilibria: Initially, we consider only those equilibria that are relevant when there is positive epistasis (ε > 0). With no recombination, a new allele at one locus will remain linked to the allele at the other locus that was on the chromosome on which it first appeared, and it cannot be joined by recombination to any other allele present when it appeared. Thus, if A_{1} appears on a chromosome initially carrying B_{1}, and there are m − 1 other alleles at the B locus, B_{2} … B_{m}, A_{1} cannot later be on a chromosome carrying B_{2} … B_{m}. Mutation at B can then create new haplotypes carrying the mutant. For example, A_{1}B_{m}_{+1} could be created by mutation from B_{1} to B_{m}_{+1}.
At a complementary equilibrium, all haplotypes are either multiA haplotypes, multiB haplotypes, or unique haplotypes. As before, I is the number of A alleles in multiA haplotypes, J is the number of B alleles in multiB haplotypes, K is the number of unique haplotypes, ϕ_{i} is the fraction of multiA haplotypes that are class i, and ρ_{j} is the fraction of multiB haplotypes that are class j.
The first step in the SSWM approximation is to find the equilibrium haplotype frequencies under selection alone. By symmetry, all the unique haplotypes will have the same equilibrium frequency, denoted by z, all class i multiA haplotypes will have the same frequency, denoted by x_{i}, and all class j multiB haplotypes will have the same frequency, denoted by y_{j}. These variables satisfy the condition that the haplotype frequencies sum to 1:
Given the selection parameters α, β, and γ, it is straightforward to show that
Invasion by new mutants: If c = 0, new haplotypes are created only by mutation, and under our assumption of the infinite alleles mutation model new haplotypes will carry a new allele at one locus. We consider in turn mutations at the A locus in each of the three kinds of haplotypes. First, consider a mutant A* to an A allele in a unique haplotype, say A_{1}B_{1}. When A*B_{1} is in very low frequency, ξ, it increases in frequency according to
Next consider a mutation to an A that is in a class j multiB haplotype. If one of the A's in the group A_{1}B_{1} … A_{j}B_{1} mutates to A*, the frequency of A*B_{1}, ξ will increase when rare according to
If the A that mutates is in a class i multiA haplotype, say A_{1}B_{1} … A_{1}B_{i}, the frequency of the new haplotype, say A*B_{1}, will increase when rare according to
Similar results apply to mutations at the B locus in unique haplotypes, multiA haplotypes, and multiB haplotypes.
If ε < 0 and there is no genetic drift, invasion by a mutation will not preserve complementarity and the method described in the next section will not work. The approximate analysis can suggest what will happen when ε < 0 but at present an analytic treatment of that case does not seem possible because the solution for the deterministic equilibrium corresponding to (4) is difficult to obtain.
Finite population size: We now develop an approximation for a large finite population of constant size in which selection is strong enough that haplotype frequencies are nearly at their deterministic equilibrium values. The overall approach is to model the population as a large Markov chain with state variables I, J, K, ϕ_{i}, and ρ_{j}. We will see later that the dynamics of the multiA and multiB haplotypes can be analyzed separately, thus simplifying the analysis considerably. For reasons discussed in the previous sections, we restrict the parameters so that 1 > α >, γ > 0 and ε = β + γ − α > 0.
We assume that the population is of constant size N and that the mutation rates at the A and B loci are u and v, respectively. The probabilities that alleles in different types of haplotypes mutate are the frequencies of those haplotypes in the population, and the probability that a new haplotype increases in frequency and becomes common is approximately twice its deterministic rate of increase when rare.
Considering first the A locus, the probability that the copy of A that mutates is in a unique haplotype is Kz. The probability that the new haplotype will increase from one copy to become common is, from (5), ~2(α − β)z. Therefore, the rate per generation at which mutations of this type create haplotypes that become common is approximately
For an A allele in a class j multiB haplotype, the net frequency of such A alleles is Jjρ_{j}y_{j} and the probability of increase from a single copy is, from (6), ~2(α − β)y_{j}. Therefore, the rate at which such haplotypes arise and become common is
We next consider transitions caused by mutations to new A alleles in class i multiA haplotypes. The frequency of such haplotypes is Iiϕ_{i}x_{i}, and if a mutation occurs, the probability that the new haplotype will increase to become common is, from (7), ~2(α + (i − 1)γ − β)x_{i}. If the new haplotype becomes common, the previously common haplotype from which it arose will become disadvantageous and quickly lost provided that ε > 0. Therefore, i will be reduced by 1 and K will be increased by 1. If i = 2, then K will be increased by 2 and I will be reduced by 1. The probability of this transition per generation is
In a finite population, each haplotype is subject to stochastic loss because of genetic drift. The rate at which such loss occurs is, as in the onelocus model, the inverse of the average time to loss. For this model, the recursion equation for the deviation, f, of any haplotype from its equilibrium frequency is the same:
To compute the net rate of loss of each haplotype per generation, we multiply by the number of different haplotypes. There are K unique haplotypes and each is in frequency z, so the net rate of loss is
For the class i multiA haplotypes, their total number is Iiϕ_{i}, and the frequency of each is x_{i}. Therefore, the probability of loss per generation is
Calculating ϕ_{i} and ρ_{j}: We can find the stationary distributions of i, ϕ_{i}, by noting that, once a class i = 2 multiA haplotype is created by mutation, subsequent changes in i can be modeled by a separate Markov chain. The value of i can increase because of mutation at the B locus, and the probability of increase per generation is
Therefore, the Markov chain for changes in the value of i for i ≥ 2 has a particularly simple form for which many analytic results are known. The values of M_{i} and Λ_{i} are such that i cannot increase indefinitely, so ultimate loss from i = 2 is certain. In that case, the theory summarized by Ewens (1979, Equation 2.137 with Λ_{i} and M_{i} replacing λ_{i} and μ_{i}) provides the sojourn times, which are the average times spent in each value of i, t_{i}_{2}, given that i = 2 initially. The stationary distribution of i, ϕ_{i}, is then
Zero flux condition: Given that ϕ_{i} and ρ_{j} have reached their stationary distributions, equilibrium values of I, J, and K will be such that the fluxes into and out of the three categories of haplotypes are balanced. For multiA haplotypes,
These three equations can be solved for I, J, and K for various values of the selection parameters, mutation rates, and population size. These equations can be shown to depend on the selection and mutation parameters scaled by N: 2Nα, 2Nβ, 2Nγ, 4Nu, and 4Nv. Although they appear difficult to solve, they have a simple structure that makes numerical analysis easy. Equations 22a, 22b and 22c form a linear homogeneous set of equations for I, J, and K. The general theory of such equations tells us that there is no nontrivial solution unless the determinant of the coefficient matrix is 0. The value of z determines the values of x_{i}, y_{j}, ϕ_{i}, and ρ_{j}, so the coefficient matrix depends only on z. The problem is then reduced to finding the value of z that makes that determinant 0. Once that value of z is found, any two of the equations plus the normalization condition, (3), form a linear inhomogeneous system for I, J, and K that is trivial to solve. The numerical analysis can be done in a few seconds using a Mathematica (Wolfram 1996) program that will be provided upon request to the author.
Numerical results for c = 0: Given the five combinations of parameters (2Nα, 2Nβ, 2Nγ, 4Nu, and 4Nv), the approximate equations described in the previous section yield I, J, and K, as well as the x_{i}, y_{j}, ϕ_{i}, and ρ_{j}. Although there are many possible combinations of parameters, and many patterns in the results can be envisioned, the model's predictions are relatively simple: one of two patterns is found. If the singlelocus homozygotes have roughly the same relative fitnesses, meaning that β and γ have approximately the same values, then almost all haplotypes will be unique and very few will be multiA or multiB haplotypes. Furthermore, the few multiA and multiB haplotypes that are present will almost all be in class 2. Some typical results are shown in Table 3, A–C. That pattern is found even if there are substantial differences in the mutation rates. We can see why by recalling the very strong dependence of rate of haplotype loss by genetic drift, ~_{xi} and ~_{yj} in Equations 17 and 18, on the haplotype frequencies. If ε > 0 and β and γ are comparable in magnitude, then (4b) and (4c) imply that x_{i} and y_{j} are less than z and decrease with i and j. That ensures that multiA and multiB haplotypes will be lost to drift at a much higher rate than will unique haplotypes, especially if i or j exceeds 2. If the mutation rates differ, that causes only a linear difference in the λ's, while there is highly nonlinear dependence of the ~'s on x_{i} and y_{j}. As a consequence, the numbers of alleles at each locus are roughly equal. This tendency does not depend on exact equality of β and γ, but only on the fact that they are of the same order of magnitude.
This explanation for the pattern found when β and γ are comparable in magnitude suggests when another pattern may be found. If γ is quite small, then (4) implies that x_{i} is nearly equal to z and decreases only slightly with increasing i. When γ is small, the condition ε > 0 requires that β then be nearly equal to α. In this case, the symmetry disappears and a substantial number of multiA haplotypes can be maintained at equilibrium. At this equilibrium, there are almost no multiB haplotypes and the number of alleles at B will exceed the number at A. A typical example is shown in Table 3D.
RECOMBINATION
We next consider small amounts of recombination. The Karlin and McGregor (1972) theory indicates that the stability of the equilibrium for c = 0 is important.
Stability of deterministic equilibria: The complementary equilibrium defined by (4) is feasible, meaning that all the haplotype frequencies are between 0 and 1 for all values of α, β, and γ between 0 and 1. But it is not necessarily stable to invasion by noncomplementary haplotypes created by recombination.
The local stability of (4) can be determined by finding whether haplotypes of different kinds will increase in frequency when rare. If even one noncomplementary haplotype can increase, the equilibrium is unstable. There are nine kinds of noncomplementary haplotypes that can be created by recombination, depending on the common haplotypes the A and B alleles are in, but the analysis of each case is not difficult. To illustrate, let ξ be the frequency of a rare haplotype containing alleles at both loci that are otherwise found in unique haplotypes. For example, if the rare haplotype is A_{1}B_{1}, and, at the equilibrium being analyzed, A_{1} is in a unique haplotype with some other B, say B_{2}, and B_{1} is in a unique haplotype with some other A, say A_{2}, it is easy to show that A_{1}B_{1}will decrease in frequency if ε > 0 and will increase in frequency if the inequality is reversed.
In a similar way, the condition ε > 0 ensures that the recombinant haplotype will decrease in frequency when the A allele is otherwise in a unique haplotype and the B allele is otherwise in a class j multiB haplotype, or when the B allele is in a unique haplotype and the A allele is in a class i multiA haplotype. The same result also holds when the A allele of the rare haplotype is in a class i multiA haplotype and the B allele is in a class j multiB haplotype.
However, local stability is not guaranteed by ε > 0 for the other noncomplementary haplotypes. For example, if the B allele of the rare haplotype is in a unique haplotype and the A allele is in a class j multiB haplotype, the rare haplotype will increase in frequency if
These two conditions imply that ε > 0 is not sufficient for (4) to be locally stable. For (23a) and (23b) to be satisfied for i ≥ 2 and j ≥ 2, both β and γ must be close to α. For example, if α = 0.05, β = γ = 0.03 (Table 3A), then the righthand sides of (23a and b) are both 1.83. If, instead, β = γ = 0.045 (Table 3C), the critical values of i and j increase to 3.08. In general, larger values of ε permit local stability of (4) for larger values of i and j, but only with very large values of α are much larger values of i and j possible in a stable equilibrium. If β is close to α but γ is quite small, the range of values of i and j for which (4) is stable is also restricted. For example, with α = 0.05, β = 0.0475, and γ = 0.005 (Table 3D), the righthand sides of (23a) and (23b) are 1.56 and 2.05, respectively.
These results suggest that the complementary equilibria predicted when c = 0 are not of biological interest because they are not locally stable. Simulation results presented below show otherwise. With very low recombination rates, recombinant haplotypes that could destabilize the equilibrium appear only very infrequently, and their rate of increase when rare is very small, meaning that they would have a very small chance of becoming common. As a consequence, complementary equilibria will be found in finite populations even if they would be unstable in infinitely large populations.
SIMULATION RESULTS
The predictions of the approximate theory were compared with those from a computer simulation of the same model. The computer model used a rejection scheme to simulate overdominant selection. In all cases shown in the tables, a population size of N = 5000 was used. For each set of parameter values, samples from the population were taken every 10,000 generations, beginning in generation 30,000 and continuing to generation 200,000. Results presented in the tables are values averaged over those 18 samples. Other simulations showed that this sampling scheme was sufficient to obtain results close to those from simulations carried out for much longer times.
In the calculation of the various statistics shown, it was necessary to choose a minimum number of copies of a haplotype counted as a common haplotype. The results shown in Tables 3 and 4 are for the entire population. Haplotypes were counted in the analysis only if >100 copies were present in the population. That number is less than the number of haplotypes expected at the deterministic equilibrium frequencies Equation 4, which were in the range 0.01–0.04 for the parameter values used in Tables 3 and 4. As discussed below, the choice of that minimum number affects the results somewhat but does not change the overall patterns found. In Table 5, samples of 253 chromosomes were chosen at random without replacement from each data set for comparison with the data shown in Table 1, and in those cases a minimum value of 2 was used, as in the analysis of the data in Table 1.
Positive epistasis: If ε > 0, the simulation results can be compared to predictions from the analytic theory. Some results are shown in Table 3. For c = 0, the simulation results agree roughly with the analytic predictions. When the threshold haplotype number is 100, there tend to be slightly more multiA and multiB haplotypes than predicted. If that threshold is raised to 200, the simulation results are slightly different. For example, in Table 3A, if a threshold value of 200 copies is used, the averages of I, J, and K in the simulations for c = 0 are 0.89, 0.39, and 8.2, respectively.
When c is of the same order of magnitude as u and v, the simulation results are still similar to analytic predictions. That is consistent with the Karlin and McGregor (1972) theory. When c is an order of magnitude larger than u and v, the simulation results no longer agree with the analytic theory. Fewer complementary haplotypes remain. Nevertheless, there is still substantial linkage disequilibrium present, many of the possible haplotypes are still not found in appreciable frequencies, and many complementary haplotypes are still present.
Negative epistasis: If ε < 0, the analytic theory no longer provides an approximate solution but it does suggest what will be seen. Even without recombination, some mutations in multiA or multiB haplotypes can generate noncomplementary haplotypes that will not be eliminated by selection. For example, if for a class 2 multiA haplotype, A_{1}B_{1} and A_{1}B_{2}, an A_{1} in A_{1}B_{1} mutates to A*, A*B_{1} has a chance of increasing in frequency to become common but A_{1}B_{1} will not be eliminated. However, A_{1}B_{1} is still liable to be lost because of drift, so complementarity can be restored. Although the solution for the deterministic equilibrium in that case cannot be obtained in a simple form, it is easy to see that the frequency of A_{1}B_{1} at the deterministic equilibrium is lower than the frequencies of either A*B_{1} or A_{1}B_{2} because it forms single homozygotes with both of them, and hence A_{1}B_{1} will have a higher rate of loss by drift. The magnitude of this effect cannot be predicted but we can expect that substantial complementarity can be found if ε < 0 and that a higher degree of complementarity will be found if ε is negative but small in absolute value. These patterns are seen in Table 4 for c = 0.
For cases with ε < 0 and c > 0, the Karlin and McGregor (1972) theory is not useful for predicting patterns found in the simulations. If ε < 0, a completely complementary deterministic equilibrium with linkage disequilibrium is unstable to invasion by recombinant haplotypes, so the KarlinMcGregor theory predicts that such an equilibrium should be structurally unstable to increases in c from 0. That is not seen in the simulation results. Extensive linkage disequilibrium is found for small values of c even when ε is relatively large in absolute value (ε = −0.04 in Table 4A). The loss of complementarity with increasing c is more rapid than in the cases with ε > 0, because selection is no longer acting to maintain it. But the observation of substantial linkage disequilibrium and relatively high complementarity does not imply that ε is positive.
Small sample sizes: The simulation results shown in Tables 3 and 4 were chosen to illustrate the patterns found in the whole population. I do not consider the general sampling problem but present results for samples of the same size as in Table 1, to determine what kinds of selection are consistent with these observations. We can restrict our attention to cases that exhibit some of the general features of the data in Table 1, particularly the difference in the numbers of alleles at the two loci, the relatively small number of haplotypes found, the high degree of complementarity, and the presence of multiA but no multiB or unique haplotypes. None of the cases in which ε > 0 and β and γ are comparable in magnitude, including those shown in Table 3, A–C, predict those patterns. The cases shown in Table 3D and in Table 4 are similar enough to provide an adequate comparison with Table 1.
Table 5 shows that many of the features found in the data in Table 1 can be reproduced by this model. Values of c of the same order of magnitude as the mutation rates give the closest agreement. If ε is positive (Table 5A) or only slightly negative (Table 5C), strong asymmetry in selection on the single homozygotes (β ⪢ γ) is required, but if ε is much <0, then a modest asymmetry in mutation rates (Table 5B, with v = 4u) produces similar results even with β = γ. In all these cases, we see differences in the numbers of alleles maintained, 0 or nearly 0 values of J, relatively small values of ϕ_{2} (indicating relatively large values of ϕ_{i} for i > 2), and values of both CF and FE that are roughly comparable to those found in the data. The simulation results differ from Table 1 by predicting nearly the same values of I and K, while in Table 1, I = 2 and K = 0.
For the data of Sauermann (1998) for the rhesus monkey, the pattern is very different. All 14 haplotypes are unique haplotypes, so K = 14, I = J = 0, and there is perfect complementarity (CF = 1) and no extra haplotypes (FE = 0). That pattern is found if ε is positive and relatively large, which means that β and γ are comparable in magnitude and each nearly as large as α, and if c is very small.
DISCUSSION AND CONCLUSIONS
This article has three goals. The first is to present new ways of quantifying the pattern of linkage disequilibrium between two multiallelic loci when there is obvious significant nonrandom association between them. The second is to develop an approximate analytic model that accounts for selection and mutation at two multiallelic loci in a finite population. The third is to determine the extent to which a simple highly symmetric model of selection can account for patterns found in data from the DQA1 and DQB1 loci in the MHC regions of humans and nonhuman primates.
Quantifying the extent of disequilibrium when there is strong balancing selection requires new approaches because it is not enough to show that there is significant disequilibrium. In many cases, linkage disequilibrium is so apparent that statistical tests are unnecessary. The problem instead is to distinguish between the many kinds of patterns that could be found. The approach taken here is based on the theoretical result that with very low recombination rates, selection, mutation, and drift will together create and maintain complementarity, with relatively few haplotypes found in appreciable numbers and many or most of them having Lewontin's D″ = 1. That theoretical result suggests that FE and CF indicate how close to complete complementarity a particular data set is. If high complementarity is found, then the numbers of alleles at one locus in haplotypes with each allele at the other locus, measured by I, J, K, ϕ_{i}, and ρ_{j}, are of interest because their values are predicted by the analytic model.
The analytic model predicts the simulation results with reasonable accuracy when there is positive epistasis (ε > 0) and low rates of recombination. When there is positive epistasis, one of two patterns is found. If the fitnesses of the singlelocus homozygotes are comparable (β and γ have similar magnitudes), then at an equilibrium, most haplotypes are unique and roughly equal numbers of alleles are found at both loci even if mutation rates differ. Only if one of the singlelocus homozygotes is much less fit than the other (β ⪢ γ, but still with ε > 0) are the results asymmetric, with an appreciable number of multiA or multiB haplotypes and a substantial difference in the numbers of alleles at the two loci.
If ε < 0, the analytic model can no longer be used but the development of that model suggests that nearly complete complementarity will be found for small enough recombination rates. Simulation results confirmed that prediction, but with two differences from the cases with ε > 0. Recombination is more effective in destroying complementarity and differences in the mutation rates can more easily lead to asymmetry of the results.
Both the analytic and simulation models are not intended to provide a complete description of DQA1 and DQB1 or any other pair of loci. The model was designed to have as few parameters as possible. Selection is highly symmetric, with the many possible genotypes grouped into only four fitness classes. There are no specific interactions among alleles, as has been demonstrated for DQA1 and DQB1 by Kwok et al. (1993), no allowance for the possibility that some alleles may carry intrinsic disadvantages, as modeled by Slatkin and Muirhead (1999), and no possibility that mutation rates to different alleles or classes of alleles at each locus might differ. The goal in comparing data to theory is to see how good the predictions of this simple model are and where apparent differences lie.
For Sauermann's (1998) data for the rhesus monkey, the theory fits the data well provided that selection is strong, the fitnesses of the singlelocus homozygotes are comparable and nearly as small as the fitnesses of the double homozygotes, and the recombination rate is low, not much larger than the mutation rates to functionally different alleles at each locus. The agreement between theory and data does not mean that the model's assumptions are valid for DQA1 and DQB1 in rhesus monkeys, only that no additional assumptions are necessary to account for the observations.
The fit of the model's predictions to the data in Table 1 is less good. Many of the features can be reproduced by the model but it does not seem possible to predict more multiA haplotypes than unique haplotypes. That difference would probably require a more complex model of selection that embodied the functional restrictions of the type found by Kwok et al. (1993). Still, the parameter values needed to achieve the degree of asymmetry found in the data all have in common that one or both of the singlelocus homozygotes must have a fitness comparable to the double homozygote. Not enough is known about the effect of these two loci on fitness to know whether that is true for DQA1 and DQB1 in humans, but at least the model does suggest that there are differences in the fitness interactions between these loci in humans and in the rhesus monkey.
For chimpanzees and gorillas, Gyllensten and Erlich (1993) found no significant linkage disequilibrium. They suggested that either the recombination rate is much higher or selection is weaker than in humans. The theory described here supports that view but shows that differences in either of these parameters would have to be by an order of magnitude or more. Small differences from the values in humans would not be sufficient to eliminate linkage disequilibrium. A difference in effective population size could not account for the patterns observed, because both the recombination rates and selection coefficients are scaled by the same factor, so the relative intensities of selection and recombination would remain the same.
These conclusions about the action of selection at DQA1 and DQB1 have to be regarded as tentative because there are too few populationlevel studies of haplotype frequencies in populations of humans and nonhuman primates to know whether those patterns are found in other populations as well. The data are discussed to illustrate the application of the theory and to show that populationlevel data can be used to test hypotheses about selection.
Acknowledgments
I thank S. McWeeney for providing information and references about linkage disequilibrium in DQA and DQB and M. W. Feldman, S. McWeeney, and G. Thomson for helpful discussions and advice. This research was supported in part by National Institutes of Health grant GM40282.
Footnotes

Communicating editor: N. Takahata
 Received July 23, 1999.
 Accepted November 19, 1999.
 Copyright © 2000 by the Genetics Society of America