THEORY

The assumptions are the same as those of the Morton-Charlesworth method (Mortonet al. 1956; Charlesworthet al. 1990) and the Deng-Lynch method (Deng and Lynch 1996, 1997; Deng 1998b). Namely, the population is assumed to be large, randomly mating, highly selfing or outcrossing, at linkage equilibrium, and at M-S balance. In addition, the fitness function is assumed to be multiplicative, which is biologically plausible (Mortonet al. 1956; Crow 1986; Craddocket al. 1995; Fu and Ritland 1996). Mutations at each locus are assumed to have constant effect s and h.

In this study, we consider variable mutation effects in the development of an estimation method for DGM parameters in natural populations. Under variable mutation effects across loci, homozygous effect s for mutations is a random variable between 0 and 1. We assume that, for a mutation, dominance coefficients h and s are functionally related so that h = h(s). This assumption is supported by the limited data and theory (Simmons and Crow 1977; Kacser and Burns 1981; Crow and Simmons 1983). We divide the domain of s, [0, 1], for new mutations into T intervals with each having a width of 1/T. Let I_k = [k/T, (k + 1)/T] denote the kth interval, and define the probability pk=P(s∈Ik),k=0,1,…,T−1. When T is sufficiently large, s and h are approximately constant within each interval but are variable across various intervals. Let U_k denote the mutation rate corresponding to mutations with an effect s falling into the interval I_k, and then U_k = Up_k.

With the assumptions we have, in outcrossing populations, the number of mutant alleles with mutation effects s falling into an interval I_k within an individual (all in the heterozygous state; Mortonet al. 1956; Deng and Lynch 1996) follows a Poisson distribution with an expectation n‒k=Uk∕hksk=Upk∕hksk (1a) (Deng and Lynch 1996, 1997). In selfing populations, the number of loci homozygous for mutant alleles with an effect s falling into an interval I_k within an individual follows a Poisson distribution with an expectation n‒k=Uk∕2sk=Upk∕2sk. (1b)

Outcrossing populations: We illustrate our experimental design and estimation method by using populations capable of selfing. The method may be extended to outcrossing populations where selfing is not feasible as in the Deng-Lynch method (Deng 1998b). The basic data structure is outcrossed parents and multiple selfed progeny from each parent (forming selfed families). Let W_o and W_s be the mean fitness in the parental and offspring generations, respectively, σ²_o the genetic variance of fitness in the parental generation, σ²_t the total genetic variance of fitness in the selfed progeny generation, σ²_s the genetic variance of the mean fitness of selfed progeny in selfing families, and cov(w_p, w_s) the covariance between the fitness of a parent (w_p) and the mean fitness of its selfed progeny (w_s). Under the above assumption for mutation effects that are variable across various intervals at different loci, as in Deng and Lynch (1996), it can be shown that the fitness moments are related to the DGM parameters as W‒o=Wmaxexp(−U) (2) σo2=W‒o2[exp(Uhs‒)−1], (3) W‒s=Wmaxexp{−(U∕4)[2+(1∕h∼)]}, (4) σs2=W‒s2{exp[(U∕4)(s‒+hs‒+s∕(4h)‒)]−1}, (5) σt2=W‒s2[exp[U(hs‒∕2+s∕(4h)‒)]−1], (6) cov(wp,ws)=W‒sW‒o{exp[(1∕4)U(2hs‒+s‒)]−1}, (7) where the parameters with overbars denote arithmetic mean properties of new DGM parameters,h˜ is the harmonic mean dominance coefficient of new mutations, and W_max is the expected fitness of a mutation-free genotype in an environment where fitness measurements are taken. W_max serves as a scaling factor so that the fitness measurement can be on any scale instead of just from 0.0 to 1.0 and also so that mean environmental effects of experiments do not influence estimation (Deng and Lynch 1996).

Among Equations 2, 3, 4, 5, 6, 7, there are only five independent equations containing six unknown parameters. By assuming one of the six parameters known in the estimation, estimators of the other parameters can be derived. This is the strategy employed in the likelihood characterization of DGM parameters when variable mutation effects are considered in estimation (Keightley 1994; Denget al. 1999; Deng and Li 2001). Here we assume that U is known in the estimation for the time being. Alternatively, an initial value of U may be estimated from other approaches (Denget al. 1999) or may be estimated by the current experimental design and data with the Deng-Lynch method (Deng and Lynch 1996; see below). (If we assume that one of the parametersh˜, s, and hs is known, similar estimation procedures can be derived for U and the rest of the other parameters.h˜ can be estimated by methods such as that of Deng 1998a.) Solving these equations jointly yields estimators of h˜_, hs‒ , and s and as h∼=12−4(y∕U),s‒=4b−2xU,hs‒=xU, (8) where x=ln(σo2W‒o2+1),y=ln(W‒sW‒o),z=ln(σs2W‒s2+1),b=ln(cov(wp,ws)W‒oW‒s+1). (9)

From these estimates, other composite parameters of DGM, such as the mean number of mutations per genome n, mutational variance V_m per generation, and mean mutation effects on fitness Uhs‒ , can be derived (Deng and Lynch 1996). The covariance of h and s for mutations cov(h, s) can be approximated, or at least an upper bound can be estimated, as cov(h,s)=hs‒−h‒s‒≤hs‒−h∼s‒. (10) This is because for any distribution h¯ ≥ h¯. Let cov(h,s)=hs‒−s‒h∼ , where cov(h, s) denotes an upper bound of cov(h, s). This offers us the first opportunity to quantify the magnitude and the sign of cov(h, s). It would be impossible to come up with analytical estimators for DGM parameters such as cov(h, s) if in the analytical derivation, variable mutation effects are not considered. This is simply because these parameters such as cov(h, s) would be zero and meaningless in an analytical estimation developed under constant mutation effects.

Selfing populations: Random pairs of highly selfing and homozygous parental genotypes (denoted as P generation) are crossed to obtain outcrossed progeny (denoted as F₁ generation). Let W¯_p and σp2 be the mean fitness and genetic variance of fitness in the P generation, respectively, W¯_F1 and σ²_F1 be the mean fitness and genetic variance of fitness in the F₁ generation, respectively, and cov(P, F₁) be the covariance between the mean fitness of the two parents and the fitness of their F₁ progeny. Under variable mutation effects across loci, the fitness moments are related to the DGM parameters as follows: W‒p=Wmaxexp(−U∕2), (11) σp2=W‒p2[exp(Us‒∕2)−1], (12) W‒F1=Wmaxexp(−U∕h‒), (13) σF12=W‒F12(exp(Uh2s‒)−1), (14) cov(P‒,F1)=W‒pW‒F1[exp(Uhs‒∕2)−1]. (15)

It should be noted that the derivation for Equations 2, 3, 4, 5, 6, 7 and 11, 12, 13, 14, 15 assumes mutation effects that are variable. The strategy is to divide the range of variable selection coefficient s (from zero to one) into infinitely small intervals so that s can be treated as constant within each of the intervals but varying across intervals in our analytical derivation. Again, there are six unknowns (U, h, hs‒ , s, h2s‒ , and W_max) in the above five equations. By assuming or estimating one of the six parameters, estimators of the other five parameters can be derived. Here, as earlier for outcrossing populations, we assume that U is known in the estimation for illustration. Alternatively, an initial value of U may be estimated from other approaches (Denget al. 1999) or may be estimated with the Deng-Lynch method from the same data and experimental design as the current estimation method (Deng and Lynch 1996). Solving these equations jointly yields estimators of h, hs‒ , and s, h‒=0.5−(y∕U),s‒=2xU,hs‒=2bU, (16) where x=ln(σp2W‒p2+1),y=ln(W‒F1W‒p),z=ln(σF12W‒F12+1),b=ln(cov(P‒,F1)W‒pW‒F1+1). (17)

In selfing populations, we can use Equation 10 to estimate cov(h, s) by the above estimates of h, hs‒ , and s, which are unbiased under variable mutation effects with a known correct U. The estimators for h and s when assuming U is known are the same as those in Deng and Lynch (1996) for selfing populations.

The above estimation developed herein does not assume any specific functional relationship between s and h and any specific distribution form for the selection coefficient s. Therefore, the estimates are robust to different unknown forms of the distribution of s and the functional relationship between s and h. This is true despite that we assume specific distributions of s and a functional relationship between s and h in the following simulation studies to investigate the statistical properties of our estimation.

SIMULATIONS AND RESULTS

As with Keightley (1994), we assume that s for mutations follows a gamma distribution, with a density function g(s)=αβsβ−1e−αs∕Γ(β), where Γ(β)=∫0∞yβ−1e−y dy. α and β are the scale and shape parameters, respectively. s¯ =β/α and σs2=β∕α2 . As in Deng and Lynch (1996), we let h = h(s) = e^-As/ 2, where A = 13, which is in rough accordance with the few available data (Gregory 1965; Mackayet al. 1992; Deng and Lynch 1996; Deng and Fu 1998). With these assumptions, the parameters h,h˜_, hs‒ , and cov(h, s) can be derived as h‒=αβ∕[2(A+α)β],h∼={(α−A)β∕[2αβ]whenα−A>00whenα−A≤0,}hs‒=βαβ∕[2(A+α)β+1],cov(h,s)=−Aβαβ−1∕[2(A+α)β+1]. These DGM parameters can be used for comparison to examine the estimated values with our estimation methods in simulations.

The simulation procedures are the same as those that have been documented extensively earlier (Deng and Lynch 1996; Deng 1998b) and are thus not elaborated here. In simulations, we assume that the fitnesses of various genotypes can be measured with little error, which is justifiable in the investigation of estimation bias and comparison of various estimation methods (Denget al. 1999). Under the assumptions for the analytical development of our estimation methods, the number of mutant alleles corresponding to an interval I_k per individual follows the Poisson distributions (Equations 1a and 1b) with p_k being determined as pk=P(s∈Ik)=1Γ(β)∫Ikαβsβ−1e−αsds,k=0,1,…,T−1. It can be shown that pk=−exp(−αk+1T)+exp(−αkT),whenβ=1;pk=−[exp(−αk+1T)](αk+1T+1)+[exp(−αkT)](αkT+1),whenβ=2 when β= 2 and when β= 0.5, pk=Erf(αk+1T)−Erf(αkT), where Erf(x)=(2∕π)∫0xe−t2dt(x>0) . Erf(x) can be approximated as Erf(x)≈1−(1+∑i=16aixi)−16 (18) (Gao 1995), where a₁ = 0.0705230784, a₂ = 0.0422820123, a₃ = 0.0092705272, a₄ = 0.0001520143, a₅ = 0.0002765672, a₆ = 0.0000430638.

To evaluate the performance of our estimation in outcrossing populations in simulations, for each set of parameters U, α, and β, K parents were sampled from the parental generation, and from each of these, M selfed progeny were produced. The fitness of an individual from the parental generation is Wo=Wmax∏k=1T(1−hksk)nk, where n_k is the number of mutation-bearing loci with their effects falling into the interval I_k in an individual, obtained by random sampling from the Poisson distribution defined above. The fitness of each selfed offspring was obtained by allowing the n_k heterozygous loci of a parent to segregate randomly into the AA, Aa, and aa genotypes with respective probabilities of 1/4, 1/2, and 1/4. Letting n_1k and n_2k (k = 1,..., T) be the numbers of heterozygous and homozygous loci containing mutations with effects falling into the interval I_k in a selfed offspring, the fitness of the selfed progeny is Ws=Wmax∏k=1T(1−hksk)n1k(1−sk)n2k. Unless otherwise specified, for each set of parameters (U, α, β, K, M), we performed 1000 simulations. We let W_max = 1 throughout, as the value of W_max does not influence DGM parameter estimation.

For selfing populations, the fitness of an individual from the parental generation is Wp=Wmax∏k=1T(1−sk)nk, where n_k is the number of mutation-bearing loci with mutation effects falling into the interval I_k in an individual, and it is obtained by random sampling from the Poisson distribution defined earlier. Each parent mates with another random parent (not in the original set of K) to produce a total of K progeny (one per family) with fitness WF1=Wmax∏k=1T(1−hksk)n1k+n2k, where n_1k and n_2k (k = 1,..., T) are the numbers of homozygous mutant loci in interval I_k in the two parents, respectively.

In the estimation Equations 8 or 16, U must be known, assumed, or estimated with other approaches first. In simulations, we experimented and examined two methods to estimate U: (1) by the Deng-Lynch method (Deng and Lynch 1996) and (2) by an empirical regression procedure introduced here. We simulated parents and their children according to variable effects for each set of given parameter values of U, α, and β, and obtained the estimatesÛ_1,ŝ₁, andĥ₁ by the Deng-Lynch method (Deng and Lynch 1996). (A circumflex indicates an estimated value throughout.) We found a strong linear relationship between the parameter values of U and the estimatesÛ₁ andŝ₁ under any fixed β. Through a series of simulations, we obtained samples under various parameter values of U, α, and fixed β-values, and we obtained estimatesÛ₁ andŝ₁ with the Deng-Lynch method under various fixed β-values. Then we fit a multiple regression model under each specific β-value, U^=a^1+b^1U^1+c^1s^1, (19) whereÛ estimates U with little bias when β is correctly assumed as shown by our simulation results not presented here. The empirical estimation is useful only when the shape parameter β can be estimated using other methods and experimental data (e.g., Keightley 1994).

The simulation results are represented by the data in Tables 1, 2, 3, 4. The ranges of the values for the parameters (such as U, h, and s) generally cover those reported earlier from classical empirical experiments (e.g., Mukaiet al. 1972; Lynchet al. 1999). Three general conclusions emerge from our simulation studies under variable mutation effects. First, when U is set to equal true values or when the estimates of U are obtained via Equation 19 by assuming a correct β-value, application of Equation 8 or 16 to both obligate selfing or outcrossing populations yields nearly unbiased estimates for the DGM parameters with small standard deviation. The estimates of U by Equation 8 have smaller mean square error despite larger standard deviation when U is set equal to the estimates obtained by regression Equation 19 than those obtained by the Deng-Lynch method. The larger standard deviation may be partly due to the fact that Equation 19 is established by empirical regression procedures that involve an additional level of sampling error for the final estimation. The estimates of s by Equation 8 in outcrossing populations have smaller sampling variance and smaller bias than those obtained directly by the Deng-Lynch method, e.g., by comparison of the estimates in rows 1 and 3 for each parameter set in Table 1. This is true even when no prior assumption is made about the magnitude of U, when U is first estimated directly with the Deng-Lynch method, and then the estimate of U is used in the current estimation method, (Equation 8) for the other DGM parameters. The estimates ofh˜ by Equation 8 have smaller or comparable sampling variance than those obtained directly by the Deng-Lynch method for h (for each parameter set, compare the estimates of the second to fourth rows with that of the first row in Table 1). The comparison of the estimation quality between the current estimation method and the Deng-Lynch method changes little with the parameter values (Table 1). When β= 0.5, the bias of the estimates of the parameters is larger than that when β= 1 and 2. This may be due to the approximation formula 18 used to compute p_k = P(s_i ∊ I_k) when β = 0.5, while the computation of p_k = P(s_i ∊ I_k) when β= 1 and 2 is exact.

Second, when U is set equal to the estimates (Û₁₎ that were obtained by the Deng-Lynch method (Deng and Lynch 1996) and that are downwardly biased, the estimates of the other DGM parameters by Equations 8 and 16 are biased with small sampling variance (Tables 1 and 2). For outcrossing populations, the estimation Equation 8 yields less biased estimates with smaller standard deviation for s than for the Deng-Lynch method (Table 1), and the estimates of s, hs‒ , cov(h, s) are upwardly biased and estimates ofh˜ are downwardly biased. The result can be understood from Equation 8, sinceÛ₁ is downwardly biased as estimated by the Deng-Lynch method. In selfing populations, Equation 16 yields the same estimates for s and h as those obtained by the Deng-Lynch method (Table 2), which is expected as pointed out earlier. The estimates of s, hs‒ , and cov(h, s) are upwardly biased and estimates of h are downwardly biased because Û₁ is downwardly biased, which can be understood from Equation 16.

Third, in outcrossing populations, the cov(h, s) is correctly estimated to be an upper bound of cov(h, s); however, the sign of cov(h, s) can sometimes be estimated to be different from that of cov(h, s). In selfing populations, cov(h, s) can always be estimated with correct sign and small estimation bias.

ROBUSTNESS ANALYSIS

In the estimation of the DGM parameters, we need a prior estimate of one of the six parameters (such as U as investigated here) based on some external knowledge obtained from other estimation approaches. The estimation bias of this parameter or the bias of an assumed value will cause estimation bias of the other parameters. Hence, we investigate the sensitivity of estimators to the departures of U from true value, using computer simulations (Figures 1 and 2). We define a relative bias rate (RBR), (estimate true value)/(true value), to measure the sensitivity of estimators to an incorrectly assumed or estimated U value. In examining the robustness of the estimator for cov(h, s), the true value used is the parameter value of cov(h, s) as defined after Equation 10 and not cov(h, s).

In simulations for the investigation of the robustness of our current estimation of the other DGM parameters, U is set equal to a given value (denoted as U_given), which ranges from 0.5U₀ to 1.5U₀ (U₀ is the true value of U). This range of the estimate of U investigated is reasonable given the magnitude of bias that is normally found with the method such as that of Deng and Lynch (1996). The changes in the mean relative bias rates (MRBR) of the estimates of the parameter values in 1000 simulations are shown in Figures 1 and 2. It can be seen that when U_given ranged from 0.7U₀ to 1.5U₀ (which means that the departure of U_given from U₀ ranged from -0.3U₀ to 0.5U₀), the MRBR of the estimates of the parameter values changed smoothly and changed little in both outcrossing and selfing populations. When U_given ranged from 0.9U₀ to 1.2U₀, the absolute values of the MRBR of the estimates of parameters [except cov(h, s) for outcrossing populations when α= 20] are <0.185 in both outcrossing and selfing populations. For outcrossing populations, when α= 20, if U_given ≤ 0.9U₀ or U_given ≥ 1.1U₀, the absolute values of the MRBR of cov(h, s) are >1.0 (Figure 1, b and d). (Note the scale difference of the y-axis in Figure 1, b and d, with the other plots in Figures 1 and 2.) Thus, even when U is estimated with some bias, if the magnitude is similar to that obtained by methods such as that of Deng and Lynch (1996), our current estimation method can generally still yield relatively robust estimates of DGM parameters (except cov(h, s) for outcrossing populations when α is as small as 20). In outcrossing populations, the MRBR changed the sign in the robustness investigation of cov(h, s) when s = 0.01 and 0.05, respectively. This is because the parameter value cov(h, s) changed the sign from negative to zero and then to positive values under the functions assumed when s changes from 0.047 to 0.048.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

—The changes in RBR of the estimates of s,h˜_, , cov(h, s) obtained by Equation 8 in outcrossing populations when U were given equal to the values that ranged from 0.5U₀ to 1.5U₀. Each data point was the mean in 1000 simulations with the following sets of parameters and β= 1.0: (a) U₀ = 1.5, s = 0.01, and α= 100; (b) U₀ = 1.5, s = 0.05, and α= 20; (c) U₀ = 0.5, s = 0.01, and α= 100; and (d) U₀ = 0.5, s = 0.05, and α= 20.

APPENDIX: ESTIMATION OF OTHER DGM PARAMETERS WHEN h IS ASSUMED OR ESTIMATED AND SOME REPRESENTATIVE SIMULATION RESULTS

Ifh˜ (in outcrossing populations) or h (in selfing populations) is known by other estimation methods or assumed at particular values on the basis of some external knowledge, based on Equations 2, 3, 4, 5, 6, 7 and 11, 12, 13, 14, 15, we have estimators for other DGM parameters as follows, the notations being the same as in the text, in outcrossing populations, U=4y2−(1∕h∼),s‒=4b−2xU,hs‒=xU. (A1) and in selfing populations, U=y(0.5−h‒),s‒=2xU,hs‒=2bU. (A2)

Simulations are performed similar to that described in the text and with the above estimation for other DGM parameters whenh˜ (in outcrossing populations) or h (in selfing populations) is known or estimated. The simulation and the experimental procedures, whenh˜ (in outcrossing populations) and h (in selfing populations) are estimated by the methods of Deng (1998a) or Mukai et al. (1972), are detailed in Deng et al. (1998) and thus are not elaborated here.

Some representative results are presented in Tables A1 and A2. It can be seen that, relative to the Deng-Lynch method, the new method developed here can estimate more parameters, such as cov(h, s) and its sign. In an outcrossing population, the sign of cov(h, s) cannot be reliably estimated. However, in selfing populations, if the h is estimated first by the Deng-Lynch method and then used in the current method, the sign of cov(h, s) can be characterized correctly.