Fine mapping of multiple interacting quantitative trait loci using combined linkage disequilibrium and linkage information

Fine mapping of multiple interacting quantitative trait loci using combined linkage disequilibrium and linkage information first-author †E-mail:slee38@une.edu.au †E-mail:slee38@une.edu.au slee38@une.edu.au School of Rural Science and Agriculture, University of New England, Armidale, New South Wales 2351, Australia School of Rural Science and Agriculture, University of New England, Armidale, New South Wales 2351, Australia Quantitative trait loci (QTL) and their additive, dominance and epistatic effects play a critical role in complex trait variation. It is often infeasible to detect multiple interacting QTL due to main effects often being confounded by interaction effects. Positioning interacting QTL within a small region is even more difficult. We present a variance component approach nested in an empirical Bayesian method, which simultaneously takes into account additive, dominance and epistatic effects due to multiple interacting QTL. The covariance structure used in the variance component approach is based on combined linkage disequilibrium and linkage (LDL) information. In a simulation study where there are complex epistatic interactions between QTL, it is possible to simultaneously fine map interacting QTL using the proposed approach. The present method combined with LDL information can efficiently detect QTL and their dominance and epistatic effects, making it possible to simultaneously fine map main and epistatic QTL. Quantitative trait loci (QTL) Combined linkage disequilibrium and linkage (LDL) information Epistatic effects INTRODUCTION 1 In general, there are multiple quantitative trait loci (QTL), and their intra- (dominance) and inter-locus interaction (epistasis) underlying phenotypes of a complex trait. If such gene interactions are ignored in a QTL analysis, estimated QTL positions and effects will not be accurate and precise. For considering complicate epistatic interactions in QTL analysis, variance component approaches based on linkage information have been widely used in natural and outbred populations (Mitchell et al., 1997 1997 ; Blangero et al., 2000 2000 ; 2001 2001 ; Purcell and Sham, 2004 2004 ). However, the use of linkage information alone may limit the power to detect dominance and epistasis, and mapping resolution based on linkage information is not high. Linkage disequilibrium (LD) in addition to linkage information (Meuwissen and Goddard, 2001 2001 ) can give useful extra information about additive, dominance and epistatic effects for small genomic regions. Moreover, a variance component approach nested in an empirical Bayesian method makes it possible to simultaneously map multiple QTL as shown in Lee and van der Werf ( 2006 2006 ). Therefore, simultaneous fine mapping of multiple epistatic QTL within a small region would be possible. The aim of this study is to investigate how much the mapping resolution improves when considering QTL interactions in fine mapping of a complex trait. The posterior QTL density is estimated using an empirical Bayesian approach based on combined linkage disequilibrium and linkage (LDL) information (Lee and van der Werf, 2006 2006 ) with three different statistical models (additive model, additive and dominance model, and additive, dominance and epistasis model). MATERIALS AND METHODS 1 Mixed linear model 2 A vector of observed phenotypes for the individuals in a pedigree is a linear function of fixed effects, polygenic effects, additive and dominance effects due to n QTL, epistatic interaction among the QTL, and residuals errors. The model can be written as (Cockerham, 1954 1954 ): $$y = X\beta + Zu + \sum\limits_{i = 1}^n {(Za_i + Zd_i )} + \sum\limits_{i = 1}^{n - 1} {\sum\limits_{j = i + 1}^n {(Za_i a_j + Za_i d_j + Zd_i a_j + Zd_i d_j )} } + e,$$ where y is a vector of phenotypic observations on the trait, β is a vector of fixed effects, u is a vector of random polygenic effects for each animal, a i and d i are vectors of additive and non-additive random effects due to the i th putative QTL, a i a j , a i d j , d i a j and d i d j are vectors of epistatic interactions between the i th and j th putative QTL, and e is a vector of random residual errors. The random effects in the model ( u , a i , d i , a i a j , a i d j , d i a j , d i d j and e ) are normally distributed with mean zero and variance A σ u 2 , G i σ a i 2 , D i σ d i 2 , G i G j σ a i a j 2 , G i D j σ a i d j 2 , D i G j σ d i a j 2 , D i D j σ d i d j 2 , and I σ e 2 , where σ 2 represents the variance of the component, A is a relationship matrix based on pedigree information, G i and D i are additive genetic and dominance relationship matrice at the i th putative QTL position, G i G j , G i D j , D i G j , and D i D j are the Hadamard product of the additive genetic and dominance relationship matrices at the i th and j th putative QTL positions, and I is an identity matrix. X and Z are incidence matrices for the fixed and random effects in the model. The associated variance covariance matrix ( V ) of all observations from the modeled is: $$V = ZAZ'\sigma _u^2 + \sum\limits_{i = 1}^n {(ZG_i Z'\sigma _{a_i }^2 + ZD_i Z'\sigma _{d_i }^2 )} + I\sigma _e^2 + {\text{ }}\sum\limits_{i = 1}^{n - 1} {\sum\limits_{j = i + 1}^n {(ZG_i G_j Z'\sigma _{a_i a_j }^2 + ZG_i D_j Z'\sigma _{a_i d_j }^2 } } + ZD_i G_j Z'\sigma _{d_i a_j }^2 + ZD_i D_j Z'\sigma _{d_i d_j }^2 ),$$ where Z ′ is the transpose of Z . The matrices G and D are constructed based on identical by descent (IBD) probabilities that are estimated using LDL information (Meuwissen and Goddard, 2000 2000 ; 2001 2001 ). For IBD estimation, a Markov chain Monte Carlo (MCMC) method is used, which is robust and efficient especially for complex pedigree, many markers and missing genotypes (Lee et al., 2005 2005 ). These G and D based on the IBD probabilities are incorporated into the QTL model selection in an empirical Bayesian approach (Casella, 2001 2001 ; Lee and van der Werf, 2006 2006 ). Reversible jump Markov chain Monte Carlo (MCMC) for multiple QTL 2 The number of QTL n , the position of each QTL ρ i ( i= 1,..., n ) and the model parameters $$\Theta = \{ \sigma _u^2 ,(\sigma _{a_i }^2 ,\,\sigma _{d_i }^2 ;{\text{ }}i = 1,...,n),(\,\sigma _{a_i a_{i + 1} }^2 ,...,\,\sigma _{a_i a_n }^2 ,\,\sigma _{a_i d_{i + 1} }^2 ,...,\,\sigma _{a_i d_n }^2 ,\sigma _{d_i a_{i + 1} }^2 ,...,\sigma _{d_i a_n }^2 ,\,\sigma _{d_i d_{i + 1} }^2 ,...,\sigma _{d_i d_n }^2 ;{\text{ }}i = 1,...,n - 1),\,{\text{ }}\sigma _e^2 \}$$ are unknown. Therefore, they should be estimated using Eq.( 3 3 ). The probability of estimated parameters given observed phenotypes is: $$pr(n,\rho ,\Theta |y) = \frac{{pr(y|n,\rho ,\Theta )pr(n,\rho ,\Theta )}}{{\sum {pr(y|n,\rho ,\Theta )pr(n,\rho ,\Theta )} }},$$ where pr ( y | n , ρ , Θ ) is the likelihood of the observed phenotypes given the estimated parameters, pr ( n , ρ , Θ ) is the joint prior probability of the estimated parameters and the denominator is summed over the probabilities of all possible parameter states. Since the computation of the denominator is not feasible due to a large parameter state space, an MCMC approach can be used for solving this problem (Lee and van der Werf, 2006 2006 ). In the process, the number of QTL and their positions are sampled from a proposal distribution. In a second step, residual maximum likelihood (REML) is used to estimate the model parameters for a given QTL model. The proposed model, variables and model parameters are accepted or rejected, according to the acceptance ratio from a reversible jump (RJ) MCMC from which the posterior QTL density is derived. Hence, a REML procedure is used nested within a Bayesian RJMCMC. When considering epistatic interaction among QTL, the number of random effects increases, e.g., the number of random effects due to n QTL is n for the additive model, 2 n for the additive and dominance model and 2 n +4 n ( n −1)/2 for the full model. However, not all epistatic components have significant effects in the model. It may not be necessary to include non-significant epistatic effects that do not improve the model likelihood. For a given set of QTL within a RJMCMC step, each epistatic component is tested and if the likelihood is not improved, the epistatic component will be removed from the model. Simulated data 2 One hundred generations of a historical population with effective size of 100 were simulated for 26 markers and 2 QTL in a 130 cM region. For the region from 10 to 20 cM and 110 to 120 cM as candidate regions, markers were densely positioned at 1 cM intervals. The potential QTL were simulated at 17.5 cM (QTL I) and 117.5 cM (QTL II). Unique alleles were assigned to QTL in generation 0, and one allele with moderate frequency (0.1~0.9) was randomly chosen to be the mutation in generation 100. In each generation, the number of male and female parents was 50 and their alleles were transmitted to descendents on the basis of Mendelian segregation using the gene-dropping method (MacCluer et al., 1986 1986 ; Meuwissen and Goddard, 2000 2000 ). The number of alleles assumed at each marker locus was 4 in generation 0 and starting allele frequencies were all at 0.25. The marker allele was mutated at a rate of 4×10 −4 per generation (Dallas, 1992 1992 ; Weber and Wong, 1993 1993 ; Ellegren, 1995 1995 ), i.e. a new allele was introduced as a mutation. Therefore, this historical population would have an equilibrium distribution of alleles in all marker loci and would generate LD among closely linked regions. Note that pedigree and genotype information was deemed not available for these 100 generations. In generation 100 and afterwards, phenotypic values for individuals were simulated as: y = μ + a QTL I + a QTL II + d QTL I + d QTL II + i QTL I×QTL II + e , where a is additive QTL effects, d is dominance QTL effects, i is the effects due to interaction between QTL I and QTL II. The population mean ( μ ) was 100, values for residuals ( e ) were from N (0, σ e 2 ) with σ e 2 =50. Two data sets were simulated. For the first data, the ratio of additive QTL variance over total phenotypic variance was 0 for the first QTL and 0.11 for the second QTL, the ratio of dominance QTL variance was 0 for both QTL, and the ratio of epistatic QTL variance was 0.2. For the second data, the ratio of additive QTL variance over total phenotypic variance was 0.15 for the first QTL and 0.06 for the second QTL, the ratio of dominance QTL variance was 0.08 for the first QTL and 0.1 for the second QTL, and the ratio of epistatic QTL variance was 0.09. For QTL mapping results, the posterior QTL density was estimated in RJMCMC LDL mapping with additive effects only (additive model), additive and dominance effects only (additive and dominance model), or additive, dominance and epistasis effects (full model). In all cases, marker genotypes and phenotypes were available for the last 2 generations (200 animals) used for analyses. RESULTS AND DISCUSSION 1 When there are complex interactions between QTL, the full model fitting additive, dominance and epistatic effects gives a higher mapping resolution than other reduced models. Fig. 1 1 shows that reduced models fitting additive effects only or additive and dominance effects only without epistatic effects could not catch the signal of the first QTL (Figs. 1a and 1b 1a and 1b ). This is probably due to the fact that the additive and dominance effects of the first QTL are negligible. However, the full model clearly maps the first and second QTL on the true position within a small region (Fig. 1c 1c ). Triangle shows the true QTL positions. The ratio of additive variance over phenotypic variance is ~0 and ~0.11 for the first and second QTL respectively. The ratio of dominance variance over phenotypic variance is close to 0 for all QTL. The ratio of epistatic variance over phenotypic variance is 0.2 Fig. 2 2 shows that the model fitting additive effects only could not detect the second QTL because the second QTL has small additive effects (Fig. 2a 2a ). However, when the model includes dominance, the second QTL can be detected (Fig. 2b 2b ). When the model includes epistasis in addition to dominance and additive effects, the mapping resolution becomes higher (Fig. 2c 2c ). Triangle shows the true QTL positions. The ratio of additive variance over phenotypic variance is ~0.15 and ~0.06 for the first and second QTL respectively. The ratio of dominance variance over phenotypic variance is ~0.08 and ~0.1 for the first and second QTL respectively. The ratio of epistatic variance over phenotypic variance is 0.09 From the results, when there are complex interactions between QTL, the model considering additive, dominance and epistatic effects can properly capture the signals of the multiple QTL. The empirical Bayesian approach with LDL information (used in this study) can help to simultaneously map multiple interacting QTL within a small region. Therefore, the present approach would be an efficient tool to study complex traits. References J Blangero , JT Williams , L Almasy . Quantitative trait locus mapping using human pedigrees . Hum Biol , 2000 . 72 ( 1 ): 35 - 62 . . Quantitative trait locus mapping using human pedigrees J Blangero , JT Williams , L Almasy . Variance component methods for detecting complex trait loci . Adv Genet , 2001 . 42 151 - 181 . . Variance component methods for detecting complex trait loci G Casella . Empirical Bayes Gibbs sampling . Biostatistics , 2001 . 2 ( 4 ): 485 - 500 . DOI: 10.1093/biostatistics/2.4.485 http://doi.org/10.1093/biostatistics/2.4.485 . Empirical Bayes Gibbs sampling CC Cockerham . An extension of the concept of partitioning hereditary variance for analysis of covariances among relatives when epistasis is present . Genetics , 1954 . 39 ( 6 ): 859 - 882 . . An extension of the concept of partitioning hereditary variance for analysis of covariances among relatives when epistasis is present JF Dallas . Estimation of microsatellite mutation rates in recombinant inbred strains of mouse . Mamm Genome , 1992 . 3 ( 8 ): 452 - 456 . DOI: 10.1007/BF00356155 http://doi.org/10.1007/BF00356155 . Estimation of microsatellite mutation rates in recombinant inbred strains of mouse H Ellegren . Mutation rates at porcine microsatellite loci . Mamm Genome , 1995 . 6 ( 5 ): 376 - 377 . DOI: 10.1007/BF00364807 http://doi.org/10.1007/BF00364807 . Mutation rates at porcine microsatellite loci SH Lee , JHJ van der Werf . Simultaneous fine mapping of multiple closely linked quantitative trait loci using combined linkage disequilibrium and linkage with a general pedigree . Genetics , 2006 . 173 ( 4 ): 2329 - 2337 . DOI: 10.1534/genetics.106.057653 http://doi.org/10.1534/genetics.106.057653 . Simultaneous fine mapping of multiple closely linked quantitative trait loci using combined linkage disequilibrium and linkage with a general pedigree SH Lee , JHJ van der Werf , B Tier . Combining the meiosis Gibbs sampler with the random walk approach for linkage and association studies with a general complex pedigree and multi marker loci . Genetics , 2005 . 171 ( 4 ): 2063 - 2072 . DOI: 10.1534/genetics.104.037028 http://doi.org/10.1534/genetics.104.037028 . Combining the meiosis Gibbs sampler with the random walk approach for linkage and association studies with a general complex pedigree and multi marker loci JW MacCluer , JL VanderBerg , B Raed , OA Ryder . Pedigree analysis by computer simulation . Zoo Biology , 1986 . 5 ( 2 ): 147 - 160 . DOI: 10.1002/zoo.1430050209 http://doi.org/10.1002/zoo.1430050209 . Pedigree analysis by computer simulation THE Meuwissen , ME Goddard . Fine mapping of quantitative trait loci using linkage disequilibria with closely linked marker loci . Genetics , 2000 . 155 ( 1 ): 421 - 430 . . Fine mapping of quantitative trait loci using linkage disequilibria with closely linked marker loci THE Meuwissen , ME Goddard . Prediction of identity by descent probabilities from marker-haplotypes . Genet Sel Evol , 2001 . 33 ( 6 ): 605 - 634 . DOI: 10.1051/gse:2001134 http://doi.org/10.1051/gse:2001134 . Prediction of identity by descent probabilities from marker-haplotypes BD Mitchell , S Ghosh , JL Schneider , G Birznieks , J Blangero . Power of variance component linkage analysis to detect epistasis . Genetic Epidemiology , 1997 . 14 ( 6 ): 1017 - 1022 . DOI: 10.1002/(SICI)1098-2272(1997)14:6<1017::AID-GEPI76>3.0.CO;2-L http://doi.org/10.1002/(SICI)1098-2272(1997)14:6<1017::AID-GEPI76>3.0.CO;2-L . Power of variance component linkage analysis to detect epistasis S Purcell , PC Sham . Epistasis in quantitative trait locus linkage analysis . Behavior Genetics , 2004 . 34 ( 2 ): 143 - 152 . DOI: 10.1023/B:BEGE.0000013728.96408.f9 http://doi.org/10.1023/B:BEGE.0000013728.96408.f9 . Epistasis in quantitative trait locus linkage analysis JL Weber , C Wong . Mutation of human short tandem repeats . Hum Mol Genet , 1993 . 2 ( 8 ): 1123 - 1128 . DOI: 10.1093/hmg/2.8.1123 http://doi.org/10.1093/hmg/2.8.1123 . Mutation of human short tandem repeats The authors thank the 3rd International Conference on Quantitative Genetics (ICQG) committee for giving an opportunity to present this study in the conference. The trip to ICQG was supported by the collaboration between University of New England and National Institute of Animal Science in Korea. 10.1631/jzus.2007.B0787 Project supported by the International Pig Improvement Company (PIC) and Sheep Genomics, Australia 2007-09-13 2007-09-16 2022-01-05 Copyright © 2007, Journal of Zhejiang University Science 2007 Journal of Zhejiang University-SCIENCE B (Biomedicine and Biotechnology) 11 08