DISTRIBUTION OF EFFECT SIZES

Surprisingly, there is little theory about the suite of genetic differences that are likely to be observed among recently diverged taxa. While we have long known how to calculate the probability of fixation of a single mutation and how this depends on its effect on fitness, few studies have examined the entire set of substitutions likely to result from a period of adaptation.

Using an argument presented by Gillespie (1991, p. 266), we might expect an exponential distribution of selection coefficients among alleles that have fixed within a population. Gillespie noted that if all possible alleles at a locus were ordered in terms of absolute fitness from lowest to highest, then the selection coefficient of the second most fit allele relative to the most fit allele would follow an exponential distribution when measured across all loci. (His proof is based on extreme value theory.) Therefore, if the previously most fit allele at a locus falls only to second place in the fitness ordering when the environment (internal or external) changes, then the selection coefficients of substituted alleles should follow an exponential distribution. That is, we might expect an exponential distribution of selection coefficients among alleles that have fixed within a population if evolutionary change is accomplished primarily by the successive replacement of the two most fit alleles. The proof assumes, however, that genetic interactions are absent and that allelic fitnesses follow the same distribution at each locus, which is unlikely given functional differences among genes.

More recently, Orr (1998) studied the process of adaptation in a geometrical model first proposed by Fisher (1958). In this model, a fitness landscape in multiple dimensions has a single fitness optimum. Initially, the population is displaced from this optimum, perhaps because of a recent shift in the environment. Each new mutation that occurs affects a number of traits but does so to varying degrees. Whether the mutation is beneficial or deleterious depends on the current state of the population and on whether the mutation brings the population closer to the fitness optimum. Orr tracked the series of mutations that were incorporated into the population as it approached the adaptive optimum. In general, adaptation toward a peak in the fitness landscape produced a series of genetic changes whose effects on phenotype followed an approximately exponential distribution, a result that held under a broad variety of plausible distributions for the effects of mutations (Orr 1998, 1999).

Finally, an exponential distribution of fitness effects may also be a reasonable outcome of selection with a moving optimum. To show this, we assume that the fitness effects of new beneficial mutations follow a gamma distribution with mean, μ, and coefficient of variation, C. The shape of the gamma distribution depends on the coefficient of variation (Figure 1). When C = 1, the gamma is equivalent to an exponential distribution. For C > 1, the distribution becomes more L-shaped, while for C < 1, the distribution approaches a bell shape. We also assume that the distribution of newly arising beneficial mutations remains approximately constant throughout the adaptive process, which is plausible if a population lags steadily behind a moving optimum or fitness threshold. The distribution of fitness effects among those beneficial mutations that fix within a population can be calculated by weighting the original gamma distribution by each allele’s probability of fixation, which, in populations of large size, is approximately twice the allele’s selective advantage. It can be shown that the distribution of beneficial alleles, conditioned on their fixation, is also gamma, but with mean μ(1 + C²) and coefficient of variation C2∕(1+C2) . Interestingly, this conditional distribution always has a coefficient of variation <1. Hence, it is either exponential or bell shaped, even if newly arising beneficial mutations have a very L-shaped distribution. Unfortunately, we know little about the distribution of fitness effects of mutations in general and of beneficial mutations in particular. The distribution of spontaneous deleterious mutations is thought to be roughly L-shaped, with estimates of C ranging from 2 to 5 (Keightley 1994; Lynchet al. 1999). If newly arising beneficial mutations follow a distribution of similar shape, then those beneficial mutations that fix will follow a gamma distribution that is nearly exponential in shape, with a coefficient of variation between 0.9 and 1.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

—Potential distributions of allelic effects. Each curve describes a gamma distribution with mean μ = 1 but with different coefficients of variation (C). The QTL underlying a particular phenotypic difference represent draws from the appropriate distribution, as illustrated by the circles under the x-axis. Only those QTL above the threshold of detection (θ = 0.8, thin vertical line) are likely to be detected (solid circles). Those below the threshold are likely to remain undetected (open circles).

The above arguments suggest that an exponential distribution might often describe the distribution of effect sizes among alleles that arise and fix during adaptive divergence. There may, however, be circumstances under which alternative distributions are plausible. For example, one might be examining a trait that has diverged as a pleiotropic response to selection on other traits or through neutral drift. In these cases, the distribution should more closely reflect the distribution of effect sizes among new mutations. On the other hand, if one of the parental lines has been subject to rapid and strong selection over a very short period of time, then only large-effect mutations will have had enough time to fix. In this case the distribution of fixed effects may be more normal in shape with a mean offset from zero. Similarly, the fact that researchers choose traits and parental lines that are particularly divergent may bias the distribution of effect sizes such that a greater proportion of allelic effects are large. The gamma distribution is a natural choice to describe these various alternatives, because its shape is so flexible. Here, we focus on the exponential distribution, both because it has theoretical support and because the results are simpler to understand. We also provide a more general derivation that uses a gamma distribution to describe the underlying effects of alleles on the trait of interest and that may allow a more flexible approach to fitting real data.

THE CASTLE-WRIGHT ESTIMATOR

Historically, one of the most widely used methods for estimating the number of loci underlying a trait difference between two lines was developed by Castle (1921) and Wright (1968). Their method is based on the amount of segregation variance observed in controlled crosses. Intuitively, if there is one or very few Mendelian factors underlying a trait, then some F₂ individuals will have the same genotypes as the parents, and the phenotypic range observed among the F₂’s should span the entire range of the parental lines, P₁ and P₂. As the number of loci contributing to the trait increases, however, F₂ individuals will more often fall near the average of the parental lines as a result of the central limit theorem. Thus the amount of segregation variance, Var(S), estimated from hybrid (F₂, F₃, etc.) and backcross generations, contains information regarding the number of alleles underlying the phenotypic difference between two parental lines. If Δz is the difference between parental lines in the trait of interest, the Castle-Wright estimator for the number of underlying loci is n^CW=(Δz‒)28Var(S). (1) Equation 1 assumes that the underlying loci are unlinked and that they have equal effects. Consequently, (1) is said to estimate the “effective” number of underlying loci, equivalent to the number that there would be if all QTL were unlinked and had effects of equal magnitude.

Numerous improvements and extensions have been made to the Castle-Wright estimator (summarized in Lynch and Walsh 1998, chapter 9). In particular, Zeng (1992) analyzed a generalized model that allowed for linkage among loci and variation in their effects. He estimated the number of underlying loci to be n^CWZ=2c‒n^CW+C2(n^CW−1)1−n^CW(1−2c‒), (2) where c is the average recombination rate between loci, and C is the coefficient of variation for the distribution describing the additive effects. Given that allelic effects are likely to vary from locus to locus, Equation 2 should more closely estimate the true number of underlying loci than does the effective number presented in (1). Equation 2 requires, however, that the distribution of effects be specified. Fortunately, it can be used with a variety of distributions describing the additive effects of the QTL, including an exponential and a gamma distribution.

Sampling variances for nˆ_CW and nˆ_CWZ have been determined by Lande (1981) and Zeng (1992). Rather than reiterate their general findings, we focus here on the sampling variance relevant to our simulation study. We assume that the parental-line means are known with negligible error and that the trait has been measured under controlled environmental conditions (i.e., the parental difference is attributable to genetic differences). In this case, the segregation variance reduces to the variance among F₂ individuals, and the variance of the Castle-Wright-Zeng estimator becomes Var(nCWZ)=2(2c‒(1+C2)n^CW)2(NF2+2)(1−n^CW(1−2c‒))4, (3) where N_F2 is the number of F₂ individuals measured. Equation 3 can be used to construct confidence limits for nˆ_CWZ under the assumption that the error in the estimator is normally distributed.

A QTL-BASED ESTIMATOR

In light of the growing number of QTL studies, we develop an alternate method, based on data generated in a QTL mapping study, to estimate the total number of loci underlying a trait difference. In a typical QTL analysis, a handful of the loci that contribute to the trait difference are detected, each with an estimated additive and dominance effect and position. Our method takes advantage of the fact that the power to detect a QTL depends on the size of its effect. Given an expected distribution of effects, we can estimate the number of loci whose effects were too small to be detected (see Figure 1).

We first need to determine how the probability of detecting a QTL depends on its effect size. This power curve is approximately logistic in shape for a simple QTL study with one marker linked to one QTL (Figure 2). With N_F2 individuals scored in the F₂ generation, the probability of detecting a QTL rises at some point from near 0 to near 1. The more F₂ individuals are examined, the more steeply the curve rises. With multiple markers, however, significance is usually assessed using a permutation test (Churchill and Doerge 1994). Without knowledge of the form of the power curve in a typical QTL study with multiple markers, we approximate the power curve with a threshold function, rising from 0 to 1 at a threshold θ. This approximation should be most accurate when many F₂ individuals are scored. Other power functions can be explored, at least numerically, using the approach developed below, although simulations indicate that the estimator developed using this approximation is reasonably accurate as long as there are enough detected QTL (more than two) to provide a good estimate for the threshold.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

—The power to detect a completely linked QTL with one marker (solid curves), based on equations 15.36b and 15.37 of Lynch and Walsh (1998). The QTL was assumed to be completely additive with effect a, such that the expected difference between homozygotes at this locus was 2a. The type I error rate, α, was set to 0.05, and the phenotypic variance among F₂’s was scaled to 1.

We assume that the phenotypic difference between the parental lines, Δz, is primarily caused by fixed allelic differences at underlying QTL. We define the additive (a) and dominance (d) effects of each allele according to Zeng (1992),

Population:

P₁ F₁ P₂

Genotype:

A₁ A₁ A₁ A₂ A₂ A₂

Average phenotype:

a d -a,

where we refer to a as the effect size or additive effect of an allele. Here, we assume that each allele contributes to the parental difference in the same direction, that is, a is always positive or always negative. This assumption is reasonable if divergence is primarily a product of selection acting in different directions within the two populations or in a novel direction in one population. (In a later section, we explore the effect of relaxing this assumption through simulations.) We let D denote the sum of the additive effects across all QTL. Under the above assumptions, D also equals half the phenotypic difference between parental lines (Δz/2). We first derive an estimator for the number of loci underlying the trait difference assuming an exponential distribution of allelic effects and given an estimate of θ (the minimum threshold of detection). We then discuss how to obtain confidence limits for this estimator and how to estimate θ from QTL data. Finally, we repeat these steps assuming allelic effects follow the more general gamma distribution.

Exponentially distributed effect sizes: As argued above, we might expect alleles underlying phenotypic divergence to be exponentially distributed. To begin, we proceed under the assumption that the additive effect sizes (a) represent draws from an exponential distribution with mean, μ, and probability density function f[x,μ]=Exp[−x∕μ]μ. (4) If there are a total of n underlying QTL, then the mean effect size will be μ = D/n. Using (4), we can determine the probability density function for detectable QTL, i.e., those that lie above the threshold θ, fd[x,μ]=f[x,μ]∫θ∞f[x,μ]dx=Exp[−(x−θ)∕μ]μ, (5) which is simply the same exponential distribution with an origin shifted to the right by an amount, θ. The mean additive effect size among detected QTL is therefore expected to equal μ + θ = D/n + θ. If the number of QTL actually detected is n_d and their average effect size is M, we can estimate the total number of loci by setting M to its expectation, D/n + θ. Rearranging this equation, the estimated number of underlying loci is n^QTL=DM−θ. (6)

Equation 6 can also be derived using a likelihood approach, which will enable us to calculate confidence limits for the number of loci underlying a trait difference. From the probability density function for observable QTL, (5), the likelihood of observing a set of QTL whose additive effect sizes are a_i (i = 1... n_d) is proportional to L[n]=∏i=1ndfd[ai,μ]=Exp[−nd(M−θ)∕μ]μnd. (7) Substituting in μ = D/n and solving for the maximum, the most likely value for the true number of underlying loci is again given by (6).

Confidence limits: Confidence limits for nˆ_QTL can be obtained using the likelihood-ratio test, which holds that the values of n whose log likelihoods lie within χ1[α]2∕2 of the maximum log likelihood comprise an approximate 1 - α confidence region. This amounts to solving the equation (lnL[n^QTL]−lnL[n])=χ1[α]22 (8) for n. From (7), the confidence limits for n equal the two roots (found by numerical solution or plotting) to the equation −nd+(M−θ)nndD−ndln[(M−θ)nD]−χ1[α]22=0, (9) where χ1[0.05]2=3.841 is used to obtain 95% confidence limits. More exact confidence limits can be found using the methods described for exponential distributions in Larsen and Marx (1985, p. 313), but numerical examples suggest that the two limits are similar. Note that this approach does not constrain estimates for n to be greater than the number of detected loci; if desired, a constrained likelihood surface can be examined (with n ≥ n_d), and the confidence limits can be adjusted accordingly.

Average effect of the undetected QTL: Because of the Beavis effect, the effects of the detected QTL are inflated and may appear to explain more of the trait difference between two parental lines than is the case. Here, we use the exponential distribution to estimate the relative importance of the undetected QTL. Specifically, we estimate the expected effect size of the undetected QTL, M_undetected, by averaging the exponential distribution within the range, 0 to θ. Recalling that the mean effect of all QTL (μ) is estimated by M - θ, we obtain the estimate Mundetected=∫0θxf[x,μ]dx∫0θf[x,μ]dx=M(1−1−τ1−Exp[−(1−τ)∕τ]), (10) where τ = (M - θ)/M. τ is a measure of the power of a QTL experiment. It ranges from 0 for very weak experiments with a threshold near the mean detected effect to 1 for a very powerful experiment with a threshold near 0. Equation 10 indicates that the average effect among the undetected QTL reaches a maximum of 23% M when τ = 0.36. For experiments with little power (τ near 0), the relative size of the undetected QTL decreases toward zero because any QTL that are detected are likely to have effects well above the average. Similarly, for powerful experiments (τ near 1), most QTL will have been detected, and the remaining ones will have a relatively small effect. For the simulation experiments described below, the average power ranged from 0.27 to 0.61 (see Table 1), where the experiments with 500 F₂’s were more powerful. Over this range, the expected effect of the undetected QTL is 17-23% of the average effect of the detected QTL (M).

As an example, in the simulation study described below with 20 underlying QTL and 500 F₂’s (see Table 1), the detected QTL appeared to explain just over 100% of the parental difference, on average, even though less than half of the underlying QTL were detected. In fact, the undetected QTL accounted for ∼20% of the parental difference. Thus, the effects of undetected QTL can be substantial even when a QTL study indicates that the entire difference between parental lines has been explained.

One can also determine the expected fraction of segregation variance accounted for by the undetected QTL. From Lynch and Walsh (1998) Equation 9.26a, the segregation variance for unlinked loci equals ¹/₂ times the sum of the squared average effects of all the QTL. The expected fraction of the segregation variance contributed by loci whose effects lie below the threshold of detection, θ, can be shown to equal 1−12Exp[−1−ττ](1+1τ2). For example, in the simulation experiments, where the average power (τ) ranged from 0.27 to 0.61 (see Table 1), the fraction of the segregation variance that we expect to be accounted for by undetected QTL ranges from 51 to 3%, decreasing rapidly as the power of the experiment increases. Although one could simply calculate the total segregation variance contributed by the observed QTL, this procedure would overestimate the importance of the detected QTL because of the Beavis effect, which causes an especially large bias in variance calculations where the effect sizes of detected QTL are squared.

Estimating the threshold of detection: We have not yet addressed how to estimate the unknown threshold, θ. Ideally, one would compute the power curve for the probability of detecting a QTL. For example, one could perform Monte Carlo simulations, placing QTL of known effects on simulated genomes using the same number of markers and individuals as in the planned experiment. Alternatively, one can use the observed QTL data to estimate θ, as follows. If there are several detected QTL, then the magnitude of the smallest observed additive effect size (a_min) will often be near the threshold and can be used as an approximation for θ. Indeed, a_min is the maximum-likelihood estimator for θ because the shape of the exponential distribution is such that the most probable value for the smallest observed QTL is at the origin of the distribution (i.e., at θ). This is clearly a biased estimate for the threshold, because the smallest observed QTL will not lie below the minimum detectable size and will generally lie above θ. To obtain a less biased estimate for θ, we use a Bayesian approach. We assume that the threshold lies somewhere between 0 and a_min but that every point within this range is equally plausible (θ has a uniform prior distribution). We then weight the prior distribution for θ by the probability that the smallest detected QTL has additive effect size a_min. Noting that the minimum of n_d draws from an exponential distribution with mean μ is itself exponential with mean μ/n_d (Feller 1971, p. 18) and recalling that detectable loci follow an exponential distribution shifted to the right by θ, the probability density function for a_min - θ is exponential with mean μ/n_d. This gives us a posterior distribution for θ, whose average value is our estimated threshold (θ^) , which can be found to solve 0=amin−θ^−M−θ^nd+aminExp[−aminnd∕(M−θ^)]1−Exp[−aminnd∕(M−θ^)]. (11) If several QTL are detected, the last term becomes negligible, and the estimate for θ from (11) approaches θ^≈aminnd−Mnd−1. (12) Equation 12 makes sense: we would expect the smallest observed QTL to lie above θ by μ/n_d (because the distribution of the smallest QTL is an exponential with parameter μ/n_d shifted to the right by θ). Hence, E[a_min] = θ + μ/n_d = θ + (E[M] - θ)/n_d, which rearranges to give (11). With few detected QTL, however, Equation 12 is subject to substantial sampling error and can become negative, which is biologically unrealistic. Equation 11, however, is always bounded between 0 and a_min.

One potential problem with using a_min and M in our estimators is that they will be overestimated in QTL studies as a result of the Beavis effect. This occurs because the effect sizes of the QTL are estimated from the same data used to detect the QTL. Those QTL that, by chance, happen to have a larger effect than expected are more likely to be detected and so lead to inflated estimates of a_i. There is currently no analytical method to correct for the Beavis effect, so its impact was tested through simulation (described below). The simulations indicate that the Beavis effect did not cause a large bias in our QTL-based estimators. We suspect that nˆ_QTL is not extremely sensitive to the Beavis effect because both M and θ tend to be inflated, but the estimator depends primarily on their difference [see (6) and (9)], which may not be as strongly biased.

Gamma-distributed effect sizes: We now describe more general results that apply when allelic effects follow a gamma distribution with mean, μ, and coefficient of variation, C. The probability density function of a gamma distribution is g[x,μ,C]=Exp[−x∕(C2μ)]x(1−C2)∕C2(C2μ)−1∕C2Γ[1∕C2], (13) where Γ[x] is the gamma function (Abromowitz and Stegun 1972). Using the same methods as were used to derive (6), we can show that the estimated number of underlying loci, nˆ_QTL, equals ΨD/(M - θ), where Ψ solves ΨΓ[1C2,(1−τ)ΨτC2]=τC2Γ[1+1C2,(1−τ)ΨτC2]. (14) Here, τ = (M - θ)/M, and Γ[x, y] is the digamma function (Abromowitz and Stegun 1972). Ψ is essentially a correction factor for (6) that must be applied when the expected distribution of underlying effect sizes is gamma. If C = 1, the gamma distribution reduces to the exponential distribution, and Ψ becomes one. For other values of C, the correction term depends only on C and τ and is illustrated in Figure 3. When C > 1, the gamma distribution is L shaped, and more loci have very small effect. These small-effect loci are missed if one mistakenly assumed an exponential distribution and used Equation 6. Therefore, (6) underestimates the number of underlying loci (Ψ > 1). The converse is true when C < 1, and the gamma distribution is bell shaped. In this case, fewer loci have small effect, and fewer loci fall below the threshold of detection. Consequently, (6) overestimates the number of underlying loci (Ψ < 1). The sensitivity of the estimator to the shape of the distribution of allelic effects depends strongly on the power of the experiment (τ). When the threshold of detection is high and τ is near 0, the estimator is very sensitive to the shape of the distribution. Conversely, when the threshold of detection is low and τ is near 1, the correction term approaches 1, because most alleles are detectable regardless of the shape of the distribution.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

—The correction factor, Ψ, as a function of the coefficient of variation of the gamma distribution, C. The correction factor when multiplied by (6) gives the estimated number of underlying loci when the distribution of allelic effects follows a gamma distribution. It depends only on C and on τ = (M - θ)/M, which varies from zero for an experiment with little power to detect a QTL to one for an experiment with high power.

Again, a_min can be used as an upwardly biased estimator for θ. Alternatively, one can correct θ using a Bayesian approach that takes into account the probability that a_min is the smallest detectable QTL when effect sizes follow a gamma distribution. This leads to an expression involving θ that must be numerically evaluated simultaneously with (14): θ^≈∫0aminxΓ[1∕C2,xΨ∕C2(M−θ^)]−nddx∫0aminΓ[1∕C2,xΨ∕C2(M−θ^)]−nddx. (15) When C = 1, (15) reduces to (11).

SIMULATIONS

To assess the accuracy of our estimator and the Castle-Wright-Zeng estimator, QTL analyses were simulated using the QTLcartographer package (version 1.13a; Bastenet al. 1996). Using the Rmap program, we created a map with 20 chromosomes, each carrying five markers spaced 10 cM apart, with the first and last markers placed at the ends of the chromosome. This map was used for all simulations. On this map, we randomly placed 2, 5, 10, 20, or 100 QTL using Rqtl. For a given number of QTL, 30 different QTL maps were generated (i.e., we constructed 30 different “experimental genomes”). In the basic set of simulations, the following assumptions were made: (1) alleles were additive; (2) all allelic effects, a, had the same sign (i.e., each QTL acted in the direction of the observed parental difference); and (3) the QTL effects were drawn from an exponential distribution with a mean of 1/4. These assumptions were then relaxed in turn. To assess the effects of dominance, we used Rqtl to produce QTL with additive effects drawn from an exponential distribution with mean 1/4 and dominance effects drawn from a Beta distribution with shape parameter set to 1 (see Bastenet al. 1996). This routine produced a dominance term that ranged from -a to a, with mean equal to 0. That is, alleles from P₁ ranged from recessive to dominant. To assess the effects of having some QTL act in the opposite direction to the observed parental difference, we let the additive effect of a QTL be positive with probability 0.85 and negative with probability 0.15. Finally, to assess the effect of the underlying distribution of effect sizes, we drew effect sizes from a gamma distribution with a coefficient of variation of 0.5 or 2.0 instead of an exponential distribution (see Figure 1). For these extensions, the number of underlying QTL was set to 20, but otherwise the experimental design was the same as for the basic simulations. In all cases, environmental variation was assumed to be negligible.

For each experimental genome generated, Rcross was used to simulate the production of 200 and then 500 F₂ individuals. This F₂ sampling procedure was repeated 10 times. This design allowed us to assess whether the estimators were more sensitive to the realized distribution of QTL within the genome (variation among “experimental genomes”) or to the specific set of F₂ individuals analyzed (variation among sampled individuals). For each set of F₂ individuals, the QTL were mapped using Zmapqtl’s interval mapping routine (model 3). A single permutation test was used to determine an approximate significance threshold for each number of QTL and F₂ sample size (5% genome-wide significance level). In general, there was not much variation in significance thresholds (data not shown). This threshold was then used in conjunction with Eqtl to estimate the location and effects of the QTL (Bastenet al. 1996). Eqtl identifies putative QTL by scanning through the output of Zmapqtl and noting all intervals with a likelihood peak that exceeds the significance threshold.

Unfortunately, interval mapping methods systematically bias the effects of chromosome regions linked to QTL (Lynch and Walsh 1998, pp. 457-458). Put simply, an interval not containing a QTL but adjacent to one containing a QTL may exceed the significance threshold because of its linkage to the QTL. This problem is endemic to all QTL analyses and no adequate solution has yet been developed (Whittakeret al. 1996; Goffinet and Mangin 1998). To reduce this problem, we screened our data set for obvious “ghost” QTL. That is, we eliminated intervals containing putative QTL whose position was within another QTL’s approximate support interval (the chromosome interval within which the likelihood of a QTL was within a factor of 10 of the peak likelihood for that QTL; see Lynch and Walsh 1998, pp. 448-449). In each case, the QTL with the weaker support was eliminated. Visscher et al. (1996) have shown that a bootstrap approach is much more accurate, but it would have been too time consuming for our analysis.

We chose interval mapping for two reasons. First, interval mapping is an established method that has been repeatedly used in empirical studies and is well explored theoretically. Second, interval mapping is much faster, allowing us to perform more extensive tests. We also suggest that interval mapping provides a conservative test of nˆ_QTL. Interval mapping produces data that are less precise than more advanced methods of QTL analysis, such as composite interval mapping and multiple interval mapping (Zeng 1994; Kaoet al. 1999); this reduced precision should result in less accurate estimates for nˆ_QTL. Furthermore, interval mapping corresponds less well to our model, which assumes that QTL above the threshold of detection are always detected. Consequently, it seems reasonable that our estimator would perform even better with data from more sophisticated methods for detecting QTL.

Our simulation data were imported into Mathematica 3.0 (Wolfram 1996; package available at www.zoology.ubc.ca/~otto/Research) for analysis. In our calculations, we used the estimators appropriate for an exponential distribution of underlying QTL effect sizes. Specifically, we used Equation 6 for the QTL-based estimator and Equation 2 for the Castle-Wright-Zeng estimator, with C = 1. For nˆ_CWZ, we initially used the average recombination rate between randomly chosen pairs of loci, c, which equals 0.48 for a genome with 20 chromosomes, each of length 40 cM (Lynch and Walsh 1998, Equation 9.3). Unfortunately, this often led to nonsensical answers, especially when the number of underlying loci was large (>5). The problem lies in the fact that the denominator in (2) is subject to substantial sampling error; relatively often, the denominator approached zero (causing severe overestimates) or became negative. With 20 underlying QTL, the nˆ_CWZ estimate for the number of underlying loci averaged 49.3 (with 500 F₂’s) and -21.8 (with 200 F₂’s) over the 300 simulations! To avoid these problems, we set c to ¹/₂ in all of our analyses, which either made little difference (for small n) or improved the estimates.

Because the QTL estimator relies on the difference between the mean effect of factors found and the minimum effect found, it can only work if at least two QTL have been identified. Therefore, we excluded QTL analyses with only one detected QTL. In experimental genomes with more than two true QTL, very few analyses were excluded. With only two true QTL, however, approximately half of the analyses had to be excluded. Furthermore, in a very small number of cases with two true QTL (6/600), two QTL were detected whose estimated effects were, by chance, equal. To avoid division-by-zero errors, the mean was increased by 5% over the minimum in these cases. (If, instead, these cases were eliminated, the performance of the estimators improved slightly over the results shown in Table 1.)