Comparison of three multitrait methods for QTL detection

A comparison of power and accuracy of estimation of position and QTL effects of three multitrait methods and one single trait method for QTL detection was carried out on simulated data, taking into account the mixture of full and half-sib families. One multitrait method was based on a multivariate function as the penetrance function (MV). The two other multitrait methods were based on univariate analysis of linear combination(s) (LC) of the traits. One was obtained by a principal component analysis (PCA) performed on the phenotypic data. The second was based on a discriminate analysis (DA). It calculates a LC of the traits at each position, maximising the ratio between the genetic and the residual variabilities due to the putative QTL. Due to its number of parameters, MV was less powerful and accurate than the other methods. In general, DA better detected QTL, but it had lower accuracy for the QTL effect estimation when the detection power was low, due to higher bias than the other methods. In this case, PCA was better. Otherwise, PCA was slightly less powerful and accurate than DA. Compared to the single trait method, power can be improved by 30% to 100% with multitrait methods.


INTRODUCTION
The use of genetic markers to detect quantitative trait loci (QTL) has been well reported during the past 15 years. Statistical methods have been developed to detect QTL on the basis of genetic maps, mostly using maximum likelihood [23] or some linearised approximations [10,16,17]. In many studies, numerous traits are recorded but the analyses are carried out trait by trait. They often suggest an influence of the same chromosomic region on several traits (e.g. [2][3][4]). But whereas methods to detect multiple QTL on different linkage groups have been explored [13,29,34], tests for pleiotropy are still relatively limited. Different methods have been proposed, often derived from single trait methods. Some authors [14,19,20] have proposed to apply multivariate analysis for multitrait QTL detection to inbred lines, and Ronin et al. [29] have extended this to half sib families. Nevertheless, their applications to a mixture of full and half sib families of livestock populations are greatly limited. Indeed, in such populations, the QTL parameters must be estimated for each potentially heterozygote individual of the population, leading to a huge increase of the number of parameters estimated.
To overcome the increase of parameters, Weller et al. [31] have proposed to synthesise most of the information in a linear combination of the traits. This was obtained performing a principal component analysis (PCA) from the phenotypic covariance matrix. Such a transformation leads to as many phenotypically independent variables as traits, analysed by a univariate method. Since the same transformation is applied whatever the tested location, it does not integrate the putative effect of the QTL. To take this into account, we developed an algorithm based on the discriminate analysis (DA) to obtain one linear combination of the traits, specific to each tested position and analysed by a univariate method.
The purpose of this study was the comparison of three multitrait QTL detection methods developped for the application to experimental designs obtained in animal breeding populations, such as pig crosses [4]. We programmed a multivariate method, a principal component analysis and a discriminate analysis as options in a QTL detection software. The power, the accuracy of parameter estimates and the computing time of the three multitrait methods were compared to those of a single trait method. The comparisons were based on simulated data for a mixture of full and half sib families.

Hypothesis
For the single trait tests, the null hypothesis H0 "there is no QTL for the trait on the linkage group" was tested versus the alternative hypothesis H1 "there is a QTL influencing the trait at the x position". For the multitrait tests, the null hypothesis H00 "there is no QTL for the traits on the linkage group" was tested versus the alternative hypothesis H01 "there is a QTL influencing at least one trait at the x position".

Likelihood ratio tests
The likelihoods under the null hypothesis are independent of the position. For single trait tests, it was noted 0 Λ l . Under H1, the likelihood is calculated for each tested position: Λ x l . The likelihood ratio test (LRT) is then LRT lx = −2 ln( 0 Λ l /Λ x l ). For multitrait tests, the likelihood under H00 is noted 0 Λ, and under H01 Λ x . The likelihood ratio test is then LRT x = −2 ln( 0 Λ/Λ x ).

The single trait method (ST)
The maximum likelihood univariate method (ST) used for single trait detection of QTL has been described in Le Roy et al. [24]. It is an interval mapping method [23] considering a mixture of half-sib and full-sib families. The population is considered as a set of n sire families (i = 1, . . . , n), with n i mates for the sire i ( j = 1, . . . , n i ) and n ij progeny for the dam ij (k = 1, . . . , n ij ). Some hypotheses have been made to improve the robustness and the time of computation, based on [5,8,27]. In particular, only the most probable sire genotype was retained in the calculation, and the likelihood was linearised within full-sib families. It could not be totally linearised because of the uncertainty of the dam genotype. Following notations by Elsen et al. [5], for each trait l (l = 1, . . . , p), under H1, the likelihood function is then written at the x location as: where M i is the marker information of the sire family i; hs i is the most probable sire i genotype for genetic markers; hd ij is a dam ij genotype having a probability greater than 0.1 for genetic markers; yp ijkl is the phenotype for the quantitative trait l of the progeny ijk.
The penetrance function f is assumed to be normal and written as: σ 2 il being the variance of the sire i family for the quantitative trait l, and where p d x ijk = (q s , q d )/ hs i , hd ij , M i is the transmission probability of the couple of haplotypes q s and q d , i.e. the probability that the progeny ijk received from its sire the chromosome segment q s at the tested location x (q s = 1 from the grand sire, q s = 2 from the grand dam) and received from its dam the chromosome segment q d at the tested location x (q d = 1 from the grandsire, q d = 2 from the grand dam); µ xq s il is the mean of performances of the offspring having received the q s segment at location x from the sire i for the quantitative trait l; µ xq d ijl is the mean of performances of the offspring having received the q d segment at location x from the dam ij for the quantitative trait l.
µ il , µ ijl , α x il and α x ijl being respectively within half-sib and within full-sib, the phenotypic means and the average effects of the QTL substitution for the quantitative trait l, µ xq s il and µ xq d ijl can be written as: Thus, for each trait l, three parameters (family variance, mean and QTL substitution effect) per male, and two (family mean and QTL substitution effect) per female are estimated. Residual variances are supposed to be independent of the QTL genotype.
Under H0, the likelihood function is independent of the genotypes: where µ il and µ ijl are the phenotypic means as defined in the previous paragraph. We thus compared: H0: α il = α ijl = 0, ∀i, j, l versus H1: ∃i or ij that leads to α il = 0 or α ijl = 0. Then, the number of degrees of freedom of the test depends on the number of the effects estimated: n + n i=1 n i + 1.

The multivariate analysis (MV)
The multivariate method integrates all the traits using a multivariate function. Before performing multitrait analyses, the traits must be standardised to avoid confusion between the scales of measure and sizes of effects. Then, the likelihood expression can be written under H01 as: (5) where yp ijk = {yp ijkl ; l = 1, . . . , p} is the vector of the standardised performances for the progeny ijk.
The penetrance function is a p-dimension multinormal function. According to notation (3), it can be written as: where Y ijk = {Y ijkl ; l = 1, . . . , p} and Y ijk is the transposed vector; VC i is the residual covariance matrix between traits for the sire family i, with VC −1 i as the inverse and |VC −1 i | the determinant of this inverse: where ρ ll (l = l , l = 1, . . . , p, l = 1, . . . , p) is the residual correlation between the traits l and l .
The analytical expression of the likelihood under H00 is easily derived from equation (4) applying a multivariate penetrance function. We thus tested: H00: α il = α ijl = 0, ∀i, j, l versus H01: ∃ i or ij such as α il = 0 or α ijl = 0. Then, the number of degrees of freedom is (n + n i=1 n i )p + 1. Except ρ ll , estimated parameters were the same as for the ST analyses, but they were estimated in one likelihood maximisation. In the following, residual correlations were assumed independent from the sire families, to limit the number of parameters. Then, for the p traits, (p − 1) p/2 residual correlation coefficients were estimated. As for the ST analyses, residual variances are estimated within the sire families and supposed to be independent of the QTL genotype. This last assumption has been validated by many authors as being more powerful when it is true, but a source of great losses of power when it is wrong [20,27,29].
Jiang and Zeng [14] and independently Korol et al. [20] used such multivariate strategies but applied to inbred lines. It was mostly validated by plant data analysis [9]. It has been shown to be more efficient compared with single trait analyses if, among the p traits considered, for at least one couple of traits (l, l ), the product P ll of the QTL effects and the residual correlation α l α l ρ ll is negative.
Compared to inbred line strategies, the number of parameters to estimate is higher, potentially limiting power and increasing computing time. Wu et al. [32], and similarly Knott and Haley [19], proposed to partly overcome the computing time problem performing a least square analysis instead of a likelihood ratio test. But it does not reduce the number of parameters to estimate. Moreover, the differences of the results between the approaches have been clearly demonstrated, particularly when analysing complex genetic architecture [15].
The two following methods were developed to deal with the increase of the number of parameters. They are based on a univariate analysis of a synthetic linear combination of the traits in which most of the information from the data is summarised. Thus, the number of degrees of freedom of each test is the same as those of single trait tests.

The principal component analysis (PCA)
Principal component analysis (PCA) applied to multitrait detection of QTL was first proposed by Weller et al. [31], applied to a marker by marker analysis of dairy cow data. As [31], we performed PCA from the phenotypic covariance matrix of the data, considered as an estimation of the residual covariance matrix. p phenotypically uncorrelated linear combinations of the traits are derived from the components of the eigenvectors of the phenotypic covariance matrix. Each eigenvalue represents the part of phenotypic variability explained by the associated principal component variable.
The principal component variable associated with the eigenvalue λ m (m = 1, . . . , p for increasing λ m ) will be denoted PCAm. A univariate analysis of each PCAm is performed as for ST analysis, with a normal distribution as the penetrance function. Familial means, residual variances and QTL substitution effects are estimated for each synthetic variable. According to equation (3), the penetrance function for PCAm is: where yq ijkm = x m yp ijk is the phenotype of the progeny ijk for PCAm; x m = {x ml ; l = 1, . . . , p} is the eigenvector associated with λ m and σ 2 im is the variance of the sire family i for PCAm.
PCA has previously been analytically analysed [26]. It is asymptotically equivalent to the multivariate method if (1) the transformation is performed from the residual covariance matrix (e.g. estimated multivariately) and (2) the test statistic at each position is based on the sum of the test statistics of each principal component variable (sPCA). Consequently, PCA performed from the residual variability has the same limits as MV concerning its ability to improve the results of ST analysis. But this implies carrying out the transformation from the residual covariance matrix of the traits, which is difficult to access. However, when the residual covariance matrix does not differ much from the phenotypic one, PCA performed from the phenotypic matrix can produce similar results than when it is performed from the residual one [26].
As an example, in the simplest situation where phenotypic correlations ρ between all couples of traits are equal, λ 1 = 1 + ( p − 1)|ρ| and other are all λ m = 1 − |ρ|, m = 2, . . . , p. The eigenvector associated with λ 1 , x 1 , is then {x 1l ; l = 1, . . . , p} with |x 1l | = | √ p/p|, the sign depending on the sign of the correlation between the traits. For other eigenvectors, x 1 x m = 0, m = 2, . . . , p. Whatever the inner structure of the traits, the linear combinations then only depend on the phenotypic variability. As previously mentioned, PCA was performed from the phenotypic matrix. We also performed sPCA as a potential alternative to MV. The amount of phenotypic variability explained by each PCAm, λ m , cannot be used to limit the number of variables to analyse. Indeed, even low λ m can represent a significative part of the residual variability, depending on the relationship between the traits. Thus, all the variables must be analysed. But the interpretation of the results is difficult, since individual variables, phenotypically defined, do not only reflect the putative QTL variability at the tested position. Consequently two (or more) variables can suggest QTL at highly linked positions [19], meaning that either two different QTL are linked or a pleiotropic QTL determined the two linear combinations. Finally, testing for correlated variables must be taken into account in the calculation of the thresholds.
The following method was aimed at dealing with these two limits: the number of variables and the transformation from the residual covariance matrix.

The discriminate analysis (DA)
The discriminate analysis (DA) is a statistical method applied to individuals split into groups due to the disparity of the effect of several variables. It is aimed at finding the linear combination z of the variables which best discriminates these groups [28]. Analytical studies have shown that z maximises the ratio z Bz/z Wz of the between group variability B and the within group variability W. z is the eigenvector associated with κ, the greatest eigenvalue of W −1 B. Applied to multitrait QTL detection, this linear combination maximises the ratio of the genetic and residual variabilities specific to the putative QTL when performed at the QTL position.
Assuming a QTL is segregating, the residual and genetic covariance matrix due to the QTL must be estimated at its location. If the QTL allele inherited from a parent (say the sire) by the progeny is known, we could distinguish two groups of progeny depending on these alleles. The residual covariance matrix R, due to any reason but the QTL, would then be estimated from the within group variability. The genetic covariance matrix G, due to the QTL, would correspond to the between group variability. z would then characterise the sire QTL allele effects. Consider a case with two traits of equal residual variances σ, a residual correlation noted ρ 12 , genetic variances respectively noted σ g1 and σ g2 determined by a pleiotropic QTL, so that the QTL genetic correlation between those traits is 1. The components z 1 and z 2 of z associated with κ would be related by the expression: Korol et al. [22] proposed a similar method but applied to inbred lines. For each interval between two genetic markers, groups of haplotypes were constructed from non recombinant progeny (ignoring double recombination), considering the origin of the haplotype they inherited for this chromosomic segment to be known. Then, they applied a transformation which maximises the ratio of the genetic variability due to the chromosomic segment (equivalent to B in the formal description) and the total variability (T). The linear combination which maximises W −1 B is known to be identical to the linear combination which maximises T −1 B [28]. The principle of this strategy is then similar to applying DA to groups of haplotypes determined from genetic intervals. However, in animal populations, informative genetic intervals are generally more than 20 cM-sized, and different from one family to another. This implies that more than 20% of the population could be eliminated from the analysis due to recombination events, whereas the number of progeny is a limiting factor in these analyses.
To estimate the genetic and residual covariance matrix, we proposed to distinguish the two groups of haplotypes inherited from the sire at each tested position: each progeny belongs to both groups, its performance for each trait weighted by the probability that it received the corresponding haplotype (deduced from the p d The estimators of the genetic and residual variabilities are then the between and within group sums of squares and products of the weighted performances. The linear combination of the traits which best discriminates the groups at the position x, z x , is then calculated as previously. z x is univariately analysed, with a normal penetrance function. The likelihood is written as for PCA, but the coefficients of z x depend on the position tested. This implies maximising the likelihood under the null hypothesis at each tested position, since the variables are position specific. Considering a mixture of half-sib and full-sib families, groups can be defined within sire (two groups per family) or within dam and sire (three or four groups per family). In intercross designs with fixed QTL alleles in grand parental breeds, all F1 carry the same heterozygous QTL genotype. Groups can then be defined at the population level. If the QTL determinism is strictly additive, it can be demonstrated that the linear combination is the same when considering two, three or four groups (Appendix). In the present study, the experimental designs followed such hypotheses (see below section 2.4). Thus, DA was implemented dividing the population into two groups according to the breed origin of the haplotype received from the sire. Moreover, since the sire and dam QTL alleles are identical, defining the linear combination of the traits which best discriminates the QTL alleles of the sire means defining the linear combination of the traits which best discriminates the QTL alleles of the dam.

Axis of comparison
ST, MV, PCA and DA were compared for their ability to detect QTL using simulated data sets. They were each successively applied to each replicate.

Threshold calculation
Thresholds of rejection of the null hypothesis were estimated through simulations. They were performed assuming that there was no QTL on the linkage group. Couples of traits were simulated with the following residual correlations: ρ 12 = −0.4, ρ ll = 0.4 for l = l , l > 2, l > 3. The values of the residual correlation were previously confirmed not to influence the threshold (data not shown). On every replicate, the maximum of the test statistic of each method obtained on the linkage group was retained. The thresholds at the 5% level were then estimated as the 5% quantiles of the empirical distributions of the test statistics obtained with 2000 simulations [12]. ST and PCA implied the analysis of several variables. The type-I errors then required correction to obtain a general 5% type-I error, due to the multiple tests performed. When analysing independent variables, a Bonferroni correction of the type I error can be performed. Its approximation corresponds to the division of the type-I error by the number p of variables analysed. But carrying out such a transformation on correlated variables leads to too stringent thresholds. The independence between the traits analysed must then be checked before correction of the type-I errors. In the cases we simulated (see below), the QTL had similar effects on each of the traits, so that each trait generally contributes similarly to the linear combination. Thus, the approximated Bonferroni correction can be applied directly considering that the p traits are independent. As references, we will present the corrected and uncorrected thresholds and the corresponding powers for both of these methods in some cases.

Power
One thousand simulations were performed with a QTL to estimate the power and the accuracy of the estimates. The maximum value of each test statistic on the linkage group was retained, with corresponding position and QTL effect estimates. Powers were estimated as the percentages of maximum of the test statistics higher than the corresponding empirical thresholds.

Accuracy of estimates
The accuracy of the estimates was assessed by considering their mean squared errors (MSE) obtained over the 1000 replicates performed to estimate power. This characteristic synthesises bias and standard deviation. In some particular situations for position estimates, both will be presented separately. To make the MSE of QTL substitution effects of DA and PCA comparable with those of ST and MV, they were corrected by the square of the average over the 1000 simulations of the sums of the p coefficients of the linear combination at the maximum of the statistic test.

Computing time
Although the programme could certainly be optimised, we compared the CPU time required for each method to analyse a data set. We averaged computing times over the 1000 replicates performed under the alternative hypothesis for each type of experimental design.

Designs of the simulations
The animals came from an intercross design. There were 10 males and 10 females in the grand parental generation (F0), 10 unrelated males in the F1 parental generation and 500 progeny in the F2 generation. Each F1 male had 50 progeny with one or two unrelated female(s).
The linkage group was 100 cM-long, with equally spaced genetic markers. Each genetic marker had five isofrequent alleles identical in both the grand parental population. The grand parental haplotypes were randomly drawn. Their transmissions were simulated following Mendelian rules.
A pleiotropic QTL was simulated at position 31 cM, which almost corresponds to the middle of an interval whatever the marker density. Its alleles were assumed to be fixed in the F0 generation. Thus, all the F1 individuals were heterozygous for the QTL. The allele transmissions were simulated jointly with the marker alleles.
The progeny phenotypes followed a completely additive model. For each trait, besides the QTL, polygenic heritability was fixed to 0.2 and residual variance to 1. The traits were multinormally distributed, with residual correlations varying from −0.9 to 0.6.

Situations simulated
We grouped the simulations into three types of comparisons: (1) Depending on the number of traits: each female had 25 progeny, the marker density was 12.5 cM, and the QTL determined from two to six traits. The substitution effect of the QTL (α) was 0.50 on each of the traits. Residual correlations were 0.4 between all couples of traits in one set of simulations. In the other set, the residual correlation between traits 1 and 2 was −0.4, with all the others being 0.4.
(2) Depending on the genetic determinism of the traits: each female had 25 progeny, the marker density was 12.5 cM. Only two traits were considered. In a first set, the QTL affected the two traits with equal substitution effects (α = 0.50 or 0.30) and a residual correlation between −0.9 and 0.6. In a second set, the QTL determined only the trait 2. The aim was to quantify the ability of multitrait methods to exploit the residual correlation to detect a single trait QTL. Two levels of substitution effects (α = 0.5 and 0.3) and two levels of residual correlations (0.6 and −0.6) were simulated.
(3) Depending on the informativity of the design: the QTL determined two traits, with equal substitution effects (α = 0.50 and 0.30) and a residual correlation of −0.4. Two levels of marker density (50 or 12.5 cM) and two types of designs (25 or 50 progeny per female) were simulated.

RESULTS AND DISCUSSION
The results will be segmented depending on the influence of the characteristics previously described in Section 2.5. Table I presents the averaged computing times. For ST and PCA, they are similar. Since adding one trait means analysing one more variable, they increased linearly with the number of traits. For DA, the computing time, as the number of variables analysed, was independent of p. For MV, it was 5 to 100 times higher than that of ST, due to the number of parameters estimated.

Thresholds
Table I also presents empirical thresholds for each method depending on the number of traits p. For ST and PCA, 5% thresholds for the analysis of one variable are given (STv and PCAv respectively), jointly with global thresholds obtained after an approximate Bonferroni correction. The thresholds of STv and PCAv were independent of the number of traits, instead of thresholds of DA which slightly increased with p. When correcting the thresholds of ST and PCA, they turned out to be similar to those of DA. This increase must be due to the increase of the design informativity as p increases.
Concerning MV, the thresholds clearly increased with the number of traits, due to the increase of the number of parameters estimated. The thresholds for sPCA were equivalent to those obtained for MV, which was expected [26]. This offers a much faster strategy to obtain thresholds for MV. Table II presents the power obtained when the product of the effects and residual correlation for the first couple of traits (P 12 ) is either negative or positive. In the two cases, the power for ST presents the same tendency, decreasing with the number of traits, as the corrected thresholds increase with p.

Power
For the multitrait methods, the power is much higher when P 12 is negative. Similar trends have already been reported and analytically analysed [14,20,22,26] for inbred populations. Best power is obtained with multitrait methods when, at least for one couple of traits (l, l ), P ll is negative. This means that pleiotropic QTL are best detected when the genetic correlation (due to the QTL) is opposite to the residual correlation.
The power of MV was reduced by the increase of the number of parameters due to the increase of p. The power of sPCA was similar but generally slightly lower than that of MV, due to the transformation. The power of DA was the greatest, and tended to increase with p, whereas the power of other methods mostly decreased. When P 12 was positive, the power of PCA was slightly lower than that of DA (by 15% on average). When P 12 was negative, it greatly decreased with p. This was directly due to the construction of the principal component transformation. Since PCA is carried out from the phenotypic data, the PCA1 in those examples explained the greatest residual variability only when all the residual correlations were in the same direction. When some differed from the general direction, the analysis of PCA1 missed those directions, implying a loss of power. This was illustrated by the power of PCA for two traits when P 12 was negative, obtained from PCA2. In this case, it explains the lowest phenotypic variability but the greatest residual variability. Table II presents both bias and standard deviations of the QTL position estimates. Their means are respectively 0.073 and 0.025 for ST. The ranking of the multitrait methods was similar to that obtained for the power, except for PCA when P 12 was positive: it was more accurate than any other methods, representing a 25 to 50% improvement compared to DA bias, and decreased with the number of traits analysed. With ST, the position estimates are usually biased when the power is low when analysing a mixture of full and half sib families [18,27]. In such situations, multitrait methods hugely improve the position estimate accuracy.

Accuracy of QTL effect estimates
The MSE of the QTL effect estimates are summarised in Table II. They are 0.03 on average for ST, where the square of the bias represents less than 1% on average. Ranking of the multitrait methods was different from the previous cases. For MV, the MSE of the QTL effect estimates are lower than those of ST, with a difference of 0.0035 on average. The square of its bias represents 1.6% of the MSE. For PCA, they are half of those of ST. The trends depending on p and the sign of P 12 were similar to those of the MSE of the position estimates. For this method the square of the bias represents 0.16% of the MSE on average. For DA, the MSE are similar or greater than for ST when P 12 is positive and lower than the other methods when it is negative. This was due to an increase of the bias of the QTL effect estimate with the number of traits in the first case, from 0.018 for two traits to 0.075 for six traits (0.017 on average for ST). Indeed, the coefficients of the linear combination of DA tend to be null when the power is low, leading to a very small mean of the sums of the coefficients used to correct the MSE of DA to make them comparable to those of ST and MV (not shown).
For a mixture of full and half sib families, multitrait methods can then greatly improve the QTL detection compared to single trait strategies, as shown for inbred lines [14,20,22,26]. The comparisons demonstrated the good ability to detect QTL for the methods based on linear combination of the traits and the related gain of computing time when analysing several traits. Since the number of correlated traits in such studies can often be greater than 10 (i.e. [3]), this is strategically crucial. Table III summarises the power of the methods depending on the residual correlation and the equal substitution effects α for two traits (point 2 described in section 2.5). When α = 0.50, the power of ST is 21.9% on average. It decreases to 6.5% when α = 0.30. The decrease of power for multitrait methods previously commented in relation with the respective directions of the genetic and residual correlations was observed. When α = 0.50, these powers decrease from almost 100% when the residual correlation is −0.9 to 25% when it is 0.6, similar to the power of ST.

Power
We thus tested the ability of multitrait methods to improve the power of ST when ρ 12 = 0.4 in more favourable conditions, with twice the number of progeny in the design. On the one hand, we doubled the number of sires and dams, and on the other hand, we doubled the number of progeny per dam. In the first case, the power of ST was 37.0%, and it was improved by multitrait methods to 53.5%, 74.3% and 68.4% respectively for MV, DA and PCA.    But since this strategy is practically more realistic than the second one, its consequences on computing times were dramatic, with 18.5, 12, 76.12 and 19 seconds per implementation respectively for ST, DA, MV and PCA. In the second case, the power of ST was 60.3%, and all multitrait methods improved it (79.8%, 91.6% and 86.8% respectively for MV, DA and PCA). The computing times were similar to those of the case of 500 progeny.
Comparing power with α = 0.5 or 0.3 led to ratios of 2 to 3.5. But even when the power of ST was similar to the type-I error, multitrait methods could improve it if P 12 was negative. Ranking of the methods was similar to the one previously described when comparing the power for two traits: DA was more powerful than PCA, which was better than MV, which was more powerful than sPCA.
In Table III, the power for PCA came either from PCA1 or PCA2 depending on the simulated residual correlation. As previously shown, when the phenotypic correlation was positive, for ρ 12 = 0 to 0.6, the principal component allowing the detection of the QTL was PCA1. Otherwise, it was PCA2.

Accuracy of QTL position estimates
The average of the MSE for the QTL position estimates when α = 0.50 (0.30) was 0.066 (0.115) for ST (Tab. III). Ranking of the methods was similar to the one previously observed in Section 3.1 concerning the accuracy of position estimates. Trends depending on the QTL effects and the residual correlation were similar to those underlined for the power. The use of multitrait methods reduced both position bias and sampling variation compared to ST.

Accuracy of QTL effect estimates
The MSE for the QTL effect estimates are presented in Table III. For ST and MV, the MSE were independent of the correlation, respectively 0.0297 and 0.0270 on average when α = 0.50 and 0.0312 and 0.0287 on average when α = 0.30. The square of the bias represented for each of them 1, 0.91, 0.90 and 0.89% respectively. Concerning DA and PCA, the MSE increased when the residual correlation increased. For PCA, these MSE were lower than for the other methods, representing from 13 to 90% of the MSE of ST. The estimation was 10% more precise when α = 0.50 than when α = 0.30. For DA, when the residual correlation was negative, the MSE of the QTL effect estimates were equivalent to those of PCA. When the residual correlation was positive, they became higher. When α = 0.50, they were still lower than those of ST and MV, but when the effects were lower, the MSE of DA became larger than the others. As previously mentioned, the difference between DA and PCA came from a larger increase of the bias for DA with an increasing residual correlation than for PCA (not shown).

Using residual correlation to detect single trait QTL
A practically common use of the multitrait methods is presented in Table IV. The QTL only determined the second trait. The residual correlation between the two traits was 0.6 or −0.6. Also in such situations, using some of the multitrait methods improved the power compared to ST. The most powerful method was DA, improving the power by 40%. PCA2 allowed a better detection than PCA1, whatever the residual correlation. The power of PCA was lower than the power of the other multitrait methods. Indeed, since the QTL was no longer pleiotropic, the axis representing the variability due to the QTL was no longer the axis explaining the greatest phenotypic variability. In contrast to previous situations, sPCA power was better here than MV, due to the big excess of the number of parameters estimated in the latter method. We also compared the accuracy of the QTL position estimation (Tab. V). Results similar to those of power were obtained.
The big improvement of ST results in such cases seemed to be limited to the detection of QTL with a relatively high effect. These QTL must have already been detected by ST. But multitrait methods can state more precisely the parameter estimates. Table V synthesises the thresholds obtained depending on the marker density and the number of progeny per female. The most important factor was the number of progeny per female, since it directly modifies the number of parameters to estimate. Thresholds for marker densities equal to 12.5 and 50 cM are greatly different due to a various amount of information on the transmission pattern. As noted previously, corrected thresholds of ST and PCA are similar to the thresholds of DA, and thresholds for sPCA are equivalent to those of MV.

Thresholds
The comparison of the empirical and theoretical thresholds (χ 2 df ,5% ) for large samples (Tab. V) shows great differences between practical and asymptotical frameworks, even for sparse genetic maps. This underlines the difficulty to approximate such thresholds, increased by the mixture of different distributions in the likelihood for our study. Therefore, fast methods are required to estimate empirical thresholds accurately.

Power, accuracy of estimates of the QTL position, and accuracy of QTL effect estimates
As expected from single trait methods [18], power and mean squared error of the QTL position were improved by higher marker density and number of  progeny per female (Tab. VI). Improving ST QTL detection with low power using multitrait methods was efficient only if the genetic map was relatively dense or the residual and genetic correlations were in favourable relationships (not shown). These characteristics did not modify the ranking of the multitrait methods observed previously.
The number of progeny per female did not greatly influence the MSE of the QTL effect estimates. On the contrary, using a relatively dense genetic map reduced the MSE by a factor greater than 20 compared with a sparse map. In the case of a relatively sparse map, MV reduced the MSE of the QTL effect from 25 to 38% compared with ST. PCA reduced them by almost 30 to 50%. For DA it reduced the MSE by 20 to 65%. Once again, the greatest improvement of the MSE of the QTL effects was generally obtained from DA. But when the power was low, PCA was more accurate.

CONCLUSIONS
We compared three multitrait methods with a single trait method applied to a mixture of full and half sib families. One multitrait method was based on a multivariate analysis and was assumed to be greatly limited by its important number of parameters to estimate compared to the population size. The others rested on a univariate analysis of linear combinations of the traits.
Our results for the single trait method were in agreement with those of previous studies (see [10,11,18,27]) on half sib families or a mixture of full and half sib families. Particularly, such methods are precise for the estimation of QTL effects, but their estimation of the QTL position is frequently not very accurate, especially when the power is low.
The impact of the relationship between the residual and genetic correlations for each pair of traits has been previously reported when analysing inbred populations with MV and PCA [14,20,22,26]. It is also the major limit of the enhancement of the detection of a pleiotropic QTL by any multitrait method in our study. But the most practical situation consists in adding a trait with equivalent or lower effect from the QTL, eventually undetected by the ST methods. The effects and residual correlation due to the QTL are then generally unknown. The improvement from multitrait methods can then not be quantified before the analysis.
We confirmed the assumed limit of the multivariate method applied to a mixture of full and half sib families. Its high number of parameters greatly limits the power of detection and the accuracy of the position estimation of the method, and the computing time. This last point is mostly important when estimating thresholds through simulations or permutations, or when analysing big populations (more than 40 females seems unaffordable for systematic detection). We showed, however, that thresholds obtained from sPCA are The multitrait methods behaved similarly in terms of power and MSE of the position estimates. These trends were consistent with those of previous studies from inbred lines [14,20,22,26]. DA, since it is QTL-specific, is generally the most powerful and precise method for the estimation of QTL location. But it is not the most accurate for the QTL effect when the detection power is low. This seemed to be due to the correction factor applied to the MSE of QTL effects for DA and PCA in order to compare them to those of ST and MV. Indeed, when no QTL was detected, the coefficients of the linear combination for DA were very low. Thus, the correcting factor was small, increasing the MSE.
Moreover, to exploit the residual correlation between two traits to detect a single trait QTL, DA is also the best method. It makes it easy to select traits actually contributing to the detection of a QTL. Criteria to select traits could be comparing the maximum of the test statistics of DA with p and p − 1 traits with the χ 2 1df , or with some information criterion like the Akaike Information Criterion (AIC, [1]). Compared to the solution proposed by Korol et al. [22] to perform 10 000 permutations to test the significance of the contribution of each trait, it would save much time by limiting the number of implementations.
As for PCA and DA, the estimated effects cannot be directly transformed to effects from the QTL on each trait, we suggest a two-step procedure: (1) the use of DA (or PCA if the genetic structure is known to be simple, with only one QTL determining the traits for example) to rapidly identify the chromosomic segment of interest and the traits involved; and (2) the implementation of MV restricted to this small region to get estimates of effects and a residual covariance matrix.
In this study, the simulations were performed with fixed alleles in the grand parental breeds. The robustness of linear transformations based on the overall population variability has been confirmed elsewhere [6]. Moreover, we did not consider the problem of missing phenotypic data when performing multitrait analysis. One strategy could consist in replacing the missing data by the weighted mean of the phenotypes of individuals carrying the same haplotype at the studied position for the trait. Another could consist in removing the individual from the linear combination estimation and integrating it in the test statistic. But since no study has tested these proposals, it should be more secure to remove the individuals from the analysis, even if this means withdrawing a lot of progeny missing only one phenotype over all the traits.

APPENDIX
The data were assumed coming from a cross between grand parental breeds fixed for the QTL alleles. The determinism of the traits was assumed additive. We will demonstrate that linear combinations obtained from DA are identical when calculated from G = 2, 3 or 4 groups. For the G groups, the linear combination was obtained from the eigenvector z G associated to the greatest eigenvalue κ G of W −1 G B G , where W G is the residual covariance matrix (the within group covariance matrix), B G is the genetic covariance matrix (the between group covariance matrix). z G is also the eigenvector associated to the greatest eigenvalue of T −1 B G [28], where T is the total covariance matrix. Since T is independent of the number of groups, we will focus on B G . It can be estimated as G h=1 P h g h g h , where g h is the vector of the means of the performances of the group h and P h is the sum of the weights of the individuals of the group h, with G h=1 P h = 1. Haplotypes of progeny are considered as known, thus P h is the proportion of individuals in the group h. α l is the substitution effect of the QTL alleles on trait l. +α l corresponds to the effect of the QTL genotype 1/1, whereas −α l corresponds to 2/2. a = {α l ; l = 1, . . . , p}, and A = aa .
For two groups of haplotypes, transmitted by the sires, g q corresponds to the QTL allele q transmitted from the sire. g 1 = a/2 = −g 2 , which means that g 1 g 1 = g 2 g 2 = A/4. With P 1 = P 2 = 0.5, we then have B 2 = A/4.