QTL Mapping Using Multiple Markers Simultaneously

Methods for detecting and locating a single QTL within a 10 cM region of DNA using 10 equally spaced SNP markers were compared. The QTL was assumed to be bi-allelic and located between markers 5 and 6. Monte Carlo simulation of a granddaughter design with 30 sires and 400 sons was used. Linear regression nested within sire using either two or four marker haplotypes at a time was used. In addition, the scoring of haplotype transmissions from sire to sons were varied in three ways. Another method assumed linkage disequilibrium and estimated haplotype interval effects for all intervals simultaneously. Other variables compared were the ratio of QTL variance to total genetic variance, the number of generations of historical recombination, and frequencies of marker alleles. Empirical power was dependent on the scoring method in the linear regression method. Fourmarker haplotypes gave slightly higher empirical power than two-marker haplotypes. Reducing the proportion of QTL variance decreased empirical power. Empirical power was greater for 25 generations of historical recombination over 100 generations. Empirical power was lower when marker allele frequencies were 0.8 compared to 0.5. The linkage disequilibrium model gave results similar to those of the linear regression model.


INTRODUCTION
SNP are commonly used for detection and localization of QTL for complex traits. Several algorithms for identifying haplotypes using SNPs have been developed [6,10] . In practice, SNP genotypes will be available without phase information and the most frequent haplotypes need to be determined [1,8] . When there are enough offspring genotypes, as in half-sib designs in dairy cattle, then sire haplotypes can be reconstructed with near certainty. Determining the haplotype of the sire that is inherited by the son is difficult if the sire is heterozygous.
Detection and localization of QTL in livestock populations have used single markers, multiple markers, and haplotype-based methods analyzed with least squares or maximum likelihood procedures. Probabilities of QTL alleles being identical by descent (IBD) have been used with variance component estimation approaches [2,3] . Haplotype and linkage disequilibrium (LD) studies covering small chromosomal regions assume that sufficient generations have passed such that in small regions, the marker and QTL alleles tend to segregate together [5] .
Linear regression has been used with genotypes of single and multiple markers [4] . This approach uses linkage information and follows segregation of markers and QTL within sire families, but the use of haplotypes and linkage disequilibrium information could provide better precision in locating QTLs.
The objective of this study was to estimate QTL location within a group of markers using all marker haplotypes simultaneously. The method is assumes linkage disequilibrium and the needs to know the transmission of haplotypes from sire to sons. In a simulated granddaughter design the method was compared to linear regression with single and multiple marker methodology.

MATERIALS AND METHODS
The Simulation Method: A genomic structure of ten bi-allelic SNPs (each 1 cM from the next) with a single bi-allelic QTL located between the 5 th and 6 th markers was assumed. Initial allele frequencies, p, were either 0.5 or 0.8 for both markers and QTL. The base population was assumed to be in Hardy-Weinberg equilibrium at each locus and animals unrelated. Initial haplotypes for each sire and dam were generated independently. The number of sires and dams per generation were 30 and 400, respectively, giving an effective population size of 100. Progeny were generated using the Haldane mapping function (no interference) to determine crossover events between markers and QTL in deriving sire and dam haplotypes. Each mating produced 4 progeny. Each generation was discrete with new males and females randomly selected from the group of progeny.
Linkage disequilibrium was created over either 25 or 100 generations of random matings. In the last generation, 10 males were randomly selected and mated to 1000 females to produce 100 sons per sire. Each son was assumed to have 100 progeny. Progeny Yield Deviations (DYD) were created for each son. A DYD is the average phenotype of the progeny adjusted for systematic environmental effects and the genetic values of the dams [7]. e BV DYD son + = 5 . 0 (1) where BV son is the breeding value of the son, and e is a residual effect with mean zero and variance equal to residual variance divided by number of offspring (100 in this study). The trait was assumed to have a heritability of 0.3 and σ P =100.
Parameters allowed to vary were the ratio of QTL variance to total additive genetic variance (k=0.1 or 0.05), the allele frequencies of the SNPs and QTL (p=0.5 or 0.8), and the number of generations of historical recombination (g=100 or 25).

Linear Regression Method: The model was
where y ij is the DYD of son j of sire i, x ij indicates the transmission score of the haplotype from sire i to son j, β i is the within sire family regression coefficient, s i is the polygenic effect of the sire, and e ij is the residual effect.
The assumption was made that only sires and sons were genotyped for the markers. The scores for x ij normally used by Knott et al. (1996) are shown in Table 1. Note that progeny of homozygous sires are not used. This method uses linkage information and segregation of marker alleles within heterozygous families. In this study, two different sets of transmission scores were used and included progeny of both heterozygous and homozygous sires. The two sets of scores differed from each other by two numbers as shown in Table 1. When a sire is homozygous, then his progeny have the same haplotype, and so the transmission score in Set 1 was set to 1. With Set 2, the transmission score was set to 2 for homozygous sires with homozygous sons. Haplotypes were created using either two contiguous markers at a time (9 sets of two) denoted as method LR-2, or four contiguous markers at a time (7 sets of four) denoted as method LR-4. The analyses proceeded from one marker set to the next, and F-ratios were used to determine the location of the QTL. Transmission scores were applied to the haplotype of the set of markers.

Simultaneous
where y ijkl is the DYD of the l th son of sire k, µ is the overall mean, h ij is the j th haplotype of i th contiguous pair (i=1 to 9 pairs, j=1 to 4 possible haplotypes within pairs, 36 effects), s k is the polygenic effect of the sire, e ijkl is a residual effect, 2 e σ and 2 h σ were the same for all marker set. Haplotypes were constructed based on each pair of contiguous markers, giving 9 nonoverlapping haplotype intervals. For each interval there were 4 possible haplotypes that could be inherited from sire to son, and therefore, a total of 36 haplotype effects. If the haplotype coming from the sire could not be determined, the design matrix contained all zeros for that interval for that son. Otherwise there was a one corresponding to the haplotype that was inherited from the sire in the design matrix. In matrix notation, the model is e s h y + + + = µ (4) and the mixed model equations are To locate the QTL, the sum of absolute values of estimated haplotype effects within each haplotype pair was calculated, The middle of the haplotype interval with the largest H i was taken as the location of the QTL, and the ij ĥ within the largest H i indicated the magnitude of the QTL effects.
The ratio of residual to haplotype variances was fixed at 1.0 for all pairs of haplotypes. However, the ratios could differ between pairs, such that the ratio would be larger for pairs of markers that do not contain a QTL. Thus, simulations were repeated where the ratio for marker pairs without the QTL was 10, and for the marker pair containing the QTL was one. No attempts were made to estimate the ratio in this study.
The number of marker pairs could also have an effect on QTL detection, and so the simulations were repeated using 20 markers (with all variance ratios equal to 1) within a 10 cM chromosome segment. The QTL was placed in the middle of the markers. Thus, there were 76 haplotype effects to be estimated in this scenario.
Comparison Statistics: Let Q be the true position of a QTL and Q be the estimated position. The estimated position was the mid point of the marker haplotype interval that gave the largest F-ratio or the largest H. Bias of estimation was where n is the number of replicates, and the mean squared error was (9) The critical values of the test statistic corresponding to 5% type I error were determined empirically by simulating 1000 replicates based on the null hypothesis of no QTL segregating on the chromosome segment. Another 1000 replicates were simulated under the alternative hypothesis, with the QTL situated in the center of the markers, and the empirical power was computed as the proportion of replicates in which the test statistic exceeded the critical value.
Scenarios: Four scenarios of parameter combinations were studied. Scenario 1 had the ratio of QTL variance to total genetic variance of 0.1, used 100 generations of historic recombination, and started with an allele frequency of 0.5 for all markers and the QTL. Scenario 2 had the ratio of QTL variance to total genetic variance of 0.05, and other parameters were the same as Scenario 1. Scenario 3 was the same as Scenario 1 with the number of generations of historic recombination reduced to 25. Scenario 4 was the same as Scenario 1 with the allele frequency of markers and QTL increased to 0.8.

RESULTS
Linear Regression Two-marker Intervals: Under Scenario 1 (the base combination of parameters), using Set 1 transmission scores, empirical power was 0.77, and using Set 2 the power was 0.82 ( Table 2). The only difference was in the scores assigned to a homozygous son from a homozygous sire. Precision was better with Set 2 (1.33 cM versus 1.37 cM for Set 1). However, the probability of a type I error was smaller with Set 1.
Results are shown in Fig. 1 for the two sets of indicator variables and for the Knott et al. (1996) scores for comparison. With scores of Knott et al. (1996), the QTL position could not be detected because the average of 1000 replicates of F-values was below the critical value (1.83).
In Scenario 2, the QTL variance was only half as large as in Scenario 1, and this reduced power, decreased precision, and increased the probability of type I error, but not significantly compared to Scenario 1. Set 2 scores gave better results in this scenario also. Smaller QTL effects gave generally poorer results for both sets of scores.
In Scenario 3, the number of generations of historical recombination was changed from 100 to 25. Empirical power increased from 0.77 to 0.95 with Set 1 indicator variables and from 0.82 to 0.98 with Set 2. Precision was below 1 cM for both sets of scores. Probabilities of type I errors were only 0.04 and 0.10 for Sets 1 and 2, respectively. With 100 generations of random mating, there was more chance for the marker (and QTL) alleles to drift to fixation, than within 25 generations. This might explain the differences observed.
In Scenario 4, allele frequencies of markers and QTLs were changed from 0.5 to 0.8. Empirical power decreased from 0.77 to 0.54 and 0.50 for Sets 1 and 2, respectively, and the QTL location was less precise at over 1.35 cM, but the probability of type I error was very low. Bias and MSE were much greater than the other scenarios. Set 2 scores had slightly smaller Bias and MSE than Set 1, but lower power and poorer precision.

Linear Regression Four-marker Intervals:
When four markers were used to establish haplotypes for the linear regression model, results are given in Table 3. Degrees of freedom were expected to be less because of fewer intervals being considered. Empirical power was the same or greater than with the two-marker intervals, but the probabilities of type I errors were slightly greater. Precision was better in all scenarios with four markers rather than two.
Comparison of scenarios followed the same pattern as for the two-marker intervals. However, the precision, bias, and MSE in scenario 4 when the allele frequencies changed from 0.5 to 0.8 were better with the fourmarker intervals than with the two-marker intervals. The best scenario was 3 having only 25 generations of historical recombination. Hi Hx 1 0.5 0      (Fig. 2). Power for this model was similar to the linear regression model using four marker intervals and Set 1 scores. Probability of type I errors was similar to the linear regression model with two marker intervals and Set 1 scores, except for Scenario 4 where the simultaneous haplotype model was higher. Bias and MSE were greater than either linear regression approach and precision was about the same as the linear regression model with two marker intervals, except for Scenario 3, which gave poorer precision with the simultaneous haplotype model. The simultaneous haplotype model has the advantage that all marker intervals are considered at one time, compared to linear regression where only one marker interval at a time is analyzed. At the same time, this is a disadvantage in that more effects are being estimated at one time. The net result is that the simultaneous haplotype model does not appear to be better than the linear regression model.

DISCUSSION
The regression model of Knott et al. (1996) was not successful in detecting QTL within small regions of 1 to 5 cM. The method works better for large regions of 10 to 20 cM. However, by regression of transmitted haplotypes instead of genotypes, as in this study, and assuming linkage disequilibrium (LD), QTL can be detected within 1 cM regions, if the QTL exists, with considerable power and precision.
In this study the positions and order of the markers need to be known. In practice, the markers will likely not be equally spaced. One possibility is to use SNP genotypes rather than haplotypes assuming linkage disequilibrium exists, and to put 5 or more SNP markers within each cM segment, such that the markers are very near or even within a QTL. The SNP genotypes then approximate the QTL genotypes. The order of the SNPs would not need to be known, nor the distances. The SNP location giving the largest estimated effects would most likely be the location closest to the QTL.
The size of QTL effect is important in accurate mapping as expected, and also the allele frequencies for the QTL. Decreases in empirical power in this study were non significant in going from 10% to 5% of total genetic variance. The decrease in power was much less than in using the simple regression method or variance component method as in [3,4].
Empirical power of QTL detection increased when the number of generations of historical recombination decreased from 100 to 25. After 100 generations of recombination in a small population, most loci would have drifted to become homozygous and the informativeness of the haplotypes would decrease. Thus, effective population size may have some importance on QTL detection. With 25 generations of recombination, less homozygosity of loci was allowed to occur and the informativeness of haplotypes was consequently greater. There was no allowance for mutations to occur in order to maintain QTL genetic variability.
The most important aspect of the regression methods in this study was the scores used in the regression model. A preliminary study comparing seven different sets of scores was made. The better two from that study were used here. Sets 1 and 2 were not significantly different in performance, but Set 2 gave slightly greater empirical power and better precision. Set 1 was very close to conditional probabilities of inheritance of the sire haplotype, dependent on whether the sire was homozygous or heterozygous for the haplotype. With the Knott et al. (1996) method, only one set of scores was given whether the sire was heterozygous or homozygous for the haplotype alleles. There could be an optimum set of scores that could be derived, which might be dependent on the frequencies of different haplotypes.
The simultaneous haplotype model assumes the existence of linkage disequilibrium (LD), and many generations of historic recombination. LD can also be created by selection, but all replacements and matings were random in this simulation. The haplotype interval effects were assumed to be random effects, and in this study a common ratio of residual to QTL variance of 1 was used for all intervals. A comparison was made where the variance ratio was set to 1 for the interval with the highest H i , and set to 10 for all other intervals (Table 5). Compared to using a ratio of 1 for all intervals, power was lower and bias was slightly greater, while precision was slightly better. The differences were small. Thus, assuming the same ratio for all intervals is easier to apply, and estimation of the appropriate variances by Bayesian or restricted maximum likelihood is not necessary.
A variance ratio of 1 means that the variance for each interval is equal to the residual variance, but in fact, the ratio should probably be very large because the amount of variance accounted for by the QTL is only 0.10 or 0.05 of the additive genetic variance. For a heritability of 0.30 and QTL effects of 10%, the ratio of residual to QTL variance would be 23.33. Thus, using a ratio of 1 is only a small step up from assuming the interval haplotype effects are fixed. The purpose of using a ratio of 1 is to remove dependencies among the possible confounding of interval haplotype effects.  Knott et al. (1996) method with linear haplotype regression model using two sets of transmission scores for two-marker intervals SNP markers could be located more closely than every 1cM and presumably the precision of estimating QTL location should be better. Another comparison was made using 20 markers along a 10 cM segment of DNA with one QTL ( Table 6). Power of detecting QTL was greater than with 10 markers. More SNP markers at closer distances were able to locate the QTL better than having SNPs only every 1 cM apart. Denser sets of markers may be optimum. Fine mapping of QTL using linear regression on haplotypes or a simultaneous haplotype model is possible using closely linked markers in a small region of a chromosome. Further studies should examine the effects of QTLs having more than two alleles, and having more than one QTLs in a small region at one time with these new methods. The simultaneous haplotype model could be extended to cover the entire length of a chromosome. The model should be able to detect multiple QTLs of varying magnitudes of effects. Xu (2003) describes a similar model for within family analyses where regressions are on marker genotypes (as random variables) for all markers on a chromosome. Using Bayesian methods, Xu (2003) was able to pinpoint the locations of QTL without a lot of noise.