Abstract
Background
Silkworm is the basis of sericultural industry and the model organism in insect genetics study. Mapping quantitative trait loci (QTLs) underlying economically important traits of silkworm is of high significance for promoting the silkworm molecular breeding and advancing our knowledge on genetic architecture of the Lepidoptera. Yet, the currently used mapping methods are not well suitable for silkworm, because of ignoring the recombination difference in meiosis between two sexes.
Results
A mixed linear model including QTL main effects, epistatic effects, and QTL × sex interaction effects was proposed for mapping QTLs in an F_{2 }population of silkworm. The number and positions of QTLs were determined by Ftest and model selection. The Markov chain Monte Carlo (MCMC) algorithm was employed to estimate and test genetic effects of QTLs and QTL × sex interaction effects. The effectiveness of the model and statistical method was validated by a series of simulations. The results indicate that when markers are distributed sparsely on chromosomes, our method will substantially improve estimation accuracy as compared to the normal chiasmate F_{2 }model. We also found that a sample size of hundreds was sufficiently large to unbiasedly estimate all the four types of epistases (i.e., additiveadditive, additivedominance, dominanceadditive, and dominancedominance) when the paired QTLs reside on different chromosomes in silkworm.
Conclusion
The proposed method could accurately estimate not only the additive, dominance and digenic epistatic effects but also their interaction effects with sex, correcting the potential bias and precision loss in the current QTL mapping practice of silkworm and thus representing an important addition to the arsenal of QTL mapping tools.
Background
Silkworm (Bombyx mori) is the basis of sericultural industry. With nearly 5000 years' domestication, silkworm has an undoubted importance in human history and is still of great value in modern economy. In addition, it is also an ideal model organism of the Lepidoptera. Because silkworm is easy to rear and could produce large amount of mutation, it is second to fruit fly as a model organism in insect genetics study. Over these years, the "old" creature is becoming a new hot spot in genetic research.
Many important traits of silkworm are complex quantitative traits, such as whole cocoon weight and cocoon shell weight, etc. The genetic variation of quantitative traits are usually controlled by a number of genes (quantitative trait loci, QTLs) with epistatic and gene × environment interactions. To locate their positions on chromosome and estimate their contribution to the variation of trait is a key step for positional cloning and followup utilization of those genes. With the development of modern molecular biology, it has become possible to dissect the genetic mechanism of quantitative trait and to identify the associated genes and their interacting network by cosegregation analysis of molecular markers in a mapping population based on specific genetic model connecting QTL genotype with a phenotype of interest.
Since the draft sequence of silkworm genome [1] was reported, genetic research of silkworm has been greatly spurred and many linkagemaps have been constructed with various molecular makers, such as random amplified polymorphic DNAs (RAPDs), amplified fragment length polymorphisms (AFLPs), selective amplification of DNA fragments (SADFs), microsatellites (also known as simple sequence repeats, SSRs) and expressed sequence tags (ESTs) [26]. Meanwhile, many statistical models were developed for mapping QTLs, such as the interval mapping (IM) method [7] and the composite interval mapping (CIM) method [8]. Along with increasing evidence supporting the claim that epistatic and QTLenvironment interactions are usually involved in the genetic variation of complex trait [911], several complicated mapping models were developed to analyze epistatic effects [1218]. Kao et al. [12] expanded CIM to multiple interval mapping (MIM) for detecting epistasis. Wang et al. [15] established a mixed linear model based composite interval mapping (MCIM) method to analyze both epistasis and QE interaction in a double haploid (DH) population. A few years later, the MCIM was extended for other designed populations [19] and improved in searching strategy and estimating QTL parameters [16]. Parallel to these frequentist methods, Bayesian approaches have also been proposed for QTL mapping of complex traits [17,2024]. Recently, there are several new Bayesian methods developed for mapping QTLs underlying dynamic traits [25], ordinal traits [26], multiple traits [2730], expression data [31], and gene frequency data [32], and for multiple inbreed lines designs [33].
Although there are a large number of QTLs reported in other species, relatively fewer QTL mapping studies are performed in silkworm [6,3437], which in part result from the fact that the current models and corresponding analysis methods are not appropriate for the silkworm. Silkworm has a particular characteristic called female achiasmata where meiosis occurs with no crossover between homologous chromosome pairs. Yet the majority of recently developed mapping tools are based on the assumption of chiasmata without considering the genetic differences between male and female. Therefore, only specific mapping populations such as a backcross (BC) population are suggested for gene mapping in order to satisfy this assumption [3,5]. By setting female silkworm as the homozygous parental lines, they can avoid the problem of achiasmata. However, a BC population contains fewer segregant types of molecular markers than F_{2 }population. As a result, genetic information is not enough to reveal additive and dominance effects, epistatic effects of QTLs and their interaction effects with environment. Therefore, it is necessary to develop a new method of QTL mapping with consideration of female achiasmata for F_{2 }population of silkworm.
In the present study, we proposed a new statistical method for QTL mapping in silkworm. The method can analyze the additive, dominance and digenic epistic effects of QTLs, as well as their interaction with sex. The effectiveness of the method is investigated by intensive Monte Carlo simulations.
Methods
Genetic model for QTL mapping of silkworm
To specify the achiasmate characteristic of silkworm, an F_{2 }population is used to illustrate our methods. Models for other mapping populations or experimental designs can be established by including or excluding relevant QTL main effects or QTL × environment interaction effects. Without loss of generality, we assume here that the phenotypic value of a quantitative trait in a F_{2 }population is controlled by additive and dominant effects of n QTLs and m digenic epistatic effects of QTLs. Because there are differences in both the genetic material and meiosis mode between male and female silkworms, sex effects should be involved in the genetic variation of silkworm traits (e.g. cocoon quality trait). As sex functions as an endogenous environment for the development of an autosomal sexspecific trait, sex is treated as a pseudo environmental covariate in our experimental design. When necessary, the interaction effects between QTLs and sex can be included in our model. Therefore, the phenotypic value of individual i with sex h (i = 1, 2,..., t_{h}; h = 1, 2) can be expressed as the following mixed linear model which is extended from those in Gao and Zhu [19] and in Yang et al. [16]:
where μ is the population mean; a_{k }and d_{k }are the additive effect and dominant effect of the kth QTL, respectively; x_{Aik }and x_{Dik }are the coefficients of QTL effects which can, when QTL genotypes are unobserved, be derived from the conditional probability of the putative QTL given the genotypes of flanking markers (flanking marker M_{+}, M_{ }of QTL Q) and the QTL position (the recombination frequency r_{M+Q},r_{MQ}), respectively; aa_{i}, ad_{i}, da_{i }and dd_{i }are the additiveadditive, additivedominance, dominanceadditive and dominancedominance epistatic effects of the lth pair of QTLs, respectively, with their coefficients x_{AAil}, x_{ADil}, x_{ADil }and x_{DDil}, which are the products of the corresponding x_{A }and x_{D}; all the above additive, dominance and epistatic effects are of our interest and thus considered as fixed; S_{h }is the sex effect of sex h, S_{h }~ (0, σ^{2}_{S}); as_{hk }and ds_{hk }are additivesex and dominancesex interaction effects, as_{hk }~ (0, σ^{2}_{AS}) and ds_{hk }~ (0, σ^{2}_{DS}), respectively; aas_{hl}, ads_{hl}, das_{hl }and dds_{hl }are interaction effects between epistasis and sex, aas_{hl }~ (0, σ^{2}_{AAS}), ads_{hl }~ (0, σ^{2}_{ADS}), das_{hl }~ (0, σ^{2}_{DAS}) and dds_{hl }~ (0, σ^{2}_{DDS}); ε_{hi }is the random residual effect, ε_{hi }~ (0, σ^{2}_{e}).
The above model can be expressed in the matrix form:
where y is the vector of the phenotypic values; b is the vector of fixed effects and X is the coefficient matrix; e_{u }is the vector of the uth random effect and U_{u }is the corresponding coefficient matrix, and R_{u }is an identity matrix for every u in this model if the components of e_{u }are independent.
Since the QTL genotype of each individual is unobservable in the real experiment, the coefficients of QTL effects are unknown. However, they could be substituted with their expectation based on the conditional probability of QTL genotype given the genotypes of flanking markers. The expected coefficients of additive and dominance effects given the flanking markers are listed in Table 1, while the coefficients of epistatic effects and QTL × sex interaction effects can be calculated by multiplying the corresponding x_{A}, x_{D}, and coefficients of sex effects (0 or 1). When the information on a portion of markers is missing, the algorithm based on transitional possibility matrix, proposed by Jiang and Zeng [38], can be applied to impute missing data. As shown in Table 1, the major difference between our achiasmate model and the traditional chiasmate model lies in calculation of the expected coefficients.
Table 1. Conditional probabilities of QTL genotypes and coefficients of additive and dominance effects for achiasmate and chiasmate F_{2 }populations
The above QTL full model, assuming the number of QTLs and their positions are known, can be used to detect the significance of QTL effects. Based on the final QTL full model after excluding insignificant QTLs, all genetic effects of QTL and QTL × sex interaction effects will be estimated by the mixed linear model approach.
Scanning for QTLs with main effects
To remedy a potential bias in both the estimated effect values and position arisen from linked QTLs and control the residual background variation, Zeng [8] developed CIM method by integrating IM with multiple regression. Zeng [8] showed that, conditional on an intermediate marker, its two flanking markers would be independent of each other in a backcross population (it also holds true in a double haploid population), assuming no crossover interference, providing the theoretical basis of Zeng's separation of multiple linked QTLs. To avoid "ghost" QTL due to the impact of linked QTLs and control the background genetic variance, as in Zeng [8], we like first to select significant markers as background markers prior to genomewide scanning for QTLs underlying silkworm traits. Pairs of adjacent markers are selected [39] and their effects are tested in the following model:
where a^{+}_{th }(a^{}_{th}) and d^{+}_{th}(d^{}_{th}) are the additive and dominance effects due to the right (left) marker of the tth marker interval in the hth sex, with corresponding coefficients ζ^{+}_{A}(ζ^{}_{A}) and ζ^{+}_{D}(ζ^{}_{D}); the other parameters have the same definition as those in model (1). ζ^{+}_{A}(ζ^{}_{A}) takes value of 1, 0, 1 when the marker genotype is MM, Mm and mm, respectively. ζ^{+}_{D}(ζ^{}_{D}) takes value of 0.5 for homozygous genotype (MM and mm) and 0.5 for heterozygote (Mm). If the marker genotype information is missing, the transitionalpossibilitymatrix algorithm will be employed to calculate their expected values. To determine which pair of adjacent markers should be selected, the Ftesting is applied and the threshold value is determined by permutations [40]. After all the marker intervals exceeding the Fcritical value are included into the model, stepwise model selection method is performed to eliminate all ghost peaks.
Once the marker intervals with significant effects are identified, genomewide onedimensional searching for QTLs can be conducted with inclusion of selected markers in the model to control the background effects from other unknown QTLs. The following model will be used to test a putative QTL k,
where mi represents the number of selected marker intervals; the other parameters are the same as defined in the models (1) and (3). The scanning is performed with a walk step, say 1 cM, within every selected marker interval. The significance of the putative QTL is tested by the following Fstatistic:
where Q denotes QTL genetic effects with coefficient matrix of X_{Q}, and M denotes maker effects with coefficient matrix of X_{M}, X_{QM }is a matrix catenated by the X_{Q }and X_{M}; n is the number of observed values, rank (X_{QM}) and rank (X_{M}) are the ranks of matrix X_{QM }and X_{M}, respectively; SSR(QM) is the extra sum of squares due to the genetic effects of the putative QTL given the inclusion of M in the model; SSE is the residual sum of squares. SSR(QM) and SSE can be calculated using Henderson III method [41]. The permutation technique is used to determine the critical value of Ftest. For all the QTLs detected to be significant at the level of 0.05, the stepwise selection is conducted to eliminate the false positive peaks.
Scanning for paired QTLs with epistatic effects
In order to detect the paired QTLs with significant interaction effects, twodimensional whole genome scanning strategy should be adopted, while, the QTL model also need to be extended to inclusion of epistatic effects of paired QTLs. Before scanning epistatic effects, we still perform marker selection procedure. All the possible marker interval pairs are tested in the following model:
where aa^{+}_{ph}(aa^{}_{ph}), ad^{+}_{ph}(ad^{}_{ph}), da^{+}_{ph}(da^{}_{ph}) and dd^{+}_{ph}(dd^{}_{ph}) denote the additiveadditive, additivedominance, dominanceadditive and dominancedominance epistatic effects within the hth sex between two right (left) markers of the pth marker interval pairs, respectively; the coefficients of epistatic effects can be calculated by the products of coefficients of marker major effects in model (3); other parameters are defined the same as those in model (4). For each paired marker intervals tested, the Fstatistic to test its extra effects is calculated using the formula (5), and the critical value to declare significance is specified by calculating Fstatistic in a series of randomly shuffling observation vector y s. All paired intervals above the critical value are then picked up as significant candidate interactions.
Suppose there are mp pairs of marker intervals selected. Within two paired marker intervals, the epistatic effects from two paired putative QTLs is tested in twodimensional searching strategy. For the lth paired putative QTLs, the following mapping model can be analyzed with inclusion of the epistatic effects of the other selected marker intervals and the main effects of the QTLs detected in onedimensional scanning,
where mp is the number of selected interval pairs, mi is the number of QTLs identified in onedimensional searching; all other parameters are defined the same as those in model (1) and model (6). Similar Ftest and selection procedure are applied.
Estimation of QTL parameters in the full model
After the number and positions of QTLs are specified, a full model consisting of all genetic effects of QTLs and their interaction effects with environment (sex) is established. The variance components of random effects can be estimated by restricted maximum likelihood (REML), the fixed effects by generalized least squares (GLS) or ordinary least squares (OLS), and the random effects by adjusted unbiased prediction (AUP). These mixedmodel estimates of QTL effects are set to be the initial values of MCMC methods [42]. The sample distributions of QTL parameters are obtained by the Gibbs sampling [16,43]. Finally, each effect is estimated by the distribution mean, while significance of an effect is tested by t statistic.
Numerical calculation of the difference in the coefficients of QTL effects between the achiasmate and the chiasmate models
As there is no genetic material exchange for female silkworm when producing gamete in meiosis, the marker frequency distribution and the conditional probabilities of QTL genotypes in silkworm F_{2 }population are substantially different from those in the normal F_{2 }population with chiasma. To demonstrate the difference in QTL detection and investigate the inappropriateness of QTL mapping of silkworm traits based on the traditional genetic model, we compared the conditional probabilities of QTL genotypes given the flanking marker genotypes under the achiasmata F_{2 }and the traditional chiasmata F_{2 }models that were calculated from Table 1.
We used two cases where r equal to 0.09 (10 cM) and 0.16 (20 cM), respectively, to evaluate the difference in the additive and dominance coefficients between two models summarized in the following steps: (1) to set r fixed and r_{1 }changed from 0 to r; (2) to calculate each coefficient in two models based on the QTL conditional probability in Table 1 for every given marker genotype; (3) to calculate the absolute difference (D_{i }= x_{ai } x_{ci}) and relative difference (R_{i }= x_{ai } x_{ci }/x_{ai}) for the ith flanking marker genotype (7 totally), x_{a }(x_{c}) is the coefficient of QTL effect in the achiasmate (chiasmate) model; and (4) to investigate the maximum and the minimum of the set of D_{i }and R_{i }(i = 1, 2,..., 7) for every r_{1}.
Simulation scenarios
To investigate the efficiency and accuracy of the proposed methods, 300 simulations were conducted with the following QTL configuration. 5 chromosomes and 7 QTLs were considered. Each chromosome had 11 molecular markers and 10 equal marker intervals of 10 cM; 7 QTLs (Q1, Q2, Q3, Q4, Q5, Q6, Q7) were scattered on 5 chromosomes (see Table 2 and Table 3 for details), wherein, three pairs of QTLs E1 (Q1Q7), E2 (Q3Q5), E3 (Q6Q7) were involved in epistatic effects while no additive effects, dominance effects, additivesex interaction effects or dominancesex interaction effects were set for Q6 and Q7 (Tables 2 and 5). Detailed information about the positions and effects is presented in Tables 2, 3, 4 and 5. The simulations were performed based on an F_{2 }population with 300 individuals in which the numbers of males and females are equal. The proportions of total variation due to genetic effects and genetic × environment interaction effects were 50.47% and 19.53%, respectively; the narrow heritability was 33.56%. Three QTLs (Q2, Q3, Q4) and all epistatic pairs were considered to interact with sex having a variance ranging from 4 to 6.24, while Q1 and Q5 had no interaction with sex.
Table 2. Estimation of QTL positions and main effects with Models I, II and III ^{a}
Table 3. Estimation of the positions and epistatic effects of the paired QTLs with Models I and II ^{a}
In simulations, three different strategies were employed to conduct a genomewide search for QTLs. The first one used the proposed model (1) (the silkworm F_{2 }model), called Model I hereafter, the second used the traditional chiasmate F_{2 }model (i.e., all coefficients in model (1) are replaced with those determined according to the genetic structure of a normal chiasmate F_{2 }population), called Model II hereafter, and the third used the reduced version of model (1) where all epistatic effects and interaction effects of epistasis with sex were excluded, called Model III hereafter. The above three strategies were used to analyze the same simulated data sets generated by the silkworm F_{2 }model with QTL effects, QTL × sex interaction effects.
We first examined the performance of the newly proposed strategy (Model I) and the traditional strategy (Model II) in mapping for silkworm traits and demonstrated the potential bias and loss of power caused by Model II. These simulation results were summarized in Tables 2 and 3. As the role of epistasis in the genetic control of complex traits has been well recognized, the comparison between Model I and Model III could offer us insight into epistasis detection. These results are listed in Tables 2 and 5.
Results
Comparison of coefficients in models for the achiasmate and the normal chiasmate F_{2 }populations
As shown by the formula of conditional probability in Table 1, we could see that each probability value in achiasmate model is approximately equal to a first order function of r_{1}, r_{2 }or r, while, each one in chiasmate model approximates to a second order function of recombination rate, suggesting that there should be considerable difference between the two models, which can potentially result in estimation bias and loss in accuracy.
The numerical examples could also illustrate this point. In the setting of r = 0.09, the absolute difference of coefficient ranged from 0 to 0.0025 for additive effect and from 0.0001 to 0.0050 for dominance effect (Figure 1A, respectively, while it did from 0 to 0.0090 for additive effect and from 0.0002 to 0.0170 for dominance effect under the setting of r = 0.16 (Figure 1B). The relative absolute difference varied in range of 0.0005 to 0.0093 for additive effect and in range of 0.0009 to 0.0145 for dominance effect when r = 0.09 (Figure 1C), while it did in range of 0.0018 to 0.0344 and in range of 0.0036 to 0.0526 when r = 0.16 (Figure 1D). It was also very clear that the maximum or minimum difference was reached when the QTL was at the middle of the marker interval.
Figure 1. Difference of coefficients in the achiamata and chiamata models. xaxis stands for the recombination rate between QTL and left flanking marker and yaxis stands for the absolute value of coefficient difference or its percentage. Max curve is formed by the maximum absolute difference (or relative difference) for every different r_{1}, while min curve by the minimum. Plot A and C were drawn at the marker interval of 10 cM (r = 0.09), and Plot B and D at the marker interval of 20 cM (r = 0.16).
According to the above comparison, it could be concluded that the traditional chiasmate F_{2 }model would lead to a biased estimation when it was applied to mapping QTLs underlying silkworm traits and our proposed method would improve QTL mapping accuracy.
Comparison between the models for the achiasmate silkworm and the normal chiasmate F_{2 }populations
Model I and the traditional Model II were used to analyze the simulated data, and were compared for their abilities in estimating the position and genetic effects of QTLs. The results suggested that Model I had better estimation accuracy (relatively smaller bias and standard deviation) in QTL position and effects than Model II (Table 2). All bias of genetic main effects from Model I were less than 5% of the true values, whereas Model II sometimes gave a larger bias, e.g., the bias of the Q5 additive effect > 10% of the true value (Table 2). Model I had a considerably larger power to detect the five QTLs (ranging from 62.7 ~ 98.3) than Model II (ranging from 49.7 ~ 94.7) (Table 2), regardless of whether the QTL has interaction effects with sex (Q2, Q3, Q4) or not (Q1, Q5) (Table 5). For all QTLs detected in the 300 simulations, the false discovery rate of Model I is 6.17%, prominently smaller than that of Model II (9.56%) with a zstatistic value of 2.063 and a probability of 0.020 by Wilcoxon twosample test. Model I also provided an estimation of QTL position closer to the true value and had smaller standard deviation (SD) than the model II did (Table 2).
We also compared the estimation accuracies of epistatic effects from Model I with those from Model II (Table 3). Although both the Model I and the Model II could estimate all epistatic effects reasonably well, Model I had relative smaller SD and greater power than Model II in detecting the three pairs of QTLs involved in epistatic interactions, wherein E1 stands for the interaction between one QTL with main effects and another without main effects, E2 stands for the interaction between both QTLs with main effects and E3 stands for the interaction between two QTLs without main effects. As for estimation of QTL position, Model I also outperformed Model II in accuracy.
Comparison of the silkworm F_{2 }models with and without epistasis
For the estimation of QTL positions, it could be found there was a slight difference in the estimated values between two models when the detected QTL has additive or dominance effects, or their interaction effects with sex, as well as in the corresponding SDs (Table 2). Two models had the same power in detecting QTL (Table 2), although Model I included more parameters of QTL than Model III. As for the paired QTLs with purely epistatic effects, they couldn't be identified by Model III, but could be detected by Model I. Such a result is expected, considering that (1) Model I and Model III used the same marker genotype information and quantitative trait values, (2) one dimensional scanning procedure did not include epistasis in Model III, and (3) the position of QTL was determined by the result of one dimensional scanning.
Although there was a small difference in position estimation, the estimates of additive and dominance effects were apparently different between two models (Table 2). Most of the estimated values in Model I were closer to the true values as compared with those in Model III. We could also see that, in Model I every estimate of additive or dominance effects had smaller SD than that in Model III, suggesting that Model I could provide more stable and more unbiased estimation (Table 2). As for QTL × sex interaction, although Model III could also give relatively accurate estimation, it was clear that the results of Model I were better than those of Model III, in terms of biasedness or SD (see Table 5).
Prediction of QTL × sex interaction effects
Based on the QTL model proposed for silkworm F_{2 }population in this article, we could obtain the unbiased prediction of random effects including additivesex, dominancesex and epistasissex interaction effects by the MCMC method. In the simulation studies, Model I provided not only the estimation of the genetic effects of QTLs including epistasis but also their interaction effects with sex (shown in Tables 4 and 5). For two kinds of QTLs involved in interaction with sex (Q2, Q3, Q4) or not (Q1, Q5), the predicted values of additive × sex and dominance × sex interaction effects were very close to the true realized values of random effects, and the corresponding SDs were acceptably low (Table 5). For the interaction between epistasis and sex (Table 4), all four types of interaction effects could be well predicted by Model I. Although slightly bigger bias and SDs were observed, compared with the results for additive or dominance × sex interaction effects, they were still acceptable.
Discussion
Crossingover is an important issue in organisms with meiosis, which can not only enrich phenotypic variation among individuals but also speed up the process of evolution. However, there are many exceptions such as Drosophila and silkworm, which are not only of great value as model organism for biology study but also important for agriculture. Drosophila and silkworm have a common characteristic in reproduction that crossingover occurs only in one sex, while there is a difference that such a phenomenon occurs in female parent for silkworm and in male parent for Drosophila. The particular action in meiosis, also called achiasmata, needs to be considered in the process of gene mapping. In order to avoid the problem on the achiasmate and use available genetic model and software for a normal chiasmate mapping population, investigators have proposed to conduct QTL analysis based on BC population of silkworm. However, there are obvious disadvantages in this solution: recombination information and diversity of genotype in the BC population are not so rich as the F_{2 }population for unraveling genetic architecture of complex traits where complex epistasis and gene × sex or environment interaction are involved. The existing studies on constructing linkage map or genetic mapping with F_{2 }populations for silkworm [5,44] chose to simply neglect such a recombination difference between two sexes because of lack of appropriate analytical method. This potentially leads to bias and loss in precision as, in silkworm F_{2 }population, every individual receives one gamete from female parent without crossingover and another one generated by potential sisterchromatid exchange from male parent, resulting in a different population structure from the normal chiasmate F_{2 }population. Thus, our proposed method that can accommodate achiasmata phenomenon and also effectively handle epistatic effects and the interaction effects of QTL and environment, represents a necessary addition to the current toolkit of QTL mapping.
Many of the widely used statistical methods and software, such as IM and CIM, do not include sex effects in the models because they are mainly designed for plants such as Arabidopsis thaliana and rice. But in animals, many quantitative traits are sex dependent and behave much differently between males and females, such as the cocoon traits of silkworm. Sex specific traits can be categorized into three types of inheritance: sexlimited, sexinfluenced (also known as sexcontrolled), and sexlinked; the former two of which are controlled by autosomal gene(s) and sex, in our term, in which there are sex and/or gene × sex interaction effects, and the latter one is caused by gene(s) carried on the sex chromosomes. It has been well documented that there are sex differences in terms of the presence/absence and locations of QTLs [45], as well as the interaction of QTL with sex [4648]. Thus, for the purpose of improving analysis power, the sex effect and the QTL × sex interaction effect should be generally included into the analysis model as a covariate to eliminate the influence from sex. In our study, as in some literature [49], the sex effect is considered in the QTL model as a random effect for the purpose of background control. Simulation results showed the proposed method could improve both statistical power to detect a QTL and estimation accuracy for genetic effects of QTL and QTL × sex interaction. We like to point out that the sex and the sex related interaction effects can also be treated as fixed ones in our model when necessary.
The method presented here is mainly to detect QTLs on nonsex chromosomes. Sex chromosomes usually play a unique role in many biological processes and phenomena, including sex determination, epigenetic chromosomewide regulation of gene expression [50]. Sex chromosomes have many different genetic features compared with autosomes and there is extra complexity in mapping of sexlinked genetically inherited traits. First, there are two categories of sex determination systems: heterogametic male (XY) and heterogametic female (ZW), and the heterogametic sex is hemizygous in which gene dosage effect or dosage compensation mechanism may occur. Second, sex chromosomes can show sexbiased transmission. Third, there may also exist random inactivation of the sex chromosome. Broman et al. [51] addressed that if the sex chromosome is treated like an autosome, a sex difference in the phenotype can lead to spurious linkage on the sex chromosome. Further, the number of degrees of freedom for the linkage test may be different for the sex chromosome than for autosomes, and so sex chromosomespecific significance threshold is required. Given the complexity of sexlinked inheritance, tailored mapping methods are needed to effectively hunt sexlinked genes. Thus, Broman et al. [51] proposed a method for mapping QTL on X chromosome in experimental crosses population. Zhang et al. [52] developed a familybased association test to detect QTLs on Xchromosome under consideration of the dosage effect due to female X chromosome inactivation. It is possible to extend the proposed method to mapping QTLs on the sex chromosomes.
The genetic variation in continuous traits is usually governed by a polygenic network system, composed of many genes with a small effect, and sometimes one or a few genes of large effect. Recently, intensive studies on quantitative variation have pointed to that epistasis is usually involved in genetic variation of quantitative traits. Strong interactions between QTLs have been observed in maize [53], soybean [54] and Drosophila [55]. In addition, QTLs with minor or no individual effect can also be involved in epistatic interaction [56]. More and more attention has been paid to molecular dissection of epistasis. Our proposed model includes not only the digenic epistatic effects, but also their interaction effects with sex. Therefore, this model can well tackle the complexity of quantitative trait in silkworm. Simulations revealed that the proposed method could present better estimation of QTL parameters no matter whether or not the epistasis and/or their interaction with sex exists.
Lastly but not least, it should be pointed out that only seven different genotypes of two QTL loci on the same chromosome can be generated in F_{2 }population because of the female achiasmate of silkworm [2], while, in the full model, eight genetic effects of a pair of QTLs (two additive effects, two dominance effects and four epistatic effects) need to be estimated. Therefore, under this situation, the proposed method could not produce unbiased estimate of all eight fixed effects. An elaborately planned design is required to effectively detect epistases between interacting loci located on the same chromosome due to the insufficient number of segregating genotypes in an achiasmate F_{2 }population. One alternative choice is excluding the higherorder genetic effects of QTL (additivedominance, dominanceadditive, dominancedominance epistatic effects) from the model. However, for the case in which two QTLs are residing on two different chromosomes, there are still nine QTL genotypes segregated in F_{2 }population of silkworm since the chromosome of female parent could be passed independently to its progeny. It explained why Model I could well estimate all epistatic effects in our simulation study. Furthermore, fortunately, the position of QTL could be estimated unbiasedly no matter whether the QTLs are distributed on the same chromosome or not, since the QTL position is distinguished based on the Fstatistic measuring the total extra effects due to tested variables in the model which is not affected by the correlation between these variables.
Conclusion
We have developed a genetic model for mapping QTL in silkworm F2 population which could analyze the additive effect, dominance effect, digenic epistatic effect and their interaction effects with sex, and correct the potential bias and precision loss in the current QTL mapping practice of silkworm, thus representing an important addition to the arsenal of QTL mapping tools.
Authors' contributions
HMX and CSW developed methods, performed simulation analyses, interpreted results and drafted the manuscript. YTT, ZHZ and YFS helped performed simulation studies. XYL directed the development, design, simulations, interpretation of results, and drafted manuscript. All authors read and approved the manuscript.
Acknowledgements
This work was supported in part by the National Basic Research Program of China (973 Program) 2007CB109007, the National Special Program for Breeding New Transgenic Variety 2009ZX08009004B, the National Natural Science Foundation 30571405, and the National Institutes of Health Grant DA025095.
References

Xia QY, Zhou ZY, Lu C, Cheng DJ, Dai FY, Li B, Zhao P, Zha XF, Cheng TC, Chai CL, et al.: A draft sequence for the genome of the domesticated silkworm (Bombyx mori).
Science 2004, 306(5703):19371940. PubMed Abstract  Publisher Full Text

Nagaraju J, Goldsmith MR: Silkworm genomics  progress and prospects.

Tan YD, Wan CL, Zhu YF, Lu C, Xiang ZH, Deng HW: An amplified fragment length polymorphism map of the silkworm.
Genetics 2001, 157(3):12771284. PubMed Abstract  PubMed Central Full Text

Yamamoto K, Nohata J, KadonoOkuda K, Narukawa J, Sasanuma M, Sasanuma S, Minami H, Shimomura M, Suetsugu Y, Banno Y, et al.: A BACbased integrated linkage map of the silkworm Bombyx mori.
Genome Biol 2008., 9(1) PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Yasukochi Y: A dense genetic map of the silkworm, Bombyx mori, covering all chromosomes based on 1018 molecular markers.
Genetics 1998, 150(4):15131525. PubMed Abstract  PubMed Central Full Text

Zhan SA, Huang JH, Guo QH, Zhao YP, Li WH, Miao XX, Goldsmith MR, Li MW, Huang YP: An integrated genetic linkage map for silkworms with three parental combinations and its application to the mapping of single genes and QTL.
Bmc Genomics 2009., 10 BioMed Central Full Text

Lander ES, Botstein D: Mapping Mendelian Factors Underlying Quantitative Traits Using Rflp Linkage Maps.
Genetics 1989, 121(1):185199. PubMed Abstract  PubMed Central Full Text

Zeng ZB: Precision Mapping of Quantitative Trait Loci.
Genetics 1994, 136(4):14571468. PubMed Abstract  PubMed Central Full Text

Mao YC, London NR, Ma L, Dvorkin D, Da Y: Detection of SNP epistasis effects of quantitative traits using an extended Kempthorne model.
Physiol Genomics 2006, 28(1):4652. PubMed Abstract  Publisher Full Text

Symington LS: Role of RAD52 epistasis group genes in homologous recombination and doublestrand break repair.
Microbiol Mol Biol R 2002, 66(4):630+. Publisher Full Text

Tabanao DA, Yu J, Bernardo R: Multilocus epistasis, linkage, and genetic variance in breeding populations with few parents.
Theor Appl Genet 2007, 115(3):335342. PubMed Abstract  Publisher Full Text

Kao CH, Zeng ZB, Teasdale RD: Multiple interval mapping for quantitative trait loci.
Genetics 1999, 152(3):12031216. PubMed Abstract  PubMed Central Full Text

Ljungberg K, Holmgren S, Carlborg O: Simultaneous search for multiple QTL using the global optimization algorithm DIRECT.
Bioinformatics 2004, 20(12):18871895. PubMed Abstract  Publisher Full Text

Piepho HP: A mixedmodel approach to mapping quantitative trait loci in barley on the basis of multiple environment data.
Genetics 2000, 156(4):20432050. PubMed Abstract  PubMed Central Full Text

Wang DL, Zhu J, Li ZK, Paterson AH: Mapping QTLs with epistatic effects and QTLxenvironment interactions by mixed linear model approaches.
Theor Appl Genet 1999, 99(78):12551264. Publisher Full Text

Yang J, Zhu J, Williams RW: Mapping the genetic architecture of complex traits in experimental populations.
Bioinformatics 2007, 23(12):15271536. PubMed Abstract  Publisher Full Text

Yi NJ, Shriner D, Banerjee S, Mehta T, Pomp D, Yandell BS: An efficient Bayesian model selection approach for interacting quantitative trait loci models with many effects.
Genetics 2007, 176(3):18651877. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lou XY, Chen GB, Yan L, Ma JZ, Zhu J, Elston RC, Li MD: A generalized combinatorial approach for detecting genebygene and genebyenvironment interactions with application to nicotine dependence.
Am J Hum Genet 2007, 80(6):11251137. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gao YM, Zhu J: Mapping QTLs with digenic epistasis under multiple environments and predicting heterosis based on QTL effects.
Theor Appl Genet 2007, 115(3):325333. PubMed Abstract  Publisher Full Text

Satagopan JM, Yandell YS, Newton MA, Osborn TC: A Bayesian approach to detect quantitative trait loci using Markov chain Monte Carlo.
Genetics 1996, 144(2):805816. PubMed Abstract  PubMed Central Full Text

Yi NJ, Yandell BS, Churchill GA, Allison DB, Eisen EJ, Pomp D: Bayesian model selection for genomewide epistatic quantitative trait loci analysis.
Genetics 2005, 170(3):13331344. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fang M, Jiang D, Gao HJ, Sun DX, Yang RQ, Zhang Q: A new Bayesian automatic model selection approach for mapping quantitative trait loci under variance component model.
Genetica 2009, 135(3):429437. PubMed Abstract  Publisher Full Text

Durrant C, Mott R: Bayesian Quantitative Trait Locus Mapping Using Inferred Haplotypes.
Genetics 2010, 184(3):839U375. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fang M: Bayesian shrinkage mapping of quantitative trait loci in variance component models.
Bmc Genet 2010., 11 BioMed Central Full Text

Yang RQ, Xu SZ: Bayesian shrinkage analysis of quantitative trait loci for dynamic traits.
Genetics 2007, 176(2):11691185. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yi NJ, Banerjee S, Pomp D, Yandell BS: Bayesian mapping of genomewide interacting quantitative trait loci for ordinal traits.
Genetics 2007, 176(3):18551864. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Liu JF, Liu YJ, Liu XG, Deng HW: Bayesian mapping of quantitative trait loci for multiple complex traits with the use of variance components.
Am J Hum Genet 2007, 81(2):304320. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Xu CW, Wang XF, Li ZK, Xu SZ: Mapping QTL for multiple traits using Bayesian statistics.
Genetics Research 2009, 91(1):2337. PubMed Abstract  Publisher Full Text

Yi NJ, Banerjee S: Hierarchical Generalized Linear Models for Multiple Quantitative Trait Locus Mapping.
Genetics 2009, 181(3):11011113. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zou F, Huang HW, Lee S, Hoeschele I: Nonparametric Bayesian Variable Selection With Applications to Multiple Quantitative Trait Loci Mapping With Epistasis and GeneEnvironment Interaction.
Genetics 2010, 186(1):385U600. PubMed Abstract  Publisher Full Text

Zhang W, Zhu J, Schadt EE, Liu JS: A Bayesian Partition Method for Detecting Pleiotropic and Epistatic eQTL Modules.
Plos Comput Biol 2010., 6(1) PubMed Abstract  Publisher Full Text  PubMed Central Full Text

He W, Fernando RL, Dekkers JCM, Gilbert H: A gene frequency model for QTL mapping using Bayesian inference.
Genet Sel Evol 2010., 42 PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Fang M, Jiang D, Chen X, Pu LJ, Liu SC: Bayesian analysis of genetic architecture of quantitative trait using data of crosses of multiple inbred lines.
Genetica 2008, 134(3):367375. PubMed Abstract  Publisher Full Text

Lu C, Li B, Zhao AC, Xiang ZH: QTL mapping of economically important traits in silkworm (Bombyx mori).
Sci China Ser C 2004, 47(5):477484. Publisher Full Text

Mirhoseini SZ, Rabiei B, Potki P, Dalirsefat SB: Amplified fragment length polymorphism mapping of quantitative trait loci for economically important traits in the silkworm, Bombyx mori.
J Insect Sci 2010., 10 PubMed Abstract  Publisher Full Text

Lie Z, Cheng L, Fangyin D, Shoumin F: Mapping of major quantitative trait loci for economic traits of silkworm cocoon.
Genet Mol Res 2010, 9(1):7888. PubMed Abstract  Publisher Full Text

Mirhosseini SZ, Bizhannia AR, Rabiei B, Taeb M, Siedavi AR: Identification of AFLP markers linked with cocoon weight genes in silkworm (Bombyx mori L.).

Jiang CJ, Zeng ZB: Mapping quantitative trait loci with dominant and missing markers in various crosses from two inbred lines.
Genetica 1997, 101(1):4758. PubMed Abstract  Publisher Full Text

Piepho HP, Gauch HG: Marker pair selection for mapping quantitative trait loci.
Genetics 2001, 157(1):433444. PubMed Abstract  PubMed Central Full Text

Doerge RW, Churchill GA: Permutation tests for multiple loci affecting a quantitative character.
Genetics 1996, 142(1):285294. PubMed Abstract  PubMed Central Full Text

Searle SR, Casella G, Mcculloch C: Variance components. New York: A WileyInterscience Publication, Hohn Wiley & Sons, INC; 1992.

Macedo FWM, Gianola D: Bayesian analysis of univariate mixed models with informative priors.
European Association for Animal Production, 38th Annual Meeting, Lisbon, Portugal 1987, 35.

Wang CS, Rutledge JJ, Gianola D: BayesianAnalysis of Mixed LinearModels Via Gibbs Sampling with an Application to Litter Size in Iberian Pigs.
Genet Sel Evol 1994, 26(2):91115. BioMed Central Full Text

Promboon A, Shimada T, Fujiwara H, Kobayashi M: Linkage Map of Random Amplified Polymorphic Dnas (Rapds) in the Silkworm, BombyxMori.
Genet Res 1995, 66(1):17. Publisher Full Text

Yu HR, Edderkaoui B, Cortez A, Davidson H, Wergedal J, Baylink D, Mohan S: Mapping of the chromosome 17 BMD QTL in the F2 male mice of MRL/MpJ × SJL/J.
Genetica 2009, 135(1):5966. PubMed Abstract  Publisher Full Text

Abasht B, Pittel F, Lagarrigue S, Le BihanDuval E, Le Roy P, Demeure O, Vignoles F, Simon J, Cogburn L, Aggrey S, et al.: Fatness QTL on chicken chromosome 5 and interaction with sex.
Genet Sel Evol 2006, 38(3):297311. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Farber CR, Medrano JF: Fine mapping reveals sex bias in quantitative trait loci affecting growth, skeletal size and obesityrelated traits on mouse chromosomes 2 and 11.
Genetics 2007, 175(1):349360. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Saless N, Franco GEL, Litscher S, Kattappuram RS, Houlihan MJ, Vanderby R, Demant P, Blank RD: Linkage mapping of femoral material properties in a reciprocal intercross of HcB8 and HcB23 recombinant mouse strains.
Bone 2010, 46(5):12511259. PubMed Abstract  Publisher Full Text

Zhu J, Weir BS: Diallel analysis for sexlinked and maternal effects.
Theor Appl Genet 1996, 92(1):19. Publisher Full Text

Kaiser VB, Bachtrog D: Evolution of sex chromosomes in insects.
Annu Rev Genet 2010, 44:91112. PubMed Abstract  Publisher Full Text

Broman KW, Sen S, Owens SE, Manichaikul A, SouthardSmith EM, Churchill GA: The X chromosome in quantitative trait locus mapping.
Genetics 2006, 174(4):21512158. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zhang L, Martin ER, Morris RW, Li YJ: Association Test for XLinked QTL in FamilyBased Designs.
Am J Hum Genet 2009, 84(4):431444. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lukens LN, Doebley J: Epistatic and environmental interactions for quantitative trait loci involved in maize evolution.
Genet Res 1999, 74(3):291302. Publisher Full Text

Lark KG, Chase K, Adler F, Mansur LM, Orf JH: Interactions between Quantitative Trait Loci in Soybean in Which Trait Variation at One Locus Is Conditional Upon a Specific Allele at Another.
P Natl Acad Sci USA 1995, 92(10):46564660. Publisher Full Text

Gurganus MC, Nuzhdin SV, Leips JW, Mackay TFC: Highresolution mapping of quantitative trait loci for sternopleural bristle number in Drosophila melanogaster.
Genetics 1999, 152(4):15851604. PubMed Abstract  PubMed Central Full Text

Montooth KL, Marden JH, Clark AG: Mapping determinants of variation in energy metabolism, respiration and flight in Drosophila.
Genetics 2003, 165(2):623635. PubMed Abstract  PubMed Central Full Text