TWO APPROACHES TO ESTIMATE h¯ AND σh2

Hypothesis testing and parameter estimation in statistics are two highly related topics. Unbiased and efficient estimation (estimation with a small sampling error) of a parameter generally forms a basis for a powerful test concerning that parameter. For a locus with the two alleles A and a, let the three genotypic values be, respectively: AAAaaa11−hisi1−si Then h_i < 0 implies overdominance, h_i = 0.5 additivity, 0 ≤ h_i ≤ 1 (h_i ≠ 0.5) dominance, and h_i > 1 implies underdominance. Note that throughout we use “dominant” or “dominance” to refer to the cases of 0 ≤ h_i ≤ 1 (h_i ≠ 0.5), which include both complete dominant or recessive, and partial dominant or recessive. For individual locus of quantitative traits, h_i is almost impossible to estimate, and thus h¯ is normally estimated. Even for h¯ (the arithmetic mean of h_i across loci), there currently is no method to directly estimate it. Often, h¯ is approximately estimated by Mukai's method (Mukaiet al. 1972), which is closely related to the earlier methods of Comstock and Robinson (1949) and Hayman (1954). Because the arithmetic mean is the most often used measure of the average, the investigation of the bias of Mukai's and our new method will be relative to the true h¯ . Throughout, h_i will refer to dominance coefficient at individual loci, h¯ refers to the arithmetic mean of h_i across loci, and h refers to the general term dominance coefficient or the dominance coefficient with constant h_i across loci.

Mukai's approach: This approach was developed under the assumptions that dominance is the sole mode of within-locus genetic effects, the frequency of deleterious allele is very small, the population is at Hardy-Weinberg equilibrium, and mutation effects across loci are additive. It approximately estimates h¯ by the slope of the regression of the outcrossed-progeny fitness (x) on the fitness sum (y) of the two corresponding parental homozygotes (Mukaiet al. 1972; Simmons and Crow 1977; Crow and Simmons 1983) h¯≈Cov(x,y)Var(y). (1a) σh2 can be approximately estimated as (Mukai and Yamaguchi 1974) σh2≈Var(x)Var(y)−(Cov(x,y)Var(y))2. (1b) This approach is readily applicable to highly selfing populations (Johnston and Schoen 1995), where construction of homozygotes is easy. In most outcrossing populations, however, construction of pure homozygotes is usually difficult, except in a few organisms (such as Drosophila) where it can be relatively easily achieved by special chromosomal constructs (e.g., Mukaiet al. 1972; Eaneset al. 1985). Therefore, the following approach is developed to estimate h¯ and σh2 in self-compatible outcrossing populations without homozygous lines.

New approach: A wide variety of outcrossing plants and invertebrates are capable of selfing. In such populations, if we self a random sample of genotypes and obtain a number of selfed progeny from each parent to form selfed families, then h¯ and σh2 can be estimated. The data needed are the genotypic value of the parent (w) and the mean genotypic value (z) of the selfed progeny within each selfed family. Under the same assumptions as Mukai's approach outlined above, for a diallelic locus, we define the one-locus genetic effects as in Table 1. It can be easily seen from Table 1 that, The variance oft′:Var(t′)≅2pqsi2 (2a) The Variance ofw′:Var(w′)≅2pqsi2hi2 (2b) The Covariance ofw′andt′:Cov(w′,t′)≅2pqsi2hi (2c) Then the expected regression of the mutation effects of outcrossed parent (w′) on the quantity [4*(mean mutation effects of selfed family z′) − 2*w′] is given by the ratio of the covariance to the variance, summed over all the relevant loci. This gives an approximate estimate of h¯ Cov(w′,t′)Var(t′)=Σpiqisi2hiΣpiqisi2≈h¯. (3a) As with Mukai's approach, h¯ is approximately estimated by the average of h_is at individual loci weighted by the genetic variance of the homozygotes. Although the above derivation is based on the mutational effects for ease of derivation, it can be easily shown (because w = 1 − w′ and z = 1 − z′) that the estimation can be performed on the original trait values Cov(w,t)Var(t)=Σpiqisi2hiΣpiqisi2≈h¯ (3b) where t = 4z − 2w. Under constant mutation effects, Cov(w,t)Var(t)=h. (3c) Additionally, we note that the ratio of the variance of w′ to that of t′ is Var(w′)Var(t′)=Var(w)Var(t)=Σpiqisi2hi2Σpiqisi2≈h2¯. (4a) As with the estimates of h, h¯,h2¯ is approximately estimated by the average of hi2 at individual loci weighted by the genetic variance of the homozygotes. Thus, an approximate σh2 estimate is σh2=h2¯−h¯2≈Var(w)Var(t)−(Cov(w,t)Var(t))2. (4b) In the above derivation, note that the assumption of Hardy-Weinberg equilibrium is not necessary; in fact, neither is it in Mukai's approach. One key assumption, though, is that the frequency of the rarer allele at any locus is low. This is essential in order for the approximation of Equation 2, a, b and c, to hold reasonably well, and also is true in deriving Equation 1, a and b (Mukaiet al. 1972).

SIMULATIONS

The above derivations make a number of assumptions, for example, the within-locus nonadditive genetic effects are dominant and genetic effects across loci are additive. However, there is good evidence that genes for fitness or its components usually act multiplicatively (Mortonet al. 1956; Crow 1986; Fu and Ritland 1996). To investigate the robustness of the above approaches under overdominance and the more reasonable mixed modes of dominance and overdominance at different loci, and to investigate their statistical properties, simulations are performed in which fitnesses are assumed to be multiplicative. To concentrate on studying the robustness of the approaches and the influence of the genetic process (selfing outcrossed individuals or outcrossing homozygous lines) on the estimation, all genotypicvalues are assumed to be measured accurately. In reality, this would require that each genotype be clonally replicated and assayed a very large number of times. Ignoring measurement error of genotypic values due to random environmental and developmental processes will likely inflate the sampling error of estimation, but unlikely bias the estimation and comparison of the two approaches under the same sample sizes of genotypes (Deng and Lynch 1996a).

This section is organized according to the presentation of different mutation effects across the genome, starting from the simplest case of constant effects, progressing to biologically more complex and plausible situations. Simulations will be described for each situation, respectively, and some necessary analytical results will be developed.

Constant mutation effects with dominance: Dominance (h_i) and selection (s_i) coefficients across loci are the same, that is, h_i = h, s_i = s.

Mukai's approach in outcrossing populations: Assume some random pairs of homozygotes are established from natural populations [such as with a special chromosome construct in Drosophila (Mukaiet al. 1972)]. Their fitness is W = (1 − s)ⁿ, where n is the number of mutations randomly determined from a Poisson distribution with mean U/(2hs), where U is the genomic mutation rate. This is because the mean number of mutations per genome in an outcrossing population is U/(hs), and nearly all mutations are heterozygous in state (Deng and Lynch 1996a), and homozygous lines established are expected to carry only about one-half of the total genomic mutations of the outcrossed generation. Because genome size is usually very big and the mutation rate (μ) of each locus is very small, the mutant allele frequency (μ/hs at mutation-selection balance; Crow and Kimura 1970) is also very small. It is unlikely that the homozygotes established have mutations at the same loci (Charlesworthet al. 1990; Deng and Lynch 1996a), as corroborated by our computer simulations (H.-W. Deng, unpublished data). Therefore, throughout for dominant loci under mutation-selection equilibrium, the number of mutations in an outcrossed progeny is the sum of the number of mutations of its two parental homozygotes (n_m and n_f), but all at heterozygous state. Thus, the fitness of an outcrossed progeny is W = (1 − hs)^{n_m + n_f}. Because lnW=ln(1−s)n≈−ns,(s<<1.0) (5) ln(W) therefore approximates the fitness reduction under the additive fitness function. Thus, throughout, logarithmic transformation is performed for fitness, and then Equation 1a is employed to estimate h. Throughout, in simulating Mukai's approach, 20 pairs of homozygotes and their hybrids are used.

Mukai'sapproach in selfing populations:Homozygous lines are readily obtainable. Simulations are performed as above, except that n in the fitness function W = (1 − s)ⁿ is randomly determined from a Poisson distribution with mean U/(2s) (Charlesworthet al. 1990; Deng and Lynch 1996a).

New approach: A number of randomly outcrossed parents are sampled, with each having fitness W = (1 − hs)ⁿ, where n is the number of mutations randomly determined from the Poisson distribution with mean U/(hs) (Deng and Lynch 1996a). A number of selfed progeny genotypes for each parent is obtained, with each selfed genotype being determined by allowing the n parental heterozygous loci to segregate randomly into AA, Aa, and aa classes with respective probabilities of 0.25, 0.5, and 0.25. Letting n₁ and n₂ (both resulting from random segregation) be the numbers of heterozygous and homozygous loci containing mutations in a selfed offspring, its fitness is W(n₁,n₂) = (1 − hs)^n₁ (1 − s)^n₂. Equation 3c is used to estimate h. For the new approach in outcrossing populations 20 parents are employed throughout, each having 50 selfed progeny genotypes.

Constant mutation effects with over(under)dominance: Under over(under)dominance, the key assumption for estimating h, that the frequency of one homozygote at any polymorphic locus is low, is unlikely to be valid in outcrossing populations; however, in at least two situations, it may hold and Mukai's approach should be applicable regardless of h value. The first is in highly selfing populations, where overdominance for mutations is unlikely to be responsible for the maintenance of genetic variability (Kimura and Ohta 1971). The second is when lines are obtained by generations of mutation accumulation from homozygous replicate lines. If the mutation rate per locus is low, overdominant mutants are unlikely to achieve high frequencies.

Generally, the distribution of the number of loci per genome with overdominant mutations is not clear. Simulations are performed for different distributions of the number of loci having (over)dominance mutations. The simulation procedures are the same as before for the dominance case, except that the distribution of the number of loci under (over)dominance is different. The results (Table 2) indicate little influence on the estimation by Mukai's approach under very different simulated distributions of the number of loci per genome having (over)dominant alleles. This is not unexpected, because the derivation of the approach is based on the one-locus results and within-family data and extended to multiple loci under additive mutation effects across loci. No assumption was made as to the distribution of the number of loci having (over)dominant alleles per genome. Therefore, as before (i.e., Poisson distribution of the polymorphic loci is used) except that we let the parameter h < 0 (overdominance) or h > 1 (underdominance), simulations are performed for Mukai's approach for selfing populations. For Mukai's approach using mutation accumulation lines, the principle is the same and the simulation and results are similar, and hence not presented.

Mixed dominance and overdominance:Lines from highly selfing populations or derived by mutation accumulations: It is possible that both dominance and overdominance underlie heterosis and that new mutations of either nature can occur in the genome. Some interesting questions are the following: In highly selfing populations, what is the major cause of heterosis? In the genome, what is the major type of new mutations for heterosis? Can we answer these questions from ĥ? Throughout, a circumflex (^) indicates an estimated value.

With both dominant and overdominant loci present, if fitness effects of individual loci are independent, as is the case for the multiplicative fitness function, the heterosis due to dominance and overdominance is independent. Let the mutation rate to deleterious alleles and constant mutational effects under dominance be U₁, h₁, and s₁, respectively. In large highly selfing populations, the expected heterosis (the ratio of the mean fitness of the outcrossed offspring generation W_o to that of the homozygous parental generation W_p), which is due to dominance (δ_d), is then (Charlesworthet al. 1990; Deng and Lynch 1996a) δd=exp(−U1h1)∕exp(−U1∕2)=eU1(0.5−h1) (6a)

New mutational occurrence in the genome most likely follows a Poisson distribution, whether it involves dominant or overdominant mutations. Throughout, mutations fitting the (over)dominance hypothesis will be referred to as (over)dominant mutations. In highly selfing populations, mutant alleles will be maintained by mutation-selection balance, regardless of their (over)dominance. This is because within a selfing line, frequent selfing will quickly bring any polymorphic locus into homozygous state. Under obligate selfing, different selfing lines are essentially reproductively isolated from each other, and thus overdominance will not contribute to the maintenance of genetic variability. Hence, as in the dominance case, we assume that the number of loci with overdominant mutants (n) (all in homozygous state) per genome in selfing populations is Poisson distributed with mean n¯ and constant effects h₂ and s₂. If the genomic mutation rate to overdominant (but less fit) allele a is U₂, it can be easily shown that at mutation-selection equilibrium, n¯=U2∕(2s2) (Charlesworthet al. 1990; Deng and Lynch 1996a). The expected heterosis due to overdominance (δ_o) is then δo=(∑n=0∞(1−h2s2)n2(n¯)ne−2n¯∕n!)∕(∑n=0∞(1−s2)n(n¯)ne−n‒∕n!)=(∑n=0∞(2n¯−2h2s2n¯)ne−2n¯∕n!)∕(∑n=0∞(n¯−s2n¯)ne−n∕n!)=exp(−2h2s2n¯)exp(−s2n¯)(∑n=0∞(2n¯−2h2s2n¯)ne−(2n¯−2h2s2n¯)∕n!)∕(∑n=0∞(n¯−s2n¯)ne−(n¯−s2n¯)∕n!)=exp(−2h2s2n¯)∕exp(−s2n¯)=e(1−2h2)s2n¯=eU2(1−2h2)∕2. (6b)

The total heterosis is δ = E(W_o)/E(W_p) = δ_dδ_o. The contribution to heterosis from dominance relative to overdominance can then be measured by the index α=δdδo=exp(U1(0.5−h1)−(1−2h2)s2n¯). (7) α is an important index that will be used more later. If α > 1, the primary cause of heterosis is dominance; if α < 1, it is overdominance; otherwise, dominance and overdominance contribute about equally to heterosis. The smaller the α, the larger the contribution to heterosis from overdominance.

Using Equation 7, we can determine the number of overdominant loci (N_o) when dominance and overdominance contribute equally to heterosis No=U1(0.5−h1)s2(1−2h2). (8) n¯>No,n¯<No,n¯=No , respectively, correspond to α < 1, α > 1, α = 1.

In simulations, the genome contains both dominant and overdominant loci, all at mutation-selection equilibrium. In the parental generation, the number of dominant loci in each individual is sampled from the Poisson distribution of mean U/(2s₁) (Charlesworthet al. 1990; Deng and Lynch 1996a), and that for the overdominant loci is determined from a Poisson distribution with mean n¯=[=U2∕(2s2)] . Simulations are then performed as before for Mukai's approach.

Homozygous lines constructed from outcrossing populations: In outcrossing populations, with overdominance the key assumption that the rarer allele at any polymorphic locus is of a low frequency is often invalid. Thus, both Mukai's and our new approaches should not be used. However, there have been some practice and data on estimating h¯ by Mukai's approach using homozygous lines constructed from outcrossing populations such as Drosophila (e.g., Mukaiet al. 1972; Mukai and Yamaguchi 1974; Eaneset al. 1985). Therefore, under mixed dominance and overdominance, we investigate whether those estimates are of any use and whether they can discern the predominant mode of genetic effects for heterosis.

Let U₁, h₁, s₁, h₂, s₂, δ_d, δ_o, α, and N_o be defined as before. In large outcrossing populations, upon selfing (Deng and Lynch 1996a) δd=exp(−U1)∕exp(−U1(h1+0.5)∕2h1)=eU1(1−2h1)∕(4h1). (9a) Overdominant alleles are mainly maintained by balancing selection in outcrossing populations. If there are n overdominant loci with constant effects h₂ and s₂ δo=(2(2h2−1)2−2h22s2(2h2−1)2(2h2−1)2−s2h2(2h22+h2−1))n. (9b) Then α=δdδo=(2(2h2−1)2−s2h2(2h22+h2−1)2(2h2−1)2−2h22s2(2h2−1))nexp(U1(1−2h1)4h1) (10) By Equation 9, a and b, we can determine N_o if we set α = 1 No=U1(1−2h1)4h1∗ln(2(2h2−1)2−2h22s2(2h2−1)2(2h2−1)2−s2h2(2h22+h2−1)). (11) Again, n > N_o, n < N_o, n = N_o, respectively, are equivalent to α < 1, α > 1, α = 1, which correspond to situationswhere dominance contributes to heterosis less than, more than, and equal to overdominance, respectively.

In the simulations, the homozygous lines constructed from large populations at mutation-selection equilibrium contain both dominant and overdominant loci. The number of loci homozygous for deleterious alleles (under the dominance hypothesis) is determined from a Poisson distribution of mean U/(2h₁s₁), as explained before. The alleles at the n overdominant loci are determined from random uniform variables (ξs) (from 0 to 1) with allele A being chosen if ξ ≤ h₂ − 1/2h₂ − 1 [where h₂ − 1/2h₁ − 1 is the equilibrium frequency of A allele (Crow 1986)], and allele a chosen otherwise. This simulation procedure allows identity by state for mutant alleles at the overdominant loci in different inbred lines. In outcrossed progeny, the number of the loci (all in heterozygous state) having dominant mutations is the sum of the number of mutations of its two homozygous parents; the genotypes at the n overdominant loci are determined by the overdominant alleles at these loci in its two homozygous parents. Other aspects being the same as before, simulations for Mukai's approach are then performed.

Variable mutation effects under dominance: Deleterious mutation effects across loci (s_i and h_i) are not constant. The few available data suggest that s_i has a roughly exponential distribution (Gregory 1965; Mackayet al. 1992; Keightley 1994). As in Deng and Lynch (1996a), we use the exponential distribution to model s_i f(si)=1s‒exp(−si∕s‒),(1>s>0) (12a) where s‒ is the mean of s_i. Little information exists on the distribution of h_i, but biochemical arguments suggest an inverse relationship between s_i and h_i (under dominance hypothesis), mutant alleles with larger effects tending to be more recessive (Kacser and Burns 1981). The few available data are consistent with this idea (Crow and Simmons 1983). Therefore, as in Deng and Lynch (1996a), the following function is adopted to approximate the inverse relationship between s_i and h_i hi=12exp(−13si), (12b) which is in rough accordance with the few available data (Crow and Simmons 1983; Deng and Lynch 1996a). Under the above distribution and function, the mean and the variance of h_i is h¯=1(26s‒+2)σh2=169s−2(26s‒+1)(13s‒+1)2 (12c)

To evaluate how Equations 1b and 4b perform for estimating σh2 and how Equations 1a and 3b behave in estimating h¯ under variable mutation effects, simulations are conducted under dominance. As in Deng and Lynch (1996a), we use a discrete version of the exponential distribution of Equation 12a by dividing the entire range of s_i (0, 1) into 200 classes of width 0.005. Each parental individual in a simulation is then randomly assigned a number of mutations from each of these classes by drawing from a Poisson distribution with means of Up_i/(h_is_i) in outcrossing and of Up_i/(2s_i) in selfing populations, respectively, where p_i is the density of the mutation distribution in the ith class.

Because the estimates are usually biased under the variable mutation effects, we compute their MSE (mean square error) for comparison: MSE=E(x^−E(x))2=Var(x^)+(x^¯−E(x))2 , where x^¯ stands for the estimated mean. Note that when x^¯ is unbiased, MSE is simply the variance of x^ .

The effects of lethals: The above study for variable mutation effects assumes that the genome contains no lethal mutations. This is a good assumption for selfing populations, where lethal mutations cannot survive for more than a few generations, due to frequent exposure to selection in homozygous state. In outcrossing populations, this assumption does not hold (Simmons and Crow 1977; Crow and Simmons 1983), as lethals are usually shielded from selection in heterozygous state by their low degree of dominance. With Mukai's approach, lethals will not appear in final homozygous lines constructed. During generations of inbreeding to construct homozygotes, lines homozygous for lethals will be immediately lost. Therefore, Mukai's approach gives the estimates only for mildly deleterious mutations. Hence, we only evaluate our new approach under variable mutation effects with lethals in outcrossing populations.

Lethals (s_L = 1) compose approximately 1% of the genomic mutations, and h_L for lethals is estimated to be about 0.02 (Simmons and Crow 1977; Crow and Simmons 1983). The simulations in outcrossing populations with lethals and variable mutation effects are identical in all respects to those in the previous section, but with an additional low genomic mutation rate (1%) to lethals (s_L = 1, h_L = 0.02). In simulations, selfed offspring homozygous for lethals will be excluded from analysis, as is most likely the practice in actual experiments. The selfed family means are for those offspring that do not contain homozygous lethals.