Balancing Selection at Closely Linked, Overdominant Loci in a Finite Population

Corresponding author: Montgomery Slatkin, Department of Integrative Biology, University of California, 3060 Valley Life Sciences Bldg., Berkeley, CA 94720-3140., slatkin{at}socrates.berkeley.edu (E-mail)

Communicating editor: N. TAKAHATA

	ABSTRACT

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

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 self-incompatibility 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 ; FELDMAN et 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) strong-selection, weak-mutation (SSWM) approximation, which has been shown to be useful for modeling balancing selection, mutation, and genetic drift at a single locus (SASAKI 1989 , SASAKI 1992 ; SLATKIN and MUIRHEAD 1999 ).

	Linkage disequilibrium in MHC

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

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 (CULLEN et 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 (FERNANDEZ-VINA et 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 (BEGOVICH et 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 well-defined 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

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

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 ² 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 1I use a threshold value of two copies, which means that two haplotypes, 0501-0302 and 0301-0303, 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_Ak_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 = , 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.

Multi-A, multi-B, 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 multi-A haplotype, indicating that one allele at A is associated with multiple alleles at B. The number of multi-A haplotypes will be denoted by I, and the fraction of multi-A 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 multi-B haplotype, indicating that one allele at B is associated with multiple alleles at A. The number of multi-B 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 ₂ = ₃ = 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 0101-0503, 0201-0201, and 0301-0301 haplotypes are not complementary. In SAUERMANN's (1998) study of rhesus monkeys, I = J = 0, and K = 14.

	DETERMINISTIC THEORY OF SELECTION

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

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 ; FELDMAN et 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₁B₁, ... , A₁B_m/2. A₂B_m/2+1, ... , A₂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₁B₁, ... , A_KB_K, are the same, z, and the frequencies of all recombinant haplotypes (A_lB_m for l m) are the same. There is always an equilibrium at which there is complete linkage equilibrium, z = , and it is easy to show that if c < /(4(K - 1)), there is a second equilibrium at which there is permanent linkage disequilibrium:

(1)

This result is similar to that of CHRISTIANSEN and FELDMAN 1975 in that the condition becomes more restrictive as the number of alleles increases. In both models, a slight change in c or can eliminate the equilibrium at which permanent linkage disequilibrium is maintained, and the solution becomes more sensitive to changes in c or as the number of alleles increases.

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 Karlin-McGregor theory can be applied to the equilibrium given by (1) to prove that it is locally stable for small values of c.

The Karlin-McGregor 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 Karlin-McGregor 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

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

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₁ appears on a chromosome initially carrying B₁, and there are m - 1 other alleles at the B locus, B₂ ... B_m, A₁ cannot later be on a chromosome carrying B₂ ... B_m. Mutation at B can then create new haplotypes carrying the mutant. For example, A₁B_m+1 could be created by mutation from B₁ to B_m+1.

At a complementary equilibrium, all haplotypes are either multi-A haplotypes, multi-B haplotypes, or unique haplotypes. As before, I is the number of A alleles in multi-A haplotypes, J is the number of B alleles in multi-B haplotypes, K is the number of unique haplotypes, _i is the fraction of multi-A haplotypes that are class i, and _j is the fraction of multi-B 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 multi-A haplotypes will have the same frequency, denoted by x_i, and all class j multi-B haplotypes will have the same frequency, denoted by y_j. These variables satisfy the condition that the haplotype frequencies sum to 1:

(2)

The numbers of alleles are

(3a)

and

(3b)

Given the selection parameters , ß, and , it is straightforward to show that

(4a)

(4b)

(4c)

Note that z is greater than both x_i and y_j and that x_i and y_j decrease with increasing values of their subscripts.

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₁B₁. When A*B₁ is in very low frequency, , it increases in frequency according to

(5)

which is positive under our restrictions on the values of and ß. If genetic drift is ignored, the A*B₁ haplotype would increase in frequency until a new equilibrium is reached at which A₁B₁ and A*B₁ are in equal frequency y₂, and the other frequencies are readjusted appropriately, thus reducing K by one and increasing J by one.

Next consider a mutation to an A that is in a class j multi-B haplotype. If one of the A's in the group A₁B₁ ... A_jB₁ mutates to A*, the frequency of A*B₁, will increase when rare according to

(6)

In this case, J will be unchanged but B₁ will move from class j to class j + 1, and the haplotype frequencies will change.

If the A that mutates is in a class i multi-A haplotype, say A₁B₁ ... A₁B_i, the frequency of the new haplotype, say A*B₁, will increase when rare according to

(7)

Once A*B₁ becomes common, then A₁B₁ will be at a disadvantage provided that

(8)

where we interpret (8) as > 0 for i = 2 and note that the right-hand side is negative for i > 2. Condition (8) is obtained by finding the rate of increase of A₁B_i from a very low frequency once A*B₁ and the A₁B₁ ... A₁B_i-1 have reached their new equilibrium frequencies. The general theory of one-locus models tells us that if a rare allele cannot increase in frequency when rare, it cannot be present at a stable polymorphic equilibrium (NAGYLAKI 1992 ). Therefore, if > 0, then after a new haplotype becomes common, selection will quickly reestablish a complementary equilibrium. At the new equilibrium, there is a new unique haplotype, A*B₁, thereby increasing K by 1, and A₁ is in class i - 1. If i = 2, a second unique haplotype, A₁B₂, remains, with the result that K is increased by 2 and I is decreased by 1.

Similar results apply to mutations at the B locus in unique haplotypes, multi-A haplotypes, and multi-B 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 multi-A and multi-B 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

(9a)

where denotes the probability of increase of I or J, the subscript A indicates that the increase is the result of a mutation at the A locus, and the subscript 1 indicates that the A that mutated was associated with only one B allele. The result of this transition is an increase in J by 1, a decrease in K by 1, an increase in ₂, and a corresponding decrease in _j for j > 2. Similarly,

(9b)

which is the probability of an increase in I by 1, a decrease in K by 1, an increase in ₂, and a corresponding decrease in _i for i > 2.

For an A allele in a class j multi-B haplotype, the net frequency of such A alleles is Jj_jy_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

(10)

In this case, there is no change in I, J, or K, but an increase in j by one. Similarly, the rate at which mutant B's in class i multi-A haplotypes increase and become common is approximately

(11)

We next consider transitions caused by mutations to new A alleles in class i multi-A haplotypes. The frequency of such haplotypes is Ii_ix_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

(12)

Similarly, the probability that a mutation of a B allele in a class j multi-B haplotype increases to become common is approximately

(13)

In this case, j would be reduced by 1 and K would be increased by 1 if j > 2, and K would be increased by 2 and J reduced by 1 if j = 2.

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 one-locus 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:

(14)

TAKAHATA's (1990) result for the average time to loss in a population of size N for the one-locus model tells us that the probability of loss per generation of a haplotype in frequency because of genetic drift is approximately

(15)

where S = 2N. As discussed by SLATKIN and MUIRHEAD 1999 , this formula differs by a factor of from that used by SASAKI 1989 , SASAKI 1992 . Although Sasaki's formula is correct for very large values of S, we found that Equation 15 is slightly more accurate when S is on the order of 100. The numerical results are nearly the same when Sasaki's formula is used instead of (15). Note that the exponential dependence on S/² implies that slight differences in result in very different rates of loss by genetic drift.

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

(16)

When such a loss occurs K is reduced by 1.

For the class i multi-A haplotypes, their total number is Ii_i, and the frequency of each is x_i. Therefore, the probability of loss per generation is

(17)

When such a loss occurs, i is reduced by 1; if i = 2, then I is reduced by 1 and K is increased by 1. Similarly, the rate of loss of class j multi-B haplotypes because of genetic drift is

(18)

and when this loss occurs, j is reduced by 1, and if j = 2, J is reduced by 1 and K is increased by 1.

Calculating _i and _j:
We can find the stationary distributions of i, _i, by noting that, once a class i = 2 multi-A 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

(19)

which is Equation 11 divided by I_i. The value of i can decrease either because of mutation of one of the A alleles or because of genetic drift. The probability of loss by mutation is _Ai/(I_i) and the probability of loss by drift is _xi/(I_i). Thus, the net probability of a decrease in i is the sum,

(20)

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, _i2, given that i = 2 initially. The stationary distribution of i, _i, is then

(21)

Although both _i and M_i are proportional to x_i² and hence depend on I, J, and K, the factor x_i² cancels from the resulting expression for _i, which then depends only on the selection and mutation parameters and the population size. In principle, there is no upper bound on i, but in the numerical analysis of this model described below, relatively small upper bounds can be assumed because _i decreases very rapidly with i. A similar analysis yields _j with appropriate redefinition of _j and M_j.

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 multi-A haplotypes,

(22a)

for multi-B haplotypes,

(22b)

and for unique haplotypes,

which, using (22a) and (22b), can be reduced to

(22c)

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. Equation 22a HREF="#FD22b">Equation 22bEquation 22cc 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 single-locus 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 multi-A or multi-B haplotypes. Furthermore, the few multi-A and multi-B haplotypes that are present will almost all be in class 2. Some typical results are shown in Table 3A–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 Equation 17 and Equation 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 multi-A and multi-B 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.

View this table: In this window In a new window   Table 3. Comparison of simulation and analytic results for > 0

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 multi-A haplotypes can be maintained at equilibrium. At this equilibrium, there are almost no multi-B haplotypes and the number of alleles at B will exceed the number at A. A typical example is shown in Table 3D.

	RECOMBINATION

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

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₁B₁, and, at the equilibrium being analyzed, A₁ is in a unique haplotype with some other B, say B₂, and B₁ is in a unique haplotype with some other A, say A₂, it is easy to show that A₁B₁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 multi-B haplotype, or when the B allele is in a unique haplotype and the A allele is in a class i multi-A haplotype. The same result also holds when the A allele of the rare haplotype is in a class i multi-A haplotype and the B allele is in a class j multi-B 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 multi-B haplotype, the rare haplotype will increase in frequency if

(23a)

Similarly, if the A allele of the rare haplotype is in a unique haplotype and the B allele is in a class i multi-A haplotype, then the rare haplotype will increase in frequency if

(23b)

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 right-hand 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 right-hand 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

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

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 Table 3 and Table 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 4aEquation 4b HREF="#FD4c">Equation 4c, which were in the range 0.01–0.04 for the parameter values used in Table 3 and Table 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.

View this table: In this window In a new window   Table 4. Simulation results for < 0

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 multi-A and multi-B 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 multi-A or multi-B haplotypes can generate noncomplementary haplotypes that will not be eliminated by selection. For example, if for a class 2 multi-A haplotype, A₁B₁ and A₁B₂, an A₁ in A₁B₁ mutates to A*, A*B₁ has a chance of increasing in frequency to become common but A₁B₁ will not be eliminated. However, A₁B₁ 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₁B₁ at the deterministic equilibrium is lower than the frequencies of either A*B₁ or A₁B₂ because it forms single homozygotes with both of them, and hence A₁B₁ 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 Karlin-McGregor 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 Table 3 and Table 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 multi-A but no multi-B or unique haplotypes. None of the cases in which > 0 and ß and are comparable in magnitude, including those shown in Table 3A–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 ₂ (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

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

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 single-locus 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 single-locus homozygotes is much less fit than the other (ß >> , but still with > 0) are the results asymmetric, with an appreciable number of multi-A or multi-B 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 single-locus 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 multi-A 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 single-locus 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 population-level 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 population-level 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 GM-40282.

Manuscript received July 23, 1999; Accepted for publication November 19, 1999.

	LITERATURE CITED

TOP ABSTRACT Linkage disequilibrium in MHC LINKAGE DISEQUILIBRIUM AT... DETERMINISTIC THEORY OF... SSWM APPROXIMATION RECOMBINATION SIMULATION RESULTS DISCUSSION AND CONCLUSIONS LITERATURE CITED

BEGOVICH, A. B., G. R. MCCLURE, V. C. SURAJ, R. C. HELMUTH, and N. FILDES et al., 1992  Polymorphism, recombination, and linkage disequilibrium within the HLA class II region. J. Immunol. 148:249-258[Abstract].

BENGTSSON, B. O. and G. THOMSON, 1981  Measuring the strength of associations between HLA antigens and diseases. Tissue Antigens 18:356-363[Medline].

CHRISTIANSEN, F. B. and M. W. FELDMAN, 1975  Selection in complex genetic systems. IV. Multiple alleles and interactions between two loci. J. Math. Biol. 2:179-294.

CULLEN, M., J. NOBLE, H. ERLICH, K. THORPE, and S. BECK et al., 1997  Characterization of recombination in the HLA class II region. Am. J. Hum. Genet. 60:397-407[Medline].

EWENS, W. J., 1979 Mathematical Population Genetics. Springer-Verlag, New York.

FELDMAN, M. W., R. C. LEWONTIN, I. R. FRANKLIN, and F. B. CHRISTIANSEN, 1975  Selection in complex genetic systems. III. An effect of allele multiplicity with two loci. Genetics 79:333-347[Abstract/Free Full Text].

FERNANDEZ-VIÑA, M. A., X. J. GAO, M. E. MORAES, J. R. MORAES, and I. SALATIEL et al., 1991  Alleles at four HLA class II loci determined by oligonucleotide hybridization and their associations in five ethnic groups. Immunogenetics 34:299-312[Medline].

GILLESPIE, J. H., 1984  Some properties of finite populations experiencing strong selection and weak mutation. Am. Nat. 121:691-708.

GRAHOVAC, B., R. I. SUKERNIK, C. O'HUIGIN, Z. ZALESKA-RUTCZYNSKA, and N. BLAGITKO et al., 1998  Polymorphism of the HLA class II loci in Siberian populations. Hum. Genet. 102:27-43[Medline].

GRIFFITHS, R. C., 1981  Neutral two-locus multiple alleles models with recombination. Theor. Popul. Biol. 19:169-186.

GYLLENSTEN, U. B. and H. A. ERLICH, 1993  MHC class II haplotypes and linkage disequilibrium in primates. Hum. Immunol. 36:1-10[Medline].

HILL, W. G. and A. ROBERTSON, 1968  Linkage disequilibrium in finite populations. Theor. Appl. Genet. 38:226-231.

KARLIN, S. and M. W. FELDMAN, 1970  Linkage and selection: two locus symmetric viability model. Theor. Popul. Biol. 1:39-71[Medline].

KARLIN, S. and J. L. MCGREGOR, 1972  Polymorphisms for genetic and ecological systems with weak coupling. Theor. Popul. Biol. 3:210-238[Medline].

KIMURA, M. and J. F. CROW, 1964  The number of alleles that can be maintained in a finite population. Genetics 49:725-738[Free Full Text].

KWOK, W. W., S. KOVATS, P. THURTLE, and G. T. NEPOM, 1993  HLA-DQ allelic polymorphisms constrain patterns of class II heterodimer formation. J. Immunol. 150:2263-2272[Abstract].

LEWONTIN, R. C., 1964  The interaction of selection and linkage. I. General considerations; heterotic models. Genetics 49:49-67[Free Full Text].

LEWONTIN, R. C. and K. KOJIMA, 1960  The evolutionary dynamics of complex polymorphisms. Evolution 14:458-472.

MARSH, S. G., 1998  HLA class II region sequences, 1998. Tissue Antigens 51:467-507[Medline].

NAGYLAKI, T., 1992 Introduction to Theoretical Population Genetics. Springer-Verlag, New York.

PETERSON, A. C., A. DI RIENZO, A. E. LEHESJOKI, A. DE LA CHAPELLE, and M. SLATKIN et al., 1995  The distribution of linkage disequilibrium over anonymous genome regions. Hum. Mol. Genet. 4:887-894[Abstract/Free Full Text].

SASAKI, A., 1989 Evolution of pathogen strategies. Ph.D. Thesis, Biology, Kyushu University, Fukuoka, Japan.

SASAKI, A., 1992 The evolution of host and pathogen genes under epidemiological interaction, pp. 247–263 in Population PaleoGenetics, edited by N. TAKAHATA. Japan Scientific Society Press, Tokyo.

SAUERMANN, U., 1998  DQ-haplotype analysis in rhesus macaques: implications for the evolution of these genes. Tissue Antigens 52:550-557[Medline].

SLATKIN, M. and C. A. MUIRHEAD, 1999  Overdominant alleles in a population of variable size. Genetics 152:775-781[Abstract/Free Full Text].

TAKAHATA, N., 1990  A simple genealogical structure of strongly balanced allelic lines and trans-species evolution of polymorphism. Proc. Natl. Acad. Sci. USA 87:2419-2423[Abstract/Free Full Text].

VEKEMANS, X. and M. SLATKIN, 1994  Gene and allelic genealogies at a gametophytic self-incompatibility locus. Genetics 137:1157-1165[Abstract].

WOLFRAM, S., 1996 The Mathematica Book. Wolfram Media Inc., Champaign, IL.

This article has been cited by other articles:

				J. F. Storz, M. Baze, J. L. Waite, F. G. Hoffmann, J. C. Opazo, and J. P. Hayes Complex Signatures of Selection and Gene Conversion in the Duplicated Globin Genes of House Mice Genetics, September 1, 2007; 177(1): 481 - 500. [Abstract] [Full Text] [PDF]

				D. Meyer, R. M. Single, S. J. Mack, H. A. Erlich, and G. Thomson Signatures of Demographic History and Natural Selection in the Human Major Histocompatibility Complex Loci Genetics, August 1, 2006; 173(4): 2121 - 2142. [Abstract] [Full Text] [PDF]

				K. L. Bubb, D. Bovee, D. Buckley, E. Haugen, M. Kibukawa, M. Paddock, A. Palmieri, S. Subramanian, Y. Zhou, R. Kaul, et al. Scan of Human Genome Reveals No New Loci Under Ancient Balancing Selection Genetics, August 1, 2006; 173(4): 2165 - 2177. [Abstract] [Full Text] [PDF]

				M. Hasselmann and M. Beye Signatures of selection among sex-determining alleles of the honey bee PNAS, April 6, 2004; 101(14): 4888 - 4893. [Abstract] [Full Text] [PDF]

				C. A. Muirhead, N. L. Glass, and M. Slatkin Multilocus Self-Recognition Systems in Fungi as a Cause of Trans-Species Polymorphism Genetics, June 1, 2002; 161(2): 633 - 641. [Abstract] [Full Text] [PDF]

				J. M. Comeron and M. Kreitman Population, Evolutionary and Genomic Consequences of Interference Selection Genetics, May 1, 2002; 161(1): 389 - 410. [Abstract] [Full Text] [PDF]