Email updates

Keep up to date with the latest news and content from BMC Genetics and BioMed Central.

Open Access Methodology article

A new mapping method for quantitative trait loci of silkworm

Hai-Ming Xu1, Chang-Shuai Wei1, Yun-Ting Tang2, Zhi-Hong Zhu1, Yang-Fu Sima3 and Xiang-Yang Lou14*

Author Affiliations

1 Institute of Bioinformatics, College of Agriculture and Biotechnology, Zhejiang University, Hangzhou 310029, China

2 Ningbo Institute of Technology, Zhejiang University, Ningbo 315100, China

3 College of Agriculture Science and Technology, Soochow University, Suzhou 215006, China

4 Department of Biostatistics, School of Public Health, University of Alabama at Birmingham, Birmingham, Alabama 35294, USA

For all author emails, please log on.

BMC Genetics 2011, 12:19  doi:10.1186/1471-2156-12-19


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2156/12/19


Received:5 September 2010
Accepted:28 January 2011
Published:28 January 2011

© 2011 Xu et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 F2 population of silkworm. The number and positions of QTLs were determined by F-test 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 F2 model. We also found that a sample size of hundreds was sufficiently large to unbiasedly estimate all the four types of epistases (i.e., additive-additive, additive-dominance, dominance-additive, and dominance-dominance) 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 follow-up 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 co-segregation 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 linkage-maps 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) [2-6]. 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 QTL-environment interactions are usually involved in the genetic variation of complex trait [9-11], several complicated mapping models were developed to analyze epistatic effects [12-18]. 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,20-24]. Recently, there are several new Bayesian methods developed for mapping QTLs underlying dynamic traits [25], ordinal traits [26], multiple traits [27-30], 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,34-37], 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 F2 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 F2 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 F2 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 F2 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 sex-specific 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,..., th; 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]:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/19/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/19/mathml/M1">View MathML</a>

(1)

where μ is the population mean; ak and dk are the additive effect and dominant effect of the k-th QTL, respectively; xAik and xDik 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 rM+Q,rM-Q), respectively; aai, adi, dai and ddi are the additive-additive, additive-dominance, dominance-additive and dominance-dominance epistatic effects of the l-th pair of QTLs, respectively, with their coefficients xAAil, xADil, xADil and xDDil, which are the products of the corresponding xA and xD; all the above additive, dominance and epistatic effects are of our interest and thus considered as fixed; Sh is the sex effect of sex h, Sh ~ (0, σ2S); ashk and dshk are additive-sex and dominance-sex interaction effects, ashk ~ (0, σ2AS) and dshk ~ (0, σ2DS), respectively; aashl, adshl, dashl and ddshl are interaction effects between epistasis and sex, aashl ~ (0, σ2AAS), adshl ~ (0, σ2ADS), dashl ~ (0, σ2DAS) and ddshl ~ (0, σ2DDS); εhi is the random residual effect, εhi ~ (0, σ2e).

The above model can be expressed in the matrix form:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/19/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/19/mathml/M2">View MathML</a>

(2)

where y is the vector of the phenotypic values; b is the vector of fixed effects and X is the coefficient matrix; eu is the vector of the u-th random effect and Uu is the corresponding coefficient matrix, and Ru is an identity matrix for every u in this model if the components of eu 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 xA, xD, 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 F2 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 genome-wide scanning for QTLs underlying silkworm traits. Pairs of adjacent markers are selected [39] and their effects are tested in the following model:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/19/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/19/mathml/M6">View MathML</a>

(3)

where a+th (a-th) and d+th(d-th) are the additive and dominance effects due to the right (left) marker of the t-th marker interval in the h-th 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 transitional-possibility-matrix algorithm will be employed to calculate their expected values. To determine which pair of adjacent markers should be selected, the F-testing is applied and the threshold value is determined by permutations [40]. After all the marker intervals exceeding the F-critical 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, genome-wide one-dimensional 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,

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/19/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/19/mathml/M7">View MathML</a>

(4)

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 F-statistic:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/19/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/19/mathml/M8">View MathML</a>

(5)

where Q denotes QTL genetic effects with coefficient matrix of XQ, and M denotes maker effects with coefficient matrix of XM, XQM is a matrix catenated by the XQ and XM; n is the number of observed values, rank (XQM) and rank (XM) are the ranks of matrix XQM and XM, respectively; SSR(Q|M) 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(Q|M) and SSE can be calculated using Henderson III method [41]. The permutation technique is used to determine the critical value of F-test. 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, two-dimensional 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:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/19/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/19/mathml/M9">View MathML</a>

(6)

where aa+ph(aa-ph), ad+ph(ad-ph), da+ph(da-ph) and dd+ph(dd-ph) denote the additive-additive, additive-dominance, dominance-additive and dominance-dominance epistatic effects within the h-th sex between two right (left) markers of the p-th 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 F-statistic to test its extra effects is calculated using the formula (5), and the critical value to declare significance is specified by calculating F-statistic 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 two-dimensional searching strategy. For the l-th 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 one-dimensional scanning,

<a onClick="popup('http://www.biomedcentral.com/1471-2156/12/19/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/12/19/mathml/M10">View MathML</a>

(7)

where mp is the number of selected interval pairs, mi is the number of QTLs identified in one-dimensional searching; all other parameters are defined the same as those in model (1) and model (6). Similar F-test 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 mixed-model 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 F2 population are substantially different from those in the normal F2 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 F2 and the traditional chiasmata F2 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 r1 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 (Di = |xai - xci|) and relative difference (Ri = |xai - xci |/xai) for the i-th flanking marker genotype (7 totally), xa (xc) is the coefficient of QTL effect in the achiasmate (chiasmate) model; and (4) to investigate the maximum and the minimum of the set of Di and Ri (i = 1, 2,..., 7) for every r1.

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 (Q1-Q7), E2 (Q3-Q5), E3 (Q6-Q7) were involved in epistatic effects while no additive effects, dominance effects, additive-sex interaction effects or dominance-sex 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 F2 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 genome-wide search for QTLs. The first one used the proposed model (1) (the silkworm F2 model), called Model I hereafter, the second used the traditional chiasmate F2 model (i.e., all coefficients in model (1) are replaced with those determined according to the genetic structure of a normal chiasmate F2 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 F2 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 F2 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 r1, r2 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.

thumbnailFigure 1. Difference of coefficients in the achiamata and chiamata models. x-axis stands for the recombination rate between QTL and left flanking marker and y-axis 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 r1, 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 F2 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 F2 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 z-statistic value of -2.063 and a probability of 0.020 by Wilcoxon two-sample 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 F2 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 F2 population in this article, we could obtain the unbiased prediction of random effects including additive-sex, dominance-sex and epistasis-sex 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.

Table 4. Estimation of epistasis-sex interaction effects with Model I a

Table 5. Estimation of QTL-sex interaction effects with Model I and Model III a

Discussion

Crossing-over 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 crossing-over 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 F2 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 F2 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 F2 population, every individual receives one gamete from female parent without crossing-over and another one generated by potential sister-chromatid exchange from male parent, resulting in a different population structure from the normal chiasmate F2 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: sex-limited, sex-influenced (also known as sex-controlled), and sex-linked; 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 [46-48]. 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 non-sex chromosomes. Sex chromosomes usually play a unique role in many biological processes and phenomena, including sex determination, epigenetic chromosome-wide regulation of gene expression [50]. Sex chromosomes have many different genetic features compared with autosomes and there is extra complexity in mapping of sex-linked 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 sex-biased 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 chromosome-specific significance threshold is required. Given the complexity of sex-linked inheritance, tailored mapping methods are needed to effectively hunt sex-linked 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 family-based association test to detect QTLs on X-chromosome 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 F2 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 F2 population. One alternative choice is excluding the higher-order genetic effects of QTL (additive-dominance, dominance-additive, dominance-dominance 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 F2 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 F-statistic 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 2009ZX08009-004B, the National Natural Science Foundation 30571405, and the National Institutes of Health Grant DA025095.

References

  1. 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):1937-1940. PubMed Abstract | Publisher Full Text OpenURL

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

    Curr Sci India 2002, 83(4):415-425. OpenURL

  3. 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):1277-1284. PubMed Abstract | PubMed Central Full Text OpenURL

  4. Yamamoto K, Nohata J, Kadono-Okuda K, Narukawa J, Sasanuma M, Sasanuma S, Minami H, Shimomura M, Suetsugu Y, Banno Y, et al.: A BAC-based integrated linkage map of the silkworm Bombyx mori.

    Genome Biol 2008., 9(1) PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  5. Yasukochi Y: A dense genetic map of the silkworm, Bombyx mori, covering all chromosomes based on 1018 molecular markers.

    Genetics 1998, 150(4):1513-1525. PubMed Abstract | PubMed Central Full Text OpenURL

  6. 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 OpenURL

  7. Lander ES, Botstein D: Mapping Mendelian Factors Underlying Quantitative Traits Using Rflp Linkage Maps.

    Genetics 1989, 121(1):185-199. PubMed Abstract | PubMed Central Full Text OpenURL

  8. Zeng ZB: Precision Mapping of Quantitative Trait Loci.

    Genetics 1994, 136(4):1457-1468. PubMed Abstract | PubMed Central Full Text OpenURL

  9. 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):46-52. PubMed Abstract | Publisher Full Text OpenURL

  10. Symington LS: Role of RAD52 epistasis group genes in homologous recombination and double-strand break repair.

    Microbiol Mol Biol R 2002, 66(4):630-+. Publisher Full Text OpenURL

  11. Tabanao DA, Yu J, Bernardo R: Multilocus epistasis, linkage, and genetic variance in breeding populations with few parents.

    Theor Appl Genet 2007, 115(3):335-342. PubMed Abstract | Publisher Full Text OpenURL

  12. Kao CH, Zeng ZB, Teasdale RD: Multiple interval mapping for quantitative trait loci.

    Genetics 1999, 152(3):1203-1216. PubMed Abstract | PubMed Central Full Text OpenURL

  13. Ljungberg K, Holmgren S, Carlborg O: Simultaneous search for multiple QTL using the global optimization algorithm DIRECT.

    Bioinformatics 2004, 20(12):1887-1895. PubMed Abstract | Publisher Full Text OpenURL

  14. Piepho HP: A mixed-model approach to mapping quantitative trait loci in barley on the basis of multiple environment data.

    Genetics 2000, 156(4):2043-2050. PubMed Abstract | PubMed Central Full Text OpenURL

  15. 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(7-8):1255-1264. Publisher Full Text OpenURL

  16. Yang J, Zhu J, Williams RW: Mapping the genetic architecture of complex traits in experimental populations.

    Bioinformatics 2007, 23(12):1527-1536. PubMed Abstract | Publisher Full Text OpenURL

  17. 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):1865-1877. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. Lou XY, Chen GB, Yan L, Ma JZ, Zhu J, Elston RC, Li MD: A generalized combinatorial approach for detecting gene-by-gene and gene-by-environment interactions with application to nicotine dependence.

    Am J Hum Genet 2007, 80(6):1125-1137. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  19. 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):325-333. PubMed Abstract | Publisher Full Text OpenURL

  20. 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):805-816. PubMed Abstract | PubMed Central Full Text OpenURL

  21. Yi NJ, Yandell BS, Churchill GA, Allison DB, Eisen EJ, Pomp D: Bayesian model selection for genome-wide epistatic quantitative trait loci analysis.

    Genetics 2005, 170(3):1333-1344. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  22. 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):429-437. PubMed Abstract | Publisher Full Text OpenURL

  23. Durrant C, Mott R: Bayesian Quantitative Trait Locus Mapping Using Inferred Haplotypes.

    Genetics 2010, 184(3):839-U375. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  24. Fang M: Bayesian shrinkage mapping of quantitative trait loci in variance component models.

    Bmc Genet 2010., 11 BioMed Central Full Text OpenURL

  25. Yang RQ, Xu SZ: Bayesian shrinkage analysis of quantitative trait loci for dynamic traits.

    Genetics 2007, 176(2):1169-1185. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  26. Yi NJ, Banerjee S, Pomp D, Yandell BS: Bayesian mapping of genomewide interacting quantitative trait loci for ordinal traits.

    Genetics 2007, 176(3):1855-1864. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  27. 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):304-320. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  28. Xu CW, Wang XF, Li ZK, Xu SZ: Mapping QTL for multiple traits using Bayesian statistics.

    Genetics Research 2009, 91(1):23-37. PubMed Abstract | Publisher Full Text OpenURL

  29. Yi NJ, Banerjee S: Hierarchical Generalized Linear Models for Multiple Quantitative Trait Locus Mapping.

    Genetics 2009, 181(3):1101-1113. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  30. Zou F, Huang HW, Lee S, Hoeschele I: Nonparametric Bayesian Variable Selection With Applications to Multiple Quantitative Trait Loci Mapping With Epistasis and Gene-Environment Interaction.

    Genetics 2010, 186(1):385-U600. PubMed Abstract | Publisher Full Text OpenURL

  31. 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 OpenURL

  32. 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 OpenURL

  33. 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):367-375. PubMed Abstract | Publisher Full Text OpenURL

  34. 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):477-484. Publisher Full Text OpenURL

  35. 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 OpenURL

  36. Lie Z, Cheng L, Fang-yin D, Shou-min F: Mapping of major quantitative trait loci for economic traits of silkworm cocoon.

    Genet Mol Res 2010, 9(1):78-88. PubMed Abstract | Publisher Full Text OpenURL

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

    Afr J Biotechnol 2010, 9(10):1427-1433. OpenURL

  38. Jiang CJ, Zeng ZB: Mapping quantitative trait loci with dominant and missing markers in various crosses from two inbred lines.

    Genetica 1997, 101(1):47-58. PubMed Abstract | Publisher Full Text OpenURL

  39. Piepho HP, Gauch HG: Marker pair selection for mapping quantitative trait loci.

    Genetics 2001, 157(1):433-444. PubMed Abstract | PubMed Central Full Text OpenURL

  40. Doerge RW, Churchill GA: Permutation tests for multiple loci affecting a quantitative character.

    Genetics 1996, 142(1):285-294. PubMed Abstract | PubMed Central Full Text OpenURL

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

  42. 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. OpenURL

  43. Wang CS, Rutledge JJ, Gianola D: Bayesian-Analysis of Mixed Linear-Models Via Gibbs Sampling with an Application to Litter Size in Iberian Pigs.

    Genet Sel Evol 1994, 26(2):91-115. BioMed Central Full Text OpenURL

  44. Promboon A, Shimada T, Fujiwara H, Kobayashi M: Linkage Map of Random Amplified Polymorphic Dnas (Rapds) in the Silkworm, Bombyx-Mori.

    Genet Res 1995, 66(1):1-7. Publisher Full Text OpenURL

  45. Yu HR, Edderkaoui B, Cortez A, Davidson H, Wergedal J, Baylink D, Mohan S: Mapping of the chromosome 17 BMD QTL in the F-2 male mice of MRL/MpJ × SJL/J.

    Genetica 2009, 135(1):59-66. PubMed Abstract | Publisher Full Text OpenURL

  46. Abasht B, Pittel F, Lagarrigue S, Le Bihan-Duval 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):297-311. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  47. Farber CR, Medrano JF: Fine mapping reveals sex bias in quantitative trait loci affecting growth, skeletal size and obesity-related traits on mouse chromosomes 2 and 11.

    Genetics 2007, 175(1):349-360. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  48. 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 HcB-8 and HcB-23 recombinant mouse strains.

    Bone 2010, 46(5):1251-1259. PubMed Abstract | Publisher Full Text OpenURL

  49. Zhu J, Weir BS: Diallel analysis for sex-linked and maternal effects.

    Theor Appl Genet 1996, 92(1):1-9. Publisher Full Text OpenURL

  50. Kaiser VB, Bachtrog D: Evolution of sex chromosomes in insects.

    Annu Rev Genet 2010, 44:91-112. PubMed Abstract | Publisher Full Text OpenURL

  51. Broman KW, Sen S, Owens SE, Manichaikul A, Southard-Smith EM, Churchill GA: The X chromosome in quantitative trait locus mapping.

    Genetics 2006, 174(4):2151-2158. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  52. Zhang L, Martin ER, Morris RW, Li YJ: Association Test for X-Linked QTL in Family-Based Designs.

    Am J Hum Genet 2009, 84(4):431-444. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  53. Lukens LN, Doebley J: Epistatic and environmental interactions for quantitative trait loci involved in maize evolution.

    Genet Res 1999, 74(3):291-302. Publisher Full Text OpenURL

  54. 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):4656-4660. Publisher Full Text OpenURL

  55. Gurganus MC, Nuzhdin SV, Leips JW, Mackay TFC: High-resolution mapping of quantitative trait loci for sternopleural bristle number in Drosophila melanogaster.

    Genetics 1999, 152(4):1585-1604. PubMed Abstract | PubMed Central Full Text OpenURL

  56. Montooth KL, Marden JH, Clark AG: Mapping determinants of variation in energy metabolism, respiration and flight in Drosophila.

    Genetics 2003, 165(2):623-635. PubMed Abstract | PubMed Central Full Text OpenURL