THEORY

Consider two linked loci, I and II, in a random mating population with N diploids. We consider two neutral alleles, A and a, so that there are four haplotypes, A-A, A-a, a-A, and a-a (the first letter represents the allele at locus I and the second one represents the allele at locus II). It is assumed that the mutation rate between two alleles is μ per locus per generation. The recombination rate between two loci is assumed to be r per generation. Intrachromosomal gene conversion occurs at the rate c per locus per generation; e.g., A-a changes into A-A with probability c and into a-a with the same probability. Interchromosomal gene conversion is not considered. Let the frequencies of A-A, A-a, a-A, and a-a be x₁, x₂, x₃, and x₄ (x₁ + x₂ + x₃ + x₄ = 1), respectively. Given x₁, x₂, x₃, and x₄, their expectations in the next generation are given by x1′=(1−2μ)x1+(μ+c)(x2+x3)−rD, (1a) x2′=(1−2μ−2c)x2+μ(x1+x4)+rD, (1b) x3′=(1−2μ−2c)x3+μ(x1+x4)+rD, (1c) and x4′=(1−2μ)x4+(μ+c)(x2+x3)−rD, (1d) where D = x₁x₄ – x₂x₃.

Under this model, we calculate the expectations of moments of allele frequencies using a diffusion method, which was introduced to population genetics by Kimura (1964). In equilibrium, it is known that a function, g(x₁, x₂, x₃), satisfies the equation E[L(g)]=0, (2) where L is the differential operator of the Kolmogorov backward equation (Kimura 1964; Ohta and Kimura 1969a). In this model, L(g) is given by L(g)=x1(1−x1)4N∂2g∂x12+x2(1−x2)4N∂2g∂x22=x3(1−x3)4N∂2g∂x32−x1x22N∂2g∂x1∂x2−x1x32N∂2g∂x1∂x3−x2x32N∂2g∂x2∂x3+[−2μx1+(μ+c)(x2+x3)−rD]∂g∂x1+[−2(μ+c)x2+μ(1−x2−x3)+rD]∂g∂x2+[−2(μ+c)x3+μ(1−x2−x3)+rD]∂g∂x3. (3) We can transform the three variables, x₁, x₂, and x₃, in Equation 3 into p, q, and D: p=x1+x2,q=x1+x3,D=x1x4−x2x3=x1−x12−x1x2−x1x3−x2x3. (4) p and q represent the frequencies of A in loci I and II, respectively. Then, (3) becomes L(g)=L′(g)4N, (5) and L′(g)=p(1−p)∂2g∂p2+q(1−q)∂2g∂q2+[pq(1−p)(1−q)+D(1−2p)(1−2q)−D2]∂2g∂D2+2D∂2g∂p∂q+2D(1−2p)∂2g∂p∂D+2D(1−2q)∂2g∂q∂D+[θ(1−2p)−C(p−q)]∂g∂p+[θ(1−2q)+C(p−q)]∂g∂q+[Cp(1−p)+Cq(1−q)−(2+4θ+2C+R)D]∂g∂D, (6) where θ = 4Nμ, C = 4Nc, and R = 4Nr. Without gene conversion (C = 0), this equation is the same as Equation 12 in Ohta and Kimura (1969b).

First, letting g = p and q in (2) and (6), we have E(p)=E(q)=0.5, (7) when θ ≠ 0 and C ≠ 0.

Next, letting g = p², q², pq, and D, we obtain the following four equations: 1+θ−2(1+2θ+C)E(p2)+2CE(pq)=0, (8) E(p2)=E(q2), (9) 2E(D)+CE(p2)+CE(q2)−(4θ+2C)E(pq)+θ=0, (10) and C−CE(p2)−CE(q2)−(2+4θ+2C+R)E(D)=0. (11) From (8, 9, 10, 11), we have E(p2)=E(q2)=λω, (12) E(pq)=−1+θ2C+(1+α)λCω, (13) E(D)=Cβ(1−2λω), (14) where α=2θ+C, (15a) β=2+2α+R, (15b) λ=4C2+β[2θC+2α(1+θ)], (15c) and ω=8C2+4β[α(1+α)−C2]. (15d)

Therefore, the expectations of the amounts of variation (heterozygosity) within loci I and II are given by E(hwI)=E[2p(1−p)]=1−2E(p2) (16a) and E(hwII)=E[2q(1−q)]=E(hwI), (16b) respectively, where E(p²) is given by (12). The expected amount of variation within a locus is given by the average of h_wI and h_wII: E(hw)=E(hwI)=E(hwII). (16c) Define the amount of variation between two loci, h_b, as the probability that a pair of alleles randomly chosen from each locus are different. The expectation of h_b becomes E(hb)=E[p(1−q)+(1−p)q]=1−2E(pq), (17) where E(pq) is given by (13).

In a similar way, the variances of h_wI, h_wII, and h_b are written as Var(hwI)=Var(hwII)=E(hwI)2−[E(hwI)]2=4E(p4)−8E(p3)+4E(p2)−[E(hw)]2 (18) and Var(hb)=E(hb2)−[E(hb)]2=4E(p2q2)−8E(p2q)+2E(p2)+2E(pq)−[E(hb)]2. (19) The derivations for E(p⁴), E(p³), E(p²q²), and E(p²q) are shown in the appendix. We can also obtain the covariance between h_wI and h_wII and the variance of D. That is, Cov(hwI,hwII)=4E(p2q2)−8E(p2q)+4E(pq)−E(hwI)E(hwII) (20) and Var(D)=E(D2)−[E(D)]2, (21) where the derivation for E(D²) is also in the appendix. From (18) and (20), the variance of h_w is given by Var(hw)=Var(hwI)∕2+Cov(hwI,hwII)∕2. (22)

Numerical examples for E(h_w), E(h_b), and E(D) are shown in Figure 1. Figure 1A shows the results for E(h_w) given θ = 0.01. Gene conversion increases the amount of variation within a locus. Note that E(h_w) = 0.0098 without gene conversion. When the gene conversion rate is relatively small (C = 0.1), E(h_w) is ∼1.75-fold larger than that without gene conversion, while there is almost no effect of gene conversion on E(h_w) when C = 100. Recombination also increases E(h_w) but the effect is relatively small. Figure 1B shows the results for E(h_b). Gene conversion decreases the amount of variation between two loci. The amount of variation between two loci is much bigger than that within each locus unless C is very large. When C = 100, E(h_w) and E(h_b) are almost the same. In Figure 1C, it is shown that gene conversion generates positive linkage disequilibrium. When there is no gene conversion, E(D) = 0. D is positively correlated with C, and D decreases as R increases. These results are consistent with other studies (e.g., Ohta 1982).

DATA ANALYSIS AND ESTIMATION OF PARAMETERS

Since the expectations and variances of the amounts of variation within and between two loci and linkage disequilibrium are given by functions of θ, C, and R, it may be possible to estimate these parameters from DNA polymorphism data. An estimation method is explained using the data of the Amy region in D. melanogaster as an example (see Figure 1 in Inomataet al. 1995). The data consist of the coding sequences of the proximal and distal Amy genes for nine strains. The length of coding sequence is 1482 bp. In the alignment of 18 sequences (9 sequences from the proximal gene and 9 from the distal gene), 47 sites are polymorphic, of which 37 are synonymous.

First, we estimate the amounts of variation and linkage disequilibrium for a particular site. Consider the 567th site of the Amy genes where two nucleotides, T and C, are segregating, so that there are four possible haplotypes, T-T, T-C, C-T, and C-C (the first letter represents the nucleotide in the proximal gene and the second one represents that in the distal gene). Denote the number of these haplotypes by n₁, n₂, n₃, and n₄. Estimates of heterozygosity within the proximal and distal genes are given by hwp=(n1+n2)(n3+n4)n(n−1)∕2andhwd=(n1+n3)(n2+n4)n(n−1)∕2, (23a) where n = n₁ + n₂ + n₃ + n₄. Then, we have the average of h_wp and h_wd as hw=(hwp+hwd)∕2. (23b) The amount of variation between two genes is estimated by hb=(n1+n2)(n2+n4)+(n1+n3)(n3+n4)−n2−n3n(n−1). (24) Since n₁ = 1, n₂ = 5, n₃ = 0, and n₄ = 3 at the 567th site, h_wp = 0.5, h_wd = 0.222, h_w = 0.361, and h_b = 0.639. Next, we consider the amount of linkage disequilibrium. Usually linkage disequilibrium in the sample is calculated as (n₁n₄ – n₂n₃)/n². Therefore, from Nei and Roychoudhury (1974), an estimate of linkage disequilibrium in the population may be given by δ=n1n4−n2n3n(n−1), (25) from which we have δ= 0.0417 at the 567th site.

From (23, 24, 25), we can calculate h_w, h_b, and δ for all sites of the genes and we have their averages. The averages of h_wp and h_wd correspond to π_wp and π_wd, the average numbers of pairwise differences within the proximal and distal genes per site. π_b is the average number of pairwise differences between two genes, which is the average of h_b. Let d be the average of δ. Only the data for the synonymous sites of Inomata et al. (1995) are used for the calculation because their sampling was not random. The sampling was based on the information of allozyme variation (see Inomataet al. 1995). From all synonymous sites, we have π_wp = 0.0315 and π_wd = 0.0302. Then, the average number of pairwise differences within a gene, π_w, is 0.0309. In a similar way, we have the average number of pairwise differences between two genes, π_b = 0.0452. The average of linkage disequilibrium between two genes, d, becomes 0.000452. If θ, C, and R are constant for all the sites, the expectations of π_w, π_b, and d are given by E(πw)=E(hw),E(πb)=E(hb),E(d)=E(D), (26) where E(h_w), E(h_b), and E(D) are given by (16c), (17), and (14), respectively.

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

—E(h_w), E(h_b), and E(D) given θ = 0.01. (A) Results for E(h_w). (B) Results for E(h_b). (C) Results for E(D).

Since E(π_w), E(π_b), and E(d) are given by functions of θ, C, and R, it may be possible to estimate these parameters from π_w, π_b, and d, although the equations for E(π_w), E(π_b), and E(d) are too complicated to solve for θ, C, and R. One way for the estimation is to find a set of θ, C, and R that minimizes x: x=[πw−E(πw)]2Var(πw)+[πb−E(πb)]2Var(πb)+[d−E(d)]2Var(d). (27) Although we do not have analytical expressions for Var(π_w), Var(π_b), and Var(d), we may be able to use (22), (19), and (21) for them, respectively, because these equations are used as weighting factors in (27). Note that the equations for Var(h_w), Var(h_b), and Var(D) are based on the two-locus model. The variances of π_w, π_b, and d could be smaller than those of h_w, h_b, and D, as π_w, π_b, and d are calculated from DNA sequence data. For the Amy genes of D. melanogaster, given π_w = 0.0309, π_b = 0.452, and d = 0.000452, the minimum x is obtained when θ = 0.0172, C = 1.03, and R = 66.6. Unfortunately, it is not possible to evaluate the variances of these estimates. They might depend on θ, C, R, and the sample size (n). Recombination within each gene may decrease the variances.