Multiple Interval Mapping for Quantitative Trait Loci

GENETIC MODEL

Consider m QTL, Q₁, Q₂, ···, and Q_m, in a backcross population in which there are two genotypes, Q_jQ_j and Q_jq_j, each with one-half frequency for a QTL, say Q_j. For m QTL, there are 2^m possible different QTL genotypes in the population. Cockerham's genetic model (C-H. Kao and Z-B. Zeng, unpublished results) is used to define the genetic parameters and model the relation between the genotypic value and the genetic parameters. If only up to digenic epistasis is considered, the relation between the genotypic value of individual i, G_i, and the genetic parameters can be expressed in the equation Gi=μ+∑j=1majxij+∑j<kmwjk(xijxik),i=1,⋯,2m, (1) where x_ij is coded as ½ or –½ if the genotype of Q_j is Q_jQ_j or Q _jq_j, respectively, a_j is the corresponding main effect of Q_j, and w_jk is the epistatic effect between Q_j and Q_k. The main advantage of Cockerham's model is that it possesses an orthogonal property in modeling genetic parameters.

To assist with explaining the estimation of the genetic effects in the MIM model (Equation 3), the genetic model in Equation 1 is expressed in matrix notation as Equation 2 (Scheme 1). In Equation 2, the column vector G contains the genotypic values of the 2^m possible genotypes. The subscripts of G (1 or 0) denote the homozygote or heterozygote of the QTL in the order of the first, second, third, ···, and mth QTL, respectively. The first m columns in the genetic design matrix D are the coefficients associated with the main effects of the m QTL, and the last m(m – 1)/2 columns represent the coefficients of the epistatic effects among them. Vector E contains the QTL main and epistatic effects. If there is no epistasis between some QTL, some of the columns for epistasis should be dropped out from matrix D. If higher-order epistasis is considered, the dimension of matrix D is easy to expand accordingly. The matrix D plays an important role in estimation of genetic parameters in the MIM model.

STATISTICAL MODEL OF MIM

Multiple interval mapping: Assume m QTL, Q₁, Q₂, ···, and Q_m, located at positions p₁, p₂, ···, p_m in m different marker intervals, I₁, I₂, ···, I_m, along the genome, control a quantitative trait y. Among the m QTL, some may show epistasis and some may not. The quantitative trait value for an individual, i, can be related Scheme 1 to the m putative QTL by the model yi=μ+∑j=1majxij∗+∑j≠kmδjk(wjkxij∗xik∗)+εi, (3) where μ is the mean, xij∗ is the coded variable for the genotype of Q_j, a_j and w_jk have the same definitions as those in the genetic model in Equation 1, δ_jk is an indicator variable for epistasis between Q_j and Q_k, and ϵ_i is assumed to follow N(0, σ²). Indicator variable δ_jk takes value one if Q_j and Q_k interact; otherwise its value is zero. In this model, the first summation is for the main effects of the m QTL, the second summation is for their possible epistasis, and ϵ_i is the environmental deviation. This is termed the MIM model because multiple (m) marker intervals are simultaneously used to construct multiple (m) putative QTL in the model for QTL mapping. If QTL genotypes are known, the model tells that the quantitative trait value is the sum of the QTL main effects, their possible epistatic effects, and environmental deviation, and the MIM model is a regression model. However, the putative QTL genotypes denoted by xij∗’s are usually not observed because QTL could be located in the intervals. Given observed flanking marker genotypes, the conditional distributions of QTL genotypes, xij∗’s , for QTL at specific positions, p_j's, can be inferred based on Haldane's mapping function (Haldane 1919) assuming no crossover interference (Table 1 in Kao and Zeng 1997), and the MIM model is then a normal mixture model. For each Q_j, its conditional probabilities are extracted to form a matrix Q_j (note that Q denotes QTL and Q denotes the conditional probability matrix; see Kao and Zeng 1997). The conditional probability matrices, Q_j's, j = 1, 2, ···, m, play an important role in estimation of the QTL positions in the intervals.

The MIM model is a multiple QTL model and its likelihood is a finite normal mixture. There are two problems that need to be solved for the MIM model. The first is that of parameter estimation of the finite normal mixture model. As m becomes large, the derivation of the maximum-likelihood estimates (MLEs) of the QTL effects and positions in estimation quickly becomes unwieldy. To handle the estimation problem, the general formulas derived by Kao and Zeng (1997) are used to obtain the MLEs in parameter estimation. The second problem is how to find QTL to fit into the MIM model. To select QTL for the MIM model, a stepwise model selection procedure is proposed in strategy of qtl mapping.

LIKELIHOOD OF THE MIM MODEL

In the MIM model, the genotype of each putative QTL, Q_j in interval I_j, is not observed, but its distribution can be inferred from the flanking markers of I_j based on the recombination frequency between them. For every QTL in the backcross population, the conditional probabilities of the QTL genotypes, given different flanking marker genotypes, can be found in Table 1 of Kao and Zeng (1997). To infer the joint conditional probability of the genotype of the m putative QTL, we use the property that if there is no crossing-over interference, the conditional distributions of the individual putative QTL genotypes, given the flanking marker genotypes, are independent. That is, prob(Q1,Q2,⋯,Qm∣I1,I2,⋯,Im)=∏i=1mprob(Qi∣Ii). The joint conditional probability of the m QTL is the product of the marginal conditional probabilities of individual QTL. We refer to p_ij, j = 1, 2, ···, 2^m, as the conditional probabilities of 2^m possible QTL genotypes (note that p_j's denote QTL positions and p_ij's denote the conditional probabilities). If multiple putative QTL within a single marker interval are considered, the individual and joint conditional probabilities of QTL genotypes can be also inferred directly or by a Markov chain procedure (Jiang and Zeng 1997) assuming no interference.

Given a sample with size n, the likelihood function of the MIM model for θ = (p₁, p₂, ···, p_m, a₁, ···, a_m, ···, w_jk, ···, σ²) is L(θ∣Y,X)=∏i=1n[∑j=12mpijϕ(yi−μijσ)], (4) where φ(·) is a standard normal probability density function, μ_ij's correspond to the genotypic values of the 2^m different QTL genotypes in Equation 1, and p_ij's containing information on QTL positions are the corresponding joint conditional probabilities. Statistically, this is a normal mixture model. The density of each individual is a mixture of 2^m possible normal densities with different means μ_ij's and mixing proportions p_ij's. To obtain the MLEs and the asymptotic variance-covariance matrix of the model, the general formulas of Kao and Zeng (1997), based on the expectation and maximization (EM) algorithm (Dempsteret al. 1977), are used for parameter estimation.

PARAMETER ESTIMATION

The likelihood of the MIM model is a finite normal mixture. In parameter estimation, the finite normal mixture model can be treated as an incomplete-data problem (Little and Rubin 1987) by regarding the trait and markers as observed data and the QTL as missing data. The EM algorithm can be used for obtaining the MLEs of the genetic parameters, and Louis's (1982) method can be implemented to obtain the variance-covariance matrix.

In the MIM model, when only one putative QTL (m = 1) is considered in a backcross population, the likelihood is a mixture of two normals (like IM and CIM), and four parameters need to be estimated. The derivation of the MLEs for the one putative QTL model using the EM algorithm has been provided (Carbonellet al. 1992; Zeng 1994). When arbitrary m putative QTL are considered, the likelihood is a mixture of 2^m normals, and at least 2m + 2 parameters (including mean μ, QTL positions and effects, environment variance, and epistasis) need to be estimated. The number of mixture components and parameters increases dramatically as the number of putative QTL taken into account in the model increases. Taking m = 10 as an example, the likelihood is a mixture of 1024 normals with more than 22 parameters to estimate. Therefore, one of the main difficulties with the MIM model is that the derivation of the MLEs quickly becomes unwieldy if the number of putative QTL is large, and an efficient and systematic method for parameter estimation of the MIM model is needed to avoid rederivation for each m. Here, we use the general formulas provided by Kao and Zeng (1997) for deriving the MLEs and the asymptotic variance-covariance matrix of the parameters as the estimation method of MIM. The general formulas are based on two matrices, D and Q. The matrix D is the genetic design matrix that characterizes the genetic parameters of the QTL effects, and the matrix Q is the conditional probability matrix that contains the information on QTL positions. Given the two matrices, the MLEs and the asymptotic variance-covariance matrix can be systematically obtained.

To apply the general formulas to MIM, the genetic design matrix D of the MIM model has the same first m columns as those in Equation 2 for indicating the m main QTL effects and has some or none of the last m(m – 1)/2 columns for specifying epistasis. We refer to D asa2^m × k matrix, where k is the column dimension. There are m individual conditional probability matrices, Q₁, Q₂, ···, and Q_m for the m QTL. The components of the conditional probability matrix Q_j of QTL Q_j in the interval I_j with flanking markers M_j and N_j can be found in Table 1 of Kao and Zeng (1997). For each interval, there are four possible flanking marker genotypes. Totally, there are 4^m possible flanking marker genotypes for m intervals. The joint conditional probability matrix Q then has dimension 4^m × 2^m and can be obtained by Q=Q1⊗Q2⊗⋯⊗Qm , where ⊗ denotes the Kronecker product. The 2^m mixing proportions of any individual i, p_ij's, can be found to be one of the 4^m rows in Q according to its flanking marker genotype. Given the matrices D and Q, the MLEs and the asymptotic variance-covariance matrix can be readily obtained by the general formulas.

Note that, at the tested positions p₁, p₂, ···, and p_m, the mixing proportions p_ij's in the likelihood are fixed and need not be estimated. For obtaining the MLEs of mean, environmental variance, and marginal and epistatic effects, the general equations formulate the iteration of the (t + 1) EM step as follows:

E step: Update the posterior probabilities of the 2^m possible QTL genotypes for each individual i, πij(t+1)=pijϕ((yi−μij(t))∕σ(t))∑j=12mpijϕ((yi−μij(t))∕σ(t));i=1,2,⋯,n,j=1,2,⋯,2m. (5) M step: Find θ^{(t + 1)}, which satisfies the solutions E(t+1)=r(t)−M(t)E(t) (6) μ(t+1)=1n1′[Y−Π(t)DE(t+1)] (7) σ2(t+1)=1n[(Y−1μ(t+1))′(Y−1μ(t+1))−2(Y−1μ(t+1))′Π(t)DE(t+1)+E′(t+1)V(t)E(t+1)], (8) where Π = {π_ij}_n×2m, V = {1′Π(D_i#D_j)}_k×k, r = {(Y – Xβ)′ΠD_i/1′Π(D_i#D_i)}_k×1, and M = {1′Π(D_i#D_j)/1′Π(D_i#D_i) × δ(i ≠ j)}_k×k. D_i(D_j) is the ith(jth) column of the genetic design matrix D. The notation δ(i ≠ j) is an indicator variable that takes value 1 if i ≠ j, and 0 otherwise, and # denotes Hadamard product, which is the element-by-element product of corresponding elements of two same order matrices. For more detailed procedures of the derivation see Kao and Zeng (1997). The E and M steps are iterated until a convergent criterion is satisfied. The converged values are the MLEs. The asymptotic variance-covariance matrix can also be obtained using the general formulas. The general formulas can be easily applied to obtaining the MLEs and evaluating the likelihoods for different genetic models and population structures by setting up the corresponding genetic design matrix D and conditional QTL genotype probability matrices Q_j's. Through comparisons of the likelihoods, hypotheses about the parameters of QTL can be tested by the LRT.

STRATEGY OF QTL MAPPING

For the MIM approach, the second problem that needs to be considered is how to search for QTL to fit into the MIM model. It is quite common that genetic marker data, e.g., rice (Liet al. 1997), pine (Aitkinet al. 1997), and eucalyptus (Grattapagliaet al. 1996), contain more than 100 markers in several linkage groups to cover most of the genome. A QTL is potentially located in any position of each interval. To detect QTL using the MIM model, model selection procedures are considered because all possible subset selection is not feasible. There are at least three basic model selection techniques, forward, backward, and stepwise selections, for exploring the relationship between the independent and dependent variables (Draper and Smith 1981; Kleinbaumet al. 1988; Miller 1990). Several selection criteria, such as Akaike information criterion (AIC; Akaike 1974), cross-validation (Stone 1974), predictive sample reuse (Geisser 1975), Baysian information criterion (BIC; Schwarz 1978), minimum posterior predictive loss (Gelfand and Ghosh 1998), or LRT statistic for selection of variables can be incorporated with model selection techniques to determine the final model. In QTL mapping, any criterion used has to take the genetic marker data structure, such as genome size and distribution of markers, into account. There have been studies on the connection of the LRT statistic to the data structure (see below). So far, however, these related studies lack other criteria. The stepwise selection technique with the LRT statistic as a criterion is adopted for identifying QTL here.

Critical value for claiming QTL detection: When using the LRT statistic as a criterion in model selection for QTL detection, it is very important to determine the appropriate critical value for claiming QTL detection such that correct statistical inference about QTL parameters can be made. Lander and Botstein (1989) suggested using the Bonferroni argument for the sparse-map case and Orenstein-Uhlenbeck diffusion for the dense-map case to determine the critical value. Generally, it has been pointed out that the critical value might need to be adjusted for the number and size of interval, different levels of heritability, different number of multiple linked or unlinked QTL, and linked QTL in the same or opposite direction of effects (Lander and Botstein 1989; Jansen 1993; Zeng 1994). Visscher and Haley (1996) suggested that the critical value should be reduced after a QTL of large effect has been detected. However, most of this information is not available before mapping, and consequently the answers to most of the above questions remain unknown. Churchill and Doerge (1994) therefore suggested using a permutation test for determining an appropriate critical value for specific data sets.

The above considerations on critical value are for the single-QTL model. For a multiple-QTL model, a model selection procedure is required to determine the final model. If stepwise selection is used, the final model is selected from a sequence of nested tests, and the significance level of the sequence will depend on the unknown true model (Atkinson 1980; Terasvirta and Mellin 1986). Therefore, the critical value of the multiple-QTL model depends not only on the above considerations but also on the unknown true model, and the choice of critical value for claiming QTL detection becomes even more complicated for MIM. We are not sure currently what the appropriate critical value is for the MIM model. In practice, the critical value from IM or CIM based on the Bonferroni argument may be used until the complicated issue of choosing the significance level for the multiple-QTL model is solved.

Stepwise selection procedure: The stepwise selection begins with no QTL (m = 0). QTL are then added or deleted one by one in the model. Alternatively, a group of QTL can be added or deleted together. The testing hypotheses for adding or deleting one additional QTL Q_i are H0:ai=0H1:ai≠0, (9) given other, say, k QTL in the model. In hypotheses 9, a_i denotes the effect of Q_i. A LRT statistic LRT=−2logL0L1 is used for testing the hypotheses, where L₀ and L₁ are the likelihoods of the MIM models with k and k + 1 QTL, respectively. If a group of QTL is tested, the hypothesis testing would contain several QTL effects. The stepwise model selection procedure proceeds as follows:

Step 1: Significant values for entry (SVE) and staying (SVS) of a LRT statistic are specified for adding and dropping a QTL in the MIM model. Note that SVE and SVS could be different in model selection.

Step 2: For each position on the genome covered by markers, the LRT statistic reflecting the contribution of the putative QTL to quantitative trait variation is calculated (m = 1; IM). If there are positions with LRT statistics larger than SVE, the position with the largest value will be selected and added first in the model. When m = 1, it is important to note the shape of the likelihood profile and the direction of effect change along the genome for further mapping. Note that quite often no position is found with the LRT statistic larger than SVE when m = 1 because individual QTL contribute little to the trait variation. Two alternative approaches are proposed to prevent the procedure from stopping at a very early stage.

First, when m = 1, the position with the highest LRT statistic is automatically included in the model to initiate the procedure. In our experience, when only one QTL is considered in the model (m = 1), it is quite often found that the LRT statistic of a QTL could be less than SVE. But, when multiple QTL (if any) are accumulated in the model (m > 1), the partial LRT statistics of individual QTL might become significant because more genetic variation is removed from residual variation by taking multiple QTL into account.

Second, chunkwise selection (Kleinbaumet al. 1988) can be used. For closely linked QTL with opposite effects, more than one QTL may be selected in the model as a chunk to effectively reduce the genetic residue in the model. If only one of them is selected, its contribution may not be significant because the effect is canceled out due to failure to consider the others. When m = 1 in the MIM model, the chromosome region with significant change in the directions of effect could suggest that linked QTL with opposite effects are present. Also, epistatic QTL can constitute a chunk. If QTL interact, they may not be significant if only one of them is considered, but they could be significant if they are considered together. Note that the critical value should be higher for chunkwise selection because more parameters are tested. Chunkwise selection allows the incorporation of prior knowledge and preference into the model selection procedure.

Step 3: After the first k QTL are added to the model, the MIM model with m = k + 1 QTL is considered. The position that produces the most significant partial LRT statistic at the SVE level is added into the model. After the k + 1 QTL are fitted to the model, stepwise selection checks all the QTL and deletes any QTL that does not produce a significant partial LRT statistic at the SVS level. Note that a QTL that enters at an early stage may become superfluous at a later stage in stepwise selection procedure. By the same argument, chunkwise selection (m = k + l, l > 1) can be implemented. The stepwise process ends when none of the other positions has a partial LRT statistic significant at the SVE level.

Separating linked QTL: The evidence of multiple-linked QTL clustering in a region could be suggested by the shape of the likelihood profile, for example, a likelihood profile with a wide range of significant multiple peaks, or by significant change in the direction of estimated QTL effects on a chromosome region. To separate closely linked QTL in a certain chromosome region, we can compare the likelihood of the multiple-QTL model with that of a single-QTL model in this region for separation.

Analyses of epistasis: For a backcross population, it can be shown that if epistasis is present and ignored in mapping, the estimates of main effects of epistatic QTL are asymptotically unbiased whether epistasis between QTL is considered in the model or not, and the power of the test for detecting epistatic QTL could be low (appendix). Therefore, when mapping QTL without considering epistasis in a backcross population, the positions and effects of the identified QTL could still be unbiased. For l QTL being tested, there are k = l(l – 1)/2 possible digenic epistases. For each pair of QTL Q_i and Q_j, the hypotheses for testing their epistatic effect w_ij are H0:wij=0H1:wij≠0, (10) given the l QTL in the MIM model. Again, the LRT is used to test the hypotheses. The hypotheses in Equation 10 can also be used to identify QTL with no main effect but interacting with other QTL. To choose the critical value for epistasis detection, a Bonferroni argument can be used. The critical value for rejection of H₀ is suggested as χ1,α∕k2 , where α is the overall significance level.

Fine tuning the estimates of QTL positions and effects: In the above procedures, the estimates of QTL effects and positions were obtained individually. Therefore, the model likelihood might not be at the maximum, and the model is not the final model. To obtain the MLEs of the positions and effects, a multidimensional search around the regions of the identified QTL is suggested. By doing this, QTL estimates can be fine tuned and the final model can be determined. With estimates of QTL positions and effects, other composite genetic parameters (e.g., heritability and variance components) of a quantitative trait can be estimated and response to selection can be predicted.

Construction of the confidence interval for QTL positions and effects: It is important to construct the confidence interval (C.I.) for QTL effects and positions. For example, when a particular QTL is to be transferred to a recipient, a C.I. of QTL position estimate can give us an idea about how large a chromosome segment is around the detected position to be transferred. There are several approaches to constructing a C.I. of the QTL positions and effects, including lod support interval (Lander and Botstein 1989), bootstrapping, using asymptotic standard deviation (ASD; Darvasiet al. 1993; Kao and Zeng 1997), and the methods by Dupuis and Siegmund (1999). Darvasi et al. (1993) and Kao and Zeng (1997) suggested using (p̃ – Z_(1–α/2)S_p̃, p̃ + Z_(1–α/2)S_p̃), where p̃ and S_p̃ are the estimates of QTL position and its standard deviation, to construct a C.I.

Estimation of variance components and heritability: When the final model is determined, the variance components and the heritability of the quantitative trait can be estimated. The ratio V_G/V_p, denoted by hb2 , is called the heritability of a quantitative trait in the broad sense, where V_G and V_p are the genetic and phenotypic variances. The genetic variance V_G can be estimated by the sum of squares of the final model, and the phenotypic variance V_p can be estimated by the total sum of squares. The estimate of hb2 can be approximated by the coefficient of determination R² of the MIM model h^b2=V^GV^p=Model sum of squaresTotal sum of squares=R2.

To estimate the genetic variance components, for example, the total genetic variance contributed by m QTL in the backcross population by Equation 1 is VG=∑i=1mai24+2∑i<jmDijaiaj+∑i<jmδij+∑i<jmδij(116−Dij2)wij2, (11) where D_ij is the gametic linkage disequilibrium coefficient between Q_i and Q_j (Weir 1996). The coefficient D_ij is equivalent to (1 – 2r_ij)/4, where r_ij is the recombination fraction between two QTL. In Equation 11, the first term is the genetic variance contributed by QTL Q_i. The second term, D_ija_ia_j, is the genetic covariance between two QTL due to linkage disequilibrium. The last term is contributed by epistasis. The genetic variance component contributed by Q_i is defined by ai2∕4 . However, the estimated genetic component by a^i2∕4 is biased, and this bias can be corrected by [a^i2+Var(a^i)]∕4 . The genetic covariance between Q_i and Q_j is defined by 2Da_ia_j. By the same argument, the estimated genetic covariance by 2D̂â_iâ_j is also biased and can be corrected by 2D̂ [â_iâ_j + Cov(â_i, â_j)] under the assumption that the effect and location of QTL are independent. Other genetic components can also be estimated in the same way. For an F₂ population or a backcross population with segregation distortion, the partition of genetic variance into components is presented by C-H. Kao and Z-B. Zeng (unpublished results).

Estimation of individual genotypic value and marker-assisted selection: In plant or animal breeding, individuals with high genotypic values or favorable genotypes are usually selected for producing progeny. With the estimated QTL effects and positions, the genotypic values of individuals can be estimated by Equation 1 and the favorable QTL genotypes can be determined for selection. To select individuals with large trait values, genotype AA (Aa) of nonepistatic QTL with positive (negative) effects is preferred. For QTL with epistasis, their epistatic effects must be considered in selecting the best combination of genotypes. If QTL controlling different traits are closely linked or at the same positions, traits are genetically correlated. Selecting individuals for improvement of one trait will affect the other trait due to linkage or pleiotropy. In practice, selecting individuals with the desired character for one trait will frequently accompany an undesired character for other traits. By considering circumstances such as genetic correlation between traits, the distances between markers and QTL, and the effects of QTL, the best strategies of marker-assisted selection for (multiple) trait improvement under specific purposes and requirements can be explored.

DATA ANALYSIS

Radiata pine: Radiata pine is one of the most widely planted forestry species in the Southern Hemisphere. Two elite parents were crossed to produce 134 progeny. For each progeny, random amplified polymorphic DNA (RAPD) markers were generated, and traits measured included annual brown cone number at eight years of age, diameter of stem at breast height, and branch quality score. The cone number per tree, which varied from 0 to 45, was transformed to approximate a normal distribution using a square root transformation. The quality of branches of a tree were scored on a scale from 1 (poorest) to 6 (best). The mean of several branch quality scores denoted the branch quality of a tree. A pseudotestcross strategy is used to construct a linkage map for each parent, and then a backcross model can be used for mapping QTL for each parent separately (Grattapaglia and Sederoff 1994; Grattapagliaet al. 1996). The analysis reported here is on one parent. A genetic marker map was constructed using MapMaker/EXP (Lincolnet al. 1993). The RAPD marker data contained 120 markers in 12 linkage groups and covered ∼1679.3 cM. The average spacing of the 107 marker intervals was 13.5 cM.

As mentioned in strategy of qtl mapping, the choice of critical value is a very complicated issue for the multiple-QTL model. The value depends on the marker data structure and several unknown QTL parameters (true model). In data analysis, a critical value from IM based on Bonferroni argument is used to evaluate and illustrate the MIM approach. The SVE and SVS of the LRT statistic for claiming a QTL detection at the overall α = 0.05 level were chosen as 12.12 (χ1,0.05∕1072≈χ1,0.00052) . For QTL selected as a chunk, the overall α = 0.05 level was chosen as χk,0.05∕1072 , where k is the number of tested parameters in the chunk.

QTL detection: For trait DBH, when m = 1, there is no position along the genome with an LRT statistic higher than SVE. The position with the largest LRT statistic (7.85; R² = 0.0639) was found at position [12,5,0] (0 cM away from the left marker of the fifth marker interval on the twelfth linkage group). The chromosome region between C1M3 (the third marker of the first linkage group) and C1M7 showed opposite direction of effects. At C1M3, the effect was positive (P = 0.57), while at C1M4 and C1M5, the effects were negative (P = 0.0253 and 0.4181, respectively). The genetic distance between C1M3 and C1M4 is 74.8 cM. It could suggest that there are two closely linked QTL with opposite directions of effects in this region. If only one QTL (m = 1) is fitted in the model for search, the effect can be canceled out by opposing QTL effects. QTL will be out of detection as shown by the LRT statistic profile of IM in Figure 1. Therefore, on linkage group 1, the MIM model with m = 2 selected two candidate QTL, at positions [1,3,63] and [1,4,0], as a chunk. The partial LRT statistic for fitting the two QTL in the model was 13.13 (SVE and SVS for two parameters are χ2,0.00052=15.2 ), and the model ^R2 was 0.2104. Although the LRT statistic was less than SVE, the two QTL were selected as a chunk to initiate the stepwise selection process.

The procedure restarted at m = 2 by fitting two QTL with effects of opposite directions at [1,3,63] and [1,4,0]. The partial LRT statistics were 8.034 and 8.458 for the two QTL, with estimated effects 65.65 and –73.48, respectively. Given QTL at [1,3,63] and [1,4,0] in the model, a QTL at [10,5,12] with partial LRT statistic 12.83 was selected into the model (m = 3). The partial LRT statistics became 14.89, 15.42, and 12.83, which were all larger than the SVS of 12.12, for the three QTL. The model R² was 0.3202. Given these three QTL in the model, the largest partial LRT statistic 7.40 was found at position [2,2,0]. A chunkwise selection for epistatic QTL was attempted. If the candidate QTL at [2,2,0] and [12,5,12] with epistasis were selected as a chunk (m = 5 and one epistasis, k = 6), the partial LRT statistic of the chunk would be 24.76 (compared with χ3,0.00052=17.73 ). The partial LRT statistics were 23.48, 24.39, and 8.76 for the three preselected QTL at [1,3,63], [1,4,0], and [10,5,12], respectively. The QTL at [10,5,12] became nonsignificant and, therefore, was dropped from the model. Given the four QTL [1,3,63], [1,4,0], [2,2,0], and [12,5,12] in the MIM model, no other single position had a partial LRT statistic >8.76. The chunkwise selection was implemented again to find epistatic QTL. When the candidate QTL at [5,5,0] and [10,5,12] with epistasis were considered as the third chunk, the partial LRT statistic was 19.85. Adding these two epistatic QTL into the model (m = 6 and two epistasis, k = 8), the partial LRT statistics were 19.48, 20.69, and 26.91 for QTL at positions [1,3,63], [1,4,0], and the first chunk of QTL, respectively. Given the six QTL, no other QTL were identified.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

—Results of QTL mapping. (a–e) The solid triangles denote the QTL positions localized by MIM. The size of triangle reflects the size of QTL effect.

Fine tuning the estimates of QTL position and effect: Two epistatic pairs were identified as described above; no other epistatic interaction between QTL was found. No QTL without main effect but interacting with the identified QTL were found. A multidimensional search around the detected QTL was used to fine tune the estimates of QTL parameters. The locations changed to [1,3,61], [1,4,0], [2,2,0], [5,5,0], [10,5,9], and [12,5,9]. The estimated QTL effects are shown in Table 1. QTL at positions [1,3,61], [2,2,0], [10,5,9], and [5,5,0] had positive effects, and QTL at positions [1,4,0] and [12,5,9] had negative effects. The effects of QTL at positions [1,3,61] and [1,4,0] were larger when compared with others. The model R² was 0.5226. Therefore, six identified QTL were conclusively identified in QTL mapping for the diameter trait. The partial LRT statistic profiles for each QTL are shown in Figure 1.

Epistasis: The estimated epistatic effect between QTL at positions [2,2,0] and [12,5,9] was 39.54 (partial LRT statistic 15.23), and the epistatic effect between QTL at [5,5,0] and [10,5,9] was 22.64 (partial LRT statistic 4.84). Figure 2 shows how the QTL interact. Figure 2a shows that the effect of QTL (G_BB – G_Bb) at position [12,5,9] was positive in the background of homozygote QTL (AA) at position [2,2,0], but it was negative in the heterozygote background (Aa). Figure 2b shows that the QTL at position [10,5,9] had a large effect in the background of homozygote QTL (AA) at [5,5,0], but it had a small effect in the heterozygote background (Aa).

Heritability and variance components: The broad sense heritability for tree diameter can be estimated by the R² value of the final MIM model. The R² of the model including six QTL and two epistases was 0.5226. QTL at positions [2,2,0], [5,5,0], [10,5,12], and [12,5,9] contributed ∼4.50, 1.36, 5.25, and 1.76% of the total genetic variance, respectively. The percentage of genetic variance contributed by the two linked QTL separated by 13.8 cM on the first linkage group was 76.75%. There was a negative genetic covariance between the two linked QTL. Two epistatic pairs contributed ∼10.38% to the total genetic variance.

QTL mapping for cone number and branch quality: QTL mapping was also performed on the traits of cone number and branch score. The mapping results are listed in Table 1. For cone number, seven QTL were identified (although the QTL at [1,1,3] was not significant with partial LRT statistic 9.44, we considered it as a candidate QTL). Epistasis was found between two QTL pairs using chunkwise selection. The model R² value of the MIM model fitted to the seven QTL and their epistasis was 0.5606. The two linked QTL, separated by 27.6 cM on linkage group 10, contributed 29.93% of the genetic variance. The other five QTL contributed ∼55.93% of the total genetic variance. Epistasis contributed 14.14%. For branch quality, five QTL were identified (we also considered the two QTL with partial LRT statistic values 10.37 and 10.36 at [1,4,11] and [12,5,0] as candidate QTL). No epistasis was found for QTL controlling branch score. The model R² was 0.3630. Two linked QTL, separated by 19.6 cM on linkage group 11, contributed 48.69% of the genetic variance. The remaining three QTL contributed from ∼11 to 27% of the total genetic variance.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

—Epistasis between QTL controlling tree diameter. (a) Epistasis between QTL at positions [2,2,0] and [12,5,9]. A and B denote QTL at [2,2,0] and [12,5,9], respectively. (b) Epistasis between QTL at positions [5,5,0] and [10,5,9]. A and B denote QTL at [5,5,0] and [10,5,9], respectively.

Confidence intervals of QTL positions and effects: The lod support interval and the ASD of QTL effect and position are listed in Table 1. Out of the 18 QTL detected for three traits, 9 ([2,6,0], [5,10,0], and [10,9,0] for cone number; [1,4,0], [2,2,0], and [5,5,0] for tree diameter; [2,1,0], [11,6,0], and [12,5,0] for branch score) of them were localized at the markers, and 2 ([10,5,9] and [12,5,9]) had negative ASD. Therefore, the ASD of these QTL position estimates were not available for constructing C.I.'s. The asymmetric lod support intervals are typical in this case. For example, the diameter QTL at [5,5,0] has an asymmetric lod support interval ([5,4,7], [5,5,16]). In general, the interval constructed by ASD is much narrower than the lod support interval. For example, C.I.'s constructed using four times ASD were 6.52 and 7.04 cM for the cone QTL at [6,4,18] and [12,3,2], and the lod support intervals are 59.6 cM and 14.6 cM, respectively.

Marker-assisted selection: Individuals with favorable QTL genotypes are selected as parents to produce progeny. Trees carrying all the favorable QTL genotypes were not found for each trait in the sample. Therefore, only a subset of the detected QTL was considered in selection. For tree diameter, three trees were found to carry favorable genotypes and two trees were found to carry unfavorable genotypes (consider epistasis) of the five QTL (out of the six detected QTL) at positions [1,4,0], [2,2,0], [5,5,0], [10,5,9], and [12,5,9]. The observed trait means for the two groups were 232.38 and 163.05 mm, respectively, through selection of these five diameter QTL. The estimated genotypic values of the two groups were 233.84 and 160.06 mm (Table 2). The observed and estimated values of performing selection for the other two traits on the sample based on four and five QTL are also shown in Table 2. The mapping results in Table 1 also allow us to estimate the genotypic values of certain genotypes. For example, if trees carrying all six favorable diameter QTL were selected with epistasis taken into consideration, the estimated tree diameter for those trees would be 314.17 mm and the estimated cone number would be 8.22. If trees carry all seven favorable QTL (epistasis considered) for reducing cone number, the estimated cone number would be 0.33 and the estimated tree diameter would be 196.45 mm. Consequently, the improvement of tree diameter would cause simultaneous increase in cone number, which is a reflection of the positive genetic correlation between the two traits. Generally, the estimated and observed results were quite close based on the MIM result as found in this sample.