Introduction

Paternal half-sib families are widely used as resource population for quantitative trait loci (QTL) mapping. Many methods have been proposed for such designs. Georges et al. (1995) developed a maximum likelihood method for a single-family design and implemented it to map QTL for milk production traits in the US Holstein population. A regression method of interval mapping proposed by Knott et al. (1996) has been commonly used to map QTL in half-sib families, particularly in dairy cattle (Spelman et al., 1996; Zhang et al., 1998; Heyen et al., 1999; Velmala et al., 1999; Nadesalingam et al., 2001; Plante et al., 2001; Ron et al., 2001; Freyer et al., 2002). Grignola et al. (1996a, 1996b) proposed a restricted maximum likelihood method for half-sib families, which was implemented by a number of investigators (for example, Zhang et al., 1998; Freyer et al., 2002; Liu et al., 2004). However, all of these studies are based on a single-QTL model, resulting in estimation bias due to the polygenic nature of quantitative traits. Zeng (1994) proposed a composite interval mapping approach, in which the effects of the markers that do not bracket the tested QTL are treated as a cofactor to absorb the effects of other QTL. This method was successfully used by De Koning et al. (2001) in a half-sib population. Although composite interval mapping can greatly improve the precision of QTL position estimation, it is essentially a single-QTL model. Kao et al. (1999) developed a multiple interval mapping approach, in which multiple QTL are simultaneously included in the model and the suitable multiple-QTL model is determined by sequential stepwise selection. Sequential stepwise selection is not optimal because of the dynamic changes of the null hypothesis, such that its efficiency is significantly influenced by the quality of the data (Raftery et al., 1997; Gelman et al., 2004).

With the development of the Bayesian model selection approach, the issues in multiple-QTL mapping have been solved both for inbred line crosses lines (Satagopan and Yandell, 1996; Heath, 1997; Uimari and Hoeschele, 1997; Sillanpää and Arjas, 1998; Stephens and Fisch, 1998; Gaffney, 2001; Xu, 2003; Yi et al., 2003, 2005) and for outbred populations (Meuwissen and Goddard, 2004; Liu et al., 2007; Yi and Xu, 2000). The stochastic search variable selection algorithm was first developed by George and McCulloch (1993) and used by Meuwissen and Goddard (2004) to fine map multiple QTL in outbred half-sib families. This method can simultaneously handle multiple QTL and also keeps model dimensions fixed. However, all of these authors assume that each marker interval contains one QTL and that the QTL position is fixed at the middle of marker, so that it cannot be applied for dense markers because of increased computational burden and correlations among potential QTLs. Liu et al. (2007) used the reversible jump Markov chain Monte Carlo (RJMCMC; Green, 1995) method for multitrait mapping of multiple QTL. The two methods described above are based on variance component models and use the efficient Gibbs sampler to first sample all QTL substitution effects and then update QTL variance using the sampled substitution effects. Because the parameter to estimate is large (including QTL substitution effect and variance), it converges slowly and has poor mixing.

The potential concern with RJMCMC is that it often has poor mixing due its variable dimensions (Yi, 2004; Wang et al., 2005; Liu et al., 2007; Banerjee et al., 2008) and thus is not easy to implement (Banerjee et al., 2008). To overcome this shortcoming of the RJMCMC algorithm, Godsill (2001) conceived a Bayesian composite model space approach. His method keeps the model dimensions fixed so that it mixes rapidly, and was used by Yi (2004) in QTL mapping for an allelic substitution model. In our previous work, we extended the Bayesian composite model space approach to a variance component model (Fang et al., 2009). This method also directly estimates variance components with the M–H algorithm, as Yi and Xu (2000) did, but the model dimension was unchanged. The efficiency of the method was illustrated with an outbred full-sib population and a series of simulated data. In dairy cattle, the half-sib families’ design is often used for QTL mapping. It is of interest to apply this method to multiple-QTL mapping in dairy cattle. In the present study, we modified this method for half-sib families for continuous traits as well as binary traits. The method optimally updates the covariance effect, polygenic variance and residual variance, so that the computational speed and efficiency are greatly improved compared with our previous method. Finally, we applied this method to analyse real data and a series of simulated data to fully assess its performance.

Bayesian modelling of multiple QTLs

Bayesian composite model space method

In a half-sib design, the parents of the half-sib families are randomly sampled from a large outbred population that is in Hardy–Weinberg and linkage equilibrium. Let y be an n × 1phenotypic vector, where n is the number of phenotypic observations. If y is continuously distributed, the multiple-QTL model can be expressed as

where Γ=(γ1,…γj,…γL)′ is the maximal QTL number, and γj is a binary variable indicating that the corresponding QTL effect is present (γj=1) or absent (γj=0) in the model; a=(a1, …, aj, …, aL) is an n × L matrix, and aj is an n × 1 random QTL effect vector and aj∼N(0,Θjσj2), where σj2 is the QTL variance and Θj is the identity-by-descent (IBD) matrix that can be inferred by the conditional expectation approach (Xu and Gessler, 1998); Z is an n × n identity design matrix relating records to each individual; β is a vector of covariate effects and X is the corresponding design matrix; g is denoted as an n × 1 vector of random polygenic effect and g∼N(0, A σA2), where A is the additive relationship matrix and σA2 is the polygenic additive variance; e is the vector of random error with the distribution e∼N(0,I σe2), where I is n × n identity matrix and σe2 is the residual variance. We assume that neither QTL nor polygenes have dominant effects. Because γj2=γj, the variance component model can be expressed as

Prior, likelihood and joint posterior

The prior of indicator variable γj follows a Bernoulli distribution and p(γj=1)=l0/L, where l0 is the expected QTL number and L=l0+l0√3 is the maximal QTL number (see the study by Yi et al., 2005 for details); the prior for QTL variance, polygenic variance and residual variance is assumed to follow informative inverted χ2 distributions with expressions:

and

where s (the scaled parameter) and ω (the degree of freedom) are hyperparameters; the priors of covariate effects are uniformly distributed; and the QTL position follows uniform distribution across the region. We denote marker information M, unknown variables θ={β,σ12∼σj2,γ1∼γj,σA2,σe2} and QTL position λ={λj}j=1L. Accordingly, the joint posterior probability of unknowns can be expressed as

where

and

Posterior calculation

Updating model indicator

We introduced an indicator {γj}j=1L in the variance components model, which is key to the new method. The fast M–H algorithm is used to update the indicator, which is faster than the Gibbs sampler because it can avoid calculating the denominator f(y,γj). Our algorithm first proposes a new value γj=κ, κ{0,1}, on the basis of its prior probability, and then accepts with a probability equal to min (1,r), where

Updating polygenic variance, residual variance and covariance effects

Joint updating strategy

To reduce the computational burden, we employed a joint updating strategy to generate the posterior distribution of σA2, σe2 and β via the efficient random walk Metropolis–Hastings (RWM–H; Metropolis et al., 1953; Hastings, 1970) algorithm. Specifically, new values σA2(*), σe2(*) and β(*) are proposed, which are accepted with a probability equal to min (1, r) with

where and hrβ, hrA and hre are proposal ratios for β, σA2 and σe2, respectively. The new proposal for the covariate effect β(*) is sampled from a multiple-dimension normal distribution with mean value β(0) and (co)variance Vt equal to the tuning parameter (Gelman et al., 1995; Browne, 1998). As the proposal is symmetric p(β(*)∣β(0),Vt)=p(β(0)∣β(*),Vt), the proposal ratio hrβ=1. The new proposal variance σ2(*), including σA2(*) and σe2(*), is proposed from a scaled inverted χ2 distribution with v degrees of freedom and a scaled parameter with the expectation of the current variance see Browne, 1998; Fang et al., 2009). Thus, the proposal ratio hr for variance component can be expressed as

Optimising the tuning parameter

In RWM–H analyses, the values of the degree of freedom v and the (co)variance Vt may influence the performance of the RWM–H algorithm. We introduce a general way to ascertain them in order to ensure the optimal performance of RWM–H.

When a multiple-dimension normal distribution with variance Vβ is simulated, the optimal tuning parameter of the RWM–H is suggested as Vt=(2.38/√d)2Vβ (Gelman et al., 1995), where d is the number of dimensions. However, in our joint sampling scheme, β is updated with the polygenic and the residual variance; hence, d equals the number of the covariates plus 2. In a Bayesian framework, the posterior distribution of β is normally distributed with variance Vβ=(XTV−1X)−1, where V is the unknown phenotypic (co)variance, which can be approximately estimated as V=σ′A2A+σe2I, where σ′A2 is genetic variance including polygenic and QTL variance. σ′A2 and σe2 can be determined by restricted maximum likelihood. If a binary trait is analysed, σe2=1, and thus only σ′A2 needs to be estimated. The degree of freedom is empirically set as v=60, and the whole setting of v and Vt may lead the acceptance rate falling into the range of 10∼40%. Such an acceptance rate range is optimal (Roberts and Rosenthal, 2001) because the target distribution can be efficiently explored and the autocorrelation is reduced dramatically under such conditions.

The process of updating QTL variance and QTL position is described in Appendix A In mapping QTL for a binary trait, we employ a data augmentation approach to generate the liability y, which has been illustrated at length by Yi and Xu (2000).

Applications

Real data analysis

The data from a Chinese Holstein population, with a daughter design, have been described by Chen et al. (2005). The population consisted of 26 bulls and their 2270 daughters, and the following five phenotypes of milk production traits were measured: milk yield (MY), fat yield (FY), protein yield (PY), fat percentage (FP) and protein percentage (PP). The bulls and daughters were genotyped at 14 marker loci on chromosome 6, covering a total distance of 55.7 cM. Marker names and their relative positions are shown in Figures 1 and 2. For each trait, the estimated breeding value (EBV) of each daughter was used as phenotypic data. The EBV was further standardised for convenience to choose the initial value (Jansen, 2003; Xu, 2003; Yi et al., 2005), and the standardised EBV was subject to Bayesian analyses.

Figure 1
figure 1

Profiles of the posterior probability (proposed method) and F-statistic (LS method) obtained from real data analysis. The black line and light grey line represent the BF score under the maximum QTL number L=3 and L=2, respectively, and the dark grey line represents the F-statistic from the LS method. The grey dashed horizontal line represents the threshold value of the LS method, and the black dashed horizontal line represents the standard threshold 3.2 of BF score. Figures a–e are the results of the five traits MY, FY, PY, FP and PP.

Full size image
Figure 2
figure 2

Profiles of the posterior probability from real data analysis using our method.

Full size image

Because maternal genotypes are not available in a half-sib design and treated as missing data, they are randomly sampled from possible alleles with equal probability. The Bayes factor (BF) defined by Yi et al. (2005) was used to support the existence of a QTL, and the BF threshold of 3 supports a claim of significance (Kass and Raftery, 1995). To compare the new method with the least square (LS) method, the real data were also analysed with the software QTL Express (Seaton et al., 2002).

In real data analyses, the hyperparameters ωj=ωA=ωe=3 lead to a rather non-informative prior (Fang et al., 2009); the expected QTL number l0=1 and thus the maximum QTL number L=l0+l0√3≈3. Accordingly, the prior probability of the model indicator p(γj=1)=l0/L=0.33; the prior distribution of QTL variance was truncated at 0.2 according to our prior information that the heritability explained by QTL is <20%; and the tuning parameter for QTL position was 2 cM. The MCMC ran for 21 000 rounds, the thinning interval was set to 10, and the burn-in period was 1000. As a result, the size of the posterior sample was 2000.

The profiles of the BF from the new method and the F-statistic of the LS method are plotted in Figure 1, and the profiles of posterior probability are presented in Figure 2. The profiles of F are the same as that in Chen et al. (2005) and only one QTL that affects FY was significant at the 0.01 level. However, with our method, four QTLs were detected, one affecting FY, one affecting FP and two affecting PP. Figure 1 shows that the peaks of the profiles of our method and the LS method almost overlap. The point estimates and highest posterior density regions of the QTL position from the two methods are listed in Table 1. The estimated heritabilities for the detected QTL are also given in Table 1; they are all about equal to 4%.

Table 1 The parameter estimates for the milk production traits using the proposed method
Full size table

The acceptance rate of the M–H algorithm for the joint updating of the covariance effect, the polygenic variance and the residual variance was 14.3% for MY, 15.5% for FY, 14.2% for PY, 18.9% for FP and 17.5% for PP and varied from 10 to 40%. The histograms of posterior distributions are plotted in Supplementary Figures S1–S4 of the Supplementary Material, and a typical output is shown in Figure 3, demonstrating that posterior distributions of the covariance effect, polygenic variance and residual variance are well explored. The results suggest that the new joint updating algorithm is more efficient.

Figure 3
figure 3

Profiles of the posterior probabilities of population mean (a), polygenic variance (b) and residual variance (c) from real data analysis for the trait MY.

Full size image

To test the sensitivity for the maximum QTL number (or expected QTL number), we also set the expected QTL number l0=0.73, leading to the maximum QTL number L=l0+l0√3=2 and thus the prior probability of the model indicator p(γj=1)=l0/L=0.37. The profiles of the BF are also plotted in Figure 1 and are similar to that from l0=1. The profiles of posterior probability are also similar to that from l0=1 (data not shown). The results show that the new method is not highly sensitive to L (or l0).

Simulation studies

We conducted two simulation experiments to further evaluate the performance of the method. In the first experiment, we simulated 20 independent half-sib families, each with 100 individuals, so that the final sample size was 2000. The parents of each family were randomly sampled from a large outbred population in linkage equilibrium and Hardy–Weinberg equilibrium. Eleven evenly spaced markers, each having six alleles with equal frequency, covered a 100 cM chromosome region with an average marker interval of 10 cM. Two QTL were simulated at positions 15 cM and 75 cM. The additive variances of the two QTL were σ12=0.3 and σ22=0.4, and the dominant variances were assumed to be absent. The QTL additive effects of parents were randomly sampled from N(0,σ12) and N(0,σ22). The alleles of the offspring at markers and QTL were inherited from their parents following Haldane's recombination rule. The polygenic variance was set to σA2=1.5. The residual error was sampled from a normal distribution with mean 0 and variance σe2=1.5. The overall population mean was set as 0. The phenotypic values of each offspring were the sum of the population mean, QTL effect, polygenic effect and residual error. Accordingly, the heritabilities explained by the two QTL were 8.1 and 10.8%. In this experiment, we only investigated the performance of the continuous trait, because the heritability of the simulated QTL is so low that neither the proposed method nor the LS method can give a precise estimate for a binary trait.

In the analysis, the expected QTL number l0=2, leading to the maximum QTL number L=5, and thus p(γj=1)=l0/L=0.4. The tuning parameter for QTL position was 2 cM. The thinning interval for MCMC was 10, the burn-in period was 1000 and the length of the complete chain was 21 000.

The BF from our method and the F-value from the LS method are plotted in Figure 4a and the posterior probabilities are in Figure 4b. With the new method, both simulated QTL were successfully detected. The parameter estimates are listed in Table 2, but only one QTL was found by the LS method at a 0.01 significance level. The results suggest that the new method is more powerful than the LS method. The parameter estimates are given in Table 2, and they are mostly close to their true values. The acceptance rate for the joint update of covariance, polygenic variance and residual variance was 23.3%, which also reflected the efficiency of our developed joint-updating strategy.

Figure 4
figure 4

Profiles of test statistics obtained from the simulated single-chromosome data. (a) Profiles of the posterior probability for L=5 and L=3 from the new method. (b) Profiles of BF score from the proposed method for L=5 and L=3, and profile of the F-statistic from the LS method. The true positions of the simulated QTL are represented by upward arrows. The grey dotted line and black dotted line present the threshold of the LS method at a 0.01 significance level and the BF score 3.2 of the new method, respectively.

Full size image
Table 2 The parameter estimates in the first simulated experiment
Full size table

To test the sensitivity for the maximum QTL number, we set the expected QTL number l0=1, leading to the maximum QTL number L=3 and thus p(γj=1)=l0/L=0.33. The results are plotted in Figures 4a and b and show that the profiles from L=5 and L=3 have very little difference.

To further compare the performance of the proposed method and the LS method, a receiver-operating characteristic was used to summarise the true positive rate (tp rate) and false positive rate (fp rate; Fawcett, 2006). Because of the high computational intensity and the difficulties in summarising the results of our method, we only replicated 30 experiments for the proposed method and the LS method. The threshold was divided on average into 20 points from the lowest threshold, which results in the maximum tp rate and minimum fp rate, to the highest threshold, which results in the minimum tp rate and maximum fp rate. In our study, the 20 successive thresholds ranged from 0.7 to 36 for the proposed method and from 1.7 to 7.3 for the LS method. At each threshold, the tp rate and fp rate were summarised for the two simulated QTLs together and are plotted in Figure 5. In our study, the maximum tp rate is not 1 because sometimes both methods failed to generate a clear signal of a QTL at the simulated position; thus, the QTL was treated as missing. As pointed out by Fawcett (2005), a point in receiver-operating characteristic space is better than another if it is to the northwest (tp rate is higher, fp rate is lower or both). Accordingly, our method always outperforms the LS method, as seen in Figure 5. Specifically, we summarised the results at the BF threshold of 3.0 for our method and the F threshold of 0.01 for the LS method. The statistical powers were 0.833 and 0.717 for our method and the LS method, respectively, and the fp rates were 0.033 and 0.067, showing that our method had a higher tp rate and fp rate than the LS method. The average estimates of QTL position over 30 replications were 16 (6.09) and 77.13 (6.95) for the first and second QTL using our method, and 22.25 (5.56) and 78.39 (7.21) with the LS method, where the values in parentheses are the s.d. estimates. The position estimation for the first QTL with our method is closer to the true value of 15 than that obtained using the LS method. An explanation of the results is that a multiple model is used in our method, so it tends to give less-biased estimates for QTL position. The estimates for the second QTL show no clear difference between the two methods. Moreover, the s.d. estimates with the two methods also show no clear differences.

Figure 5
figure 5

Receiver-operating characteristic (ROC) graph obtained from 30 replicated simulation experiments for the new method (grey line) and the LS method (black line).

Full size image

In the second simulation study, we focused on genome-wide mapping and simulated a large genome with length 2000 cM, with 201 evenly spaced markers covering it with a marker interval of 10 cM. Five QTLs were placed in the genome, and the corresponding positions and effects are listed in Table 3. The residual variance was set as 1.5 and the polygenic variance was absent. The total heritability was 69.4%, and the heritability explained by these QTL ranged from 8.1 to 22.4%. The family structure and other parameters were the same as in the first simulation. To simulate the phenotypic value for a binary trait, the continuous phenotypic values generated above were treated as liabilities and the threshold was set as 0, resulting in 50% incidence. Specifically, if the liability is below 0, the binary phenotype value is set to 0; otherwise it is set to 1.

Table 3 The parameter estimates in the second simulated experiment
Full size table

In the analysis, we set the expected QTL number l0=3, leading to the maximum QTL number L=8 and thus the prior probability of the model indicator p(γj=1)=l0/L=0.38. The tuning parameter for QTL position was set as 20 cM. The thinning interval of MCMC was 10, the burn-in period was 1000 and the length of the complete chain was 21 000. When binary traits were analysed, the residual variance was re-parameterised (Albert and Chib, 1993) and equal to 1.5.

For continuous trait analysis, the posterior probabilities of QTL position from the new method are plotted in Figure 6, and the BF from the new method and F-statistic from the LS method are plotted in Figure 6b. The results for binary-trait analysis are shown in Figure 7. Using our method, all simulated QTLs were identified both for continuous traits and for binary traits. However, by the LS method, one QTL was missed for continuous trait analysis and two were missed for binary trait analysis. The simultaneous modelling of multiple QTL in our method likely contributes to a gain in statistical power. The parameter estimates are given in Table 3. Generally, the position estimates are largely close to their true values, both for continuous trait analysis and for binary trait analysis, but the estimation for the second QTL in binary trait analysis deviates from its true value by ∼20 cM; the relative lower heritability of this QTL (12.2%) and the information loss in binary trait analysis could explain the biased position estimate. Nonetheless, the true positions of all QTLs are all within their highest posterior density region. The estimates of QTL variance are also given in Table 3, and they also deviate from their true values to some extent. In fact, for mapping QTLs in half-sib design, only paternal genotypes are available and the maternal genotypes are completely unknown; thus, the variance components are difficult to estimate accurately. Finally, we summarised the acceptance rate of M–H run for the joint update of covariance, polygenic variance and residual variance, which equalled 36.8 and 48.8% for continuous and binary traits, respectively, and were near their optimal regions.

Figure 6
figure 6

Profiles of statistics from the simulated genome-wide data for continuous trait. (a) Profiles of the posterior probability by the proposed method. (b) Profiles of the BF score by the proposed method and the F-statistic by the LS method. The grey horizontal line represents the thresholds for the LS method at a 0.01 significance level. The standard threshold 3.2 of BF score is used for declaring the QTL by the new method, which is not shown in the figure for clear presentation. The true positions of the simulated QTL are represented by upward arrows.

Full size image
Figure 7
figure 7

Profiles of statistics obtained from the simulated genome-wide data for a binary trait. (a) Profiles of the posterior probability by the proposed method; (b) profiles of the BF score by the proposed method and the F statistic by the LS method. The grey horizontal line represents the thresholds for the LS method at a 0.01 significance level. The standard thresholds 3.2 of BF score is used for declaring the QTL by the new method, which is not shown in the figure for clear presentation. The true positions of the simulated QTLs are represented by upward arrows.

Full size image

We also varied the expected QTL number l0 =4 and 5, leading to the maximum QTL number L=11 and 14 and p(γj=1)=l0/L=0.363 and 0.357. The profiles of posterior probability and BF are similar to that from L=8 (data was not shown). These results further illustrate that our method was not very sensitive to the expected QTL number or the maximum QTL number.

Discussion

We employed the approach reported by Yi et al. (2005) to ascertain the maximum-QTL number and found that the results were not very sensitive to the maximum QTL number. Theoretically, the maximum-QTL number may be set as any value as long as it is greater than the actual QTL number. The simplest method assumes that each marker interval contains one QTL, but the computational burden will be increased because of the large QTL number, and at the same time, the efficiency will be challenged. Xu (2007) developed a method to handle such a large model using a hierarchical prior for γj, and p(ρ)=Bata(1,1), which is useful when the number of model effect is large, especially for the investigation of QTL epistatic effects and QTL-environment interaction effects.

In our real data analysis, we found that marker BMS143 affects FP traits. This marker has also been reported by others to affect FP (Schnabel et al., 2005), MY (Gomez-Raya et al. 1998; Velmala et al., 1999) and PP (Spelman et al., 1996; Zhang et al., 1998; Velmala et al., 1999; Schnabel et al., 2005), but we only confirmed that it affects FP. Moreover, we found two extra QTL that affected PP traits. For the first extra QTL, the position estimates were close to the marker ILSTS97, which was also reported by Moisio et al. (2000) to affect the PP trait. For the second, the estimated position was close to CSN3 and CSN1, as reported by many researchers (Velmala et al., 1995; Maki-Tanila et al., 1998; Leone et al., 1998; Ikonen et al., 1999), and also close to AFR227, as reported by Velmala et al. (1999), and to BP7, as reported by Ashwell et al. (1998) and Schnabel et al. (2005), all of which were found to affect the PP trait.

The variance component model can be written as and R is a diagonal matrix with the diagonal element being the inverse of reliability for each individual 1/ri. If the difference of ri between each individual is small, the reliability may have negligible effects on the results; however, if the difference is large, it is necessary to take such information into account. Although we did not consider the reliability in real data analysis, the computational package provided by us (described below) can include such information.

In our real data analysis, we assumed that the markers of parents are in linkage equilibrium, but this assumption may not be reasonable if the marker interval is relatively small. A method that combines the information of linkage disequilibrium (LD) and linkage analyses called LDL can significantly improve the estimate precision for QTL position if the markers of parents are in linkage disequilibrium (Meuwissen and Goddard, 2000, 2001; Meuwissen and Goddard, 2004; Lee and Van der Werf, 2006). LDL also estimates variance components, and therefore the proposed multi-QTL mapping method is completely applicable if the appropriate IBD matrix is constructed.

We only consider QTL main effects in this study and ignore the interacting effects between QTL for simplicity. In our variance component model, the interactions can be included by introducing Hadamard multiplication of matrices. However, when the loci are linked, the IBD covariance components will be difficult to construct and the performance of our method will be affected if such information is ignored. Moreover, we can also extend our method to categorical trait analysis as described in Appendix B.

We investigated the efficiency of the new method in half-sib families. In practice, if the appropriate IBD matrix is constructed, the new method also can be used for a more complex population, such as the mixture of full-sib and half-sib families or more complicated mating designs. However, the method is not useful for large populations because it is infeasible to calculate the inverse of a large (co)variance matrix. One solution is to incorporate an empirical Bayesian procedure to estimate variance components (Lee and Van der Werf, 2006), such as using REML (restricted maximum likelihood) estimates, which need not calculate the inverse of (co)variance matrix.

In the new method, calculating the IBD matrix in MCMC may cost a great deal of computational time. To solve the problem, the half-stored technology was employed. We first calculated and saved the IBD matrix on a hard drive for all positions and then read the IBD for particular positions in the current iteration. As a result, computational time was dramatically reduced. All computational analyses in our research were implemented via our computational package ‘BayesMapQTL.exe’, which was built in FORTRAN language and is available upon request.