A Rank-Based Nonparametric Method for Mapping Quantitative Trait Loci in Outbred Half-Sib Pedigrees: Application to Milk Production in a Granddaughter Design

Corresponding author: Michel Georges, Department of Genetics, Faculty of Veterinary Medicine, University of Liège (B43), 20 Bd de Colonster, 4000 Liège, Belgium, michel{at}stat.fmv.ulg.ac.be (E-mail).

Communicating editor: Z-B. ZENG

	ABSTRACT

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

We describe the development of a multipoint nonparametric quantitative trait loci mapping method based on the Wilcoxon rank-sum test applicable to outbred half-sib pedigrees. The method has been evaluated on a simulated dataset and its efficiency compared with interval mapping by using regression. It was shown that the rank-based approach is slightly inferior to regression when the residual variance is homoscedastic normal; however, in three out of four other scenarios envisaged, i.e., residual variance heteroscedastic normal, homoscedastic skewed, and homoscedastic positively kurtosed, the latter outperforms the former one. Both methods were applied to a real data set analyzing the effect of bovine chromosome 6 on milk yield and composition by using a 125-cM map comprising 15 microsatellites and a granddaughter design counting 1158 Holstein-Friesian sires.

RECENT developments in DNA marker technology, such as the discovery of microsatellites (WEBER and MAY 1989 ), random amplified polymorphic DNA (RAPDs; WILLIAMS et al. 1990 ), and amplified fragment length polymorphism (AFLPs; VOS et al. 1995 ) as abundant sources of well-dispersed genetic markers, have boosted the construction of marker maps across a broad taxonomic range. Not only are such maps now available for human and model organisms such as mouse and rat but for a number of agriculturally important animal and plant species as well.

These maps increasingly are applied to locate genes underlying inheritable phenotypes of interest. Several of the most relevant phenotypes are continuously distributed quantitative traits involving multiple polygenes or quantitative trait loci (QTL), as well as nongenetic effects. Experimental back- and intercrosses are often the preferred design to map QTL. However, in a number of agriculturally important species (notably cattle and pine trees), reproductive cycles and breeding designs have led to the generation of extensive half-sib pedigrees that are readily available for QTL mapping. A well-documented example of this is the so-called granddaughter design to map genes underlying milk production in commercial cattle populations (WELLER et al. 1990 ). This design takes advantage of the numerous paternal half-brother pedigrees that exist in dairy cattle populations, generated as part of the applied progeny-test breeding design.

A number of mapping methods have been applied to such half-sib designs, including single-marker regression (e.g., COWAN et al. 1990 ), interval mapping using regression (e.g., KNOTT et al. 1996 ), and maximum likelihood methods (e.g., GEORGES et al. 1995 ). All these methods share a common assumption, namely the residual normal distribution and homoscedasiticity of the analyzed phenotypes or transformations thereof. These approaches therefore are not suitable for phenotypes that are known not to satisfy this normality assumption. Moreover, deviations from normality for traits that generally are assumed to be quasi-normally distributed are likely to affect the power and robustness of these conventional approaches.

Recently, KRUGLYAK and LANDER 1995A described a nonparametric QTL interval mapping approach based on the Wilcoxon rank-sum test applicable in experimental crosses. This method provided a robust alternative to conventional approaches, applicable to normally distributed traits with minimal loss of power and extending the scope of QTL mapping to a variety of traits not normally distributed, such as counts generated by a Poisson process, truncated data, probabilities, and qualitative data.

In this article, we describe the adaptation of this method to half-sib pedigrees in outbred populations and apply it to milk production in a granddaughter design. A computer program to implement this approach has been developed and is available from the authors upon request.

	MATERIALS AND METHODS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

A QTL interval mapping procedure based on the Wilcoxon rank-sum test—general principles:
To measure the evidence in favor of a QTL at a given map position, KRUGLYAK and LANDER 1995A define the following statistic (illustrated for an (A x B) x A backcross),

(1)

(2)

Under the null hypothesis of no QTL, Z_W is shown to behave asymptotically as a standard normal variable that reduces to a Wilcoxon rank-sum test at the marker positions.

Adaptation to outbred half-sib designs:
The method developed by KRUGLYAK and LANDER 1995A for experimental crosses was adapted to outbred half-sib designs, e.g., a founder sire mated to several dams to produce a large paternal half-sibship. The approach relies on the same Z_W(s) statistic. However, P[g_i,A(s)|g_i,L,g_i,R] (Equation 2) is now defined as the probability that progeny i has inherited QTL allele A from the founder sire—assumed to be heterozygous AB at the QTL—at map position (s) given its genotype at the left (g_i,L) and right (g_i,R) flanking markers. Only markers for which the founder sire is heterozygous are considered when computing P[g_i,A(s)|g_i,L,g_i,R]. Moreover, while the nearest flanking markers contain all information needed to compute P[g_i,A(s)|g_i,L,g_i,R] in a given interval when dealing with experimental crosses, information from more distant markers is considered in the outbred half-sib situation, when closer markers are not fully informative. This occurs in the case of missing genotypes or when the offspring has the same marker genotype as the sire, and the dam is either not genotyped or has the same heterozygous genotype as well. In the former case, part of the information can be recovered by considering marker allele frequencies in the population.

Calculation of P[g_i,A(s)|g_i,L,g_i,R] requires knowledge of the sire's marker linkage phase. In the absence of grandparental marker information, the most likely linkage phase is first estimated from the marker genotypes of the offspring. This is accomplished by calculating the likelihood of the pedigree data under the 2^x/2 possible phases (assuming x informative markers) as follows (GEORGES et al. 1995 ):

(3)

₁

_x

All marker phases are a priori considered to be equally likely; i.e., linkage equilibrium is assumed to be reached between all markers. The marker phase maximizing the likelihood of the pedigree data is considered the true one and is selected for further analysis.

As pointed out by KRUGLYAK and LANDER 1995A ,

(4)

While [1 - 2·P[g_i,A(s)|g_i,L,g_i,R]]² or the expected value of [1 - 2·P[g_i,A(s)|g_i,L,g_i,R]]² over all possible genotypes is computed easily for experimental crosses, its calculation is more cumbersome in outbred designs as it will depend on marker allele frequencies and genotype of the founder sire. The value of [1 - 2·P[g_i,A(s)|g_i,L,g_i,R]]² is therefore calculated for each half-sib pedigree by simulating all possible offspring and calculating a frequency weighted mean of [1 - 2·P[g_i,A (s)|g_i,L,g_i,R]]².

Across family analysis:
In practice, the available pedigree material is composed most often not of one half-sib pedigree but of a series of such half-sibships, such as in the granddaughter design (WELLER et al. 1990 ). In outbred populations, however, the different sibships cannot be assumed to segregate for the same QTL or even QTL alleles; i.e., one cannot assume locus and allelic homogeneity across families.

Rather than analyze the pedigrees separately, however, and reduce power by multiple testing, the individual Z_W(s) scores were squared and summed over all k families yielding a ² statistic with k degrees of freedom:

(5)

Interval mapping by regression:
The rank-sum-based approach (hereafter referred to as method RS) was compared with interval mapping by using regression (hereafter referred to as method MR for multipoint regression; KNOTT et al. 1996 ). For each half-sib family, j, phenotypes were regressed on P[g_i,A(s)|g_i,L,g_i,R], calculated as described above, yielding least-squares estimators of the y intercept, ß_0j, and the slope, ß_1j, the latter being an estimator of the QTL allele substitution effect in the corresponding family, j. The ratio

Significance thresholds:
For both the RS and MR methods, chromosome-wise significance thresholds were determined from the distribution of the test statistic over 10,000 permutations (simulated data set) or 100,000 permutations (real data set) of the phenotypes (or ranks) as suggested by CHURCHILL and DOERGE 1995 . Phenotypes were permutated within family. For each permutation, the highest value of the test statistic over the entire chromosome was retained to yield "chromosome-wise" distributions of the test statistic under the null hypothesis. For the real data set, a Bonferonni correction was applied to the chromosome-wise significance level, considering that chromosome 6 represents 1/29 of the bovine autosomes and that we analyzed the equivalent of three independent traits (SPELMAN et al. 1996 ) to obtain "experiment-wise" significance thresholds.

Simulated data set:
To test the efficacy of the proposed method, we simulated the segregation of a QTL in a granddaughter design. The pedigree material was composed of two paternal half-sib families of 100 sons, four families of 50 sons, and eight families of 25 sons, quite accurately reflecting a real data set. The 14 respective sires were considered to be unrelated.

A QTL was positioned in the center of the fourth interval of a map comprising seven markers spaced 15 recombination units apart. Markers were assumed to be polyallic markers with frequencies randomly assigned from a uniform distribution and rescaled to sum to unity, yielding a heterozygosity of

The QTL was assumed to be biallelic with frequencies p = 0.25 (Q) and q = 0.75 (q), respectively. Founder-sires therefore had an a priori probability 2pq = 0.375 to be heterozygous Qq for the QTL. Following Falconer's notation (FALCONER and MACKAY 1996) and assuming additively acting alleles, the average phenotypic values of the QQ, Qq, and qq genotypic classes were set at +a, d = 0, and -a, respectively. Assuming Hardy-Weinberg equilibrium, this yields an average effect of an allele substitution, = a, and a variance attributable to the segregation of the QTL:

The value of a was determined such that

h² was set at 9.4% for all simulations, corresponding to an a value of 0.5_P. Five scenarios were considered to model the residual variance, ²_R: (1) homoscedastic, normal residual variance, (2) heteroscedastic, normal residual variance, (3) homoscedastic, skewed, or asymetric residual variance, (4) homoscedastic, positive kurtosis or more peaked around the center than the density of the normal curve, and (5) homoscedastic, negative kurtosis or flatter around the center than the density of the normal curve.

Homoscedastic normal residual variance: Individual phenotypic values were generated as the mean of the genotypic class to which the individual belongs (QQ = a, Qq = 0, or qq = -a) plus a value drawn from a normal distribution with mean 0 and variance 1; i.e., ²_R was set at one.

Heteroscedastic normal residual variance: Individual phenotypic values were generated as the mean of the genotypic class to which the individual belongs (QQ = a, Qq = 0, or qq = -a) plus a value drawn from normal distributions with mean 0 and variances of ²_R(QQ) = 1, ²_R(Qq) = r , and ²_R(qq) = s , such that

Homoscedastic, skewed residual variance: Skewness was simulated by assuming a residual effect distributed as a chi-squared distribution with n degrees of freedom, with variance ²_R = 2n and mean n. Individual phenotypic values were generated as the mean of the genotypic class to which the individual belongs (QQ = a - n, Qq = 0 - n, or qq = -a - n) plus a value drawn from a chi-squared distribution with n degrees of freedom, obtained by summing n squared values drawn from a standard normal.

Homoscedastic, positive kurtosis: Excess of kurtosis was simulated by assuming that the residual effect was distributed as a Student's t-distribution with n degrees of freedom, with variance ²_R = and mean 0. Individual phenotypic values were generated as the mean of the genotypic class to which the individual belongs (QQ = a, Qq = 0, or qq = -a) plus a value drawn from a t-distribution with n degrees of freedom, i.e.,

Homoscedastic, negative kurtosis: Negative kurtosis was simulated by assuming that the residual effect was distributed as a hemicircular distribution with mean 0 and variance ²_R = , where r is the radius of the hemicircle. Individual phenotypic values were generated as the mean of the genotypic class to which the individual belongs (QQ = a, Qq = 0, or qq = -a) plus a value drawn from this hemicircular distribution. This was done by determining the value of t such that

Figure 1 illustrates the expected phenotypic distributions of offspring from heterozygous founder-sires, Qq, for the five examined models. Offspring are sorted in two genotypic classes depending on the QTL allele transmitted by the sire (Q or q). Each class therefore comprises two subpopulations: QQ (25%) and Qq (75%) for the Q class and Qq (25%) and qq (75%) for the q class.

View larger version (26K): In this window In a new window Download PPT slide   Figure 1. Phenotypic distributions of offspring from heterozygous Qq sires, sorted according to the QTL allele inherited from the sire (Q or q), assuming (a) a heteroscedastic normal residual variance (r = 2; s = 4); (b) a homoscedastic, skewed residual variance (28 ); (c) a homoscedastic, positively kurtosed residual variance (t5); and (d) a homoscedastic, negatively kurtosed residual variance (hemicircular residual variance). The phenotypic distributions of the q? offspring are shown (thick black lines) and compared with the corresponding distribution assuming a homoscedastic normal residual variance (thick gray lines). The corresponding distributions of the Q? offspring are shown as thin lines. Each class therefore comprises two subpopulations: QQ (25%) and Qq (75%) for the Q? class and Qq (25%) and qq (75%) for the q? class. The differences between the means of the Q? and q? populations, corresponding to the effect of the Q to q allele substitution, equal 0.5 P.

At least 200 datasets (ranging from 200 to 866) were simulated under each of the five models of residual variation and analyzed with the RS and MR methods. Permutations were used to estimate the significance levels reached for each of these analyses (CHURCHILL and DOERGE 1995 ). For each replicate, 10,000 permutations were performed and analyzed with the RS and MR methods to yield a dataset-specific, chromosome-wise distribution of the RS and MR statistics under the null-hypothesis, allowing us to measure the P value of the unpermutated data under the null hypothesis of no QTL. Average P values over replicates were calculated for each of the five models. For each model, the proportion of datasets yielding a P value less than 0.05 (=) was used to measure the corresponding power (1 - ß) of the RS and MR methods (Table 1). Within each model, we compared the relative merits of the RS vs. MR methods by applying the Wilcoxon matched pairs test on all resulting pairs of P values (HOLLANDER and WOLFE 1973 ). Within each method, the effect of the model on the power to detect the QTL was evaluated by using the Mann-Whitney U test (HOLLANDER and WOLFE 1973 ), using model 1 as reference.

Real data set:
The real data set was a Holstein-Friesian granddaughter design comprising 1158 sons distributed over 29 paternal half-sib families, partially described in SPELMAN et al. 1996 . The number of sons per family ranged from 11 to 153.

All animals were genotyped for a battery of 15 previously described (KAPPES et al. 1997 ) microsatellite markers from bovine chromosome 6 (Table 2). Genotyping was performed as described (GEORGES et al. 1995 ) or by using the "four dye-one lane" technology on an ABI373 or ABI377 sequencer.

Marker maps were built by using the CRIMAP program (LANDER and GREEN 1987 ) to determine the most likely order and the ANIMAP program to refine the most likely recombination rates between adjacent markers (GEORGES et al. 1995 ).

Information content along the marker map (KRUGLYAK and LANDER 1995B ) was measured as

QTL mapping was performed for five milk production traits: milk yield (kg), protein yield (kg), fat yield (kg), protein percentage, and fat percentage. The phenotypes used for QTL mapping were deviations of individual daughter yield deviations from the corresponding average of the parental predicted transmitting abilities (VAN RADEN and WIGGANS 1991 ).

Marker allele frequencies, required for map construction, measuring of information content, and QTL mapping, were estimated from the dam population, separately for each pedigree, as

p1 and p2 correspond to the frequencies of two alleles from the sire, while p3 is the frequency of all other alleles pooled. n_xy corresponds to the number of sons in the pedigree with genotype xy, and n equals the total number of sons in the pedigree.

	RESULTS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Simulated data:
Using the approach described in MATERIALS AND METHODS, we simulated GDDs segregating for a QTL explaining a fixed 9.4% of the phenotypic variance (corresponding to a = 0.5_P) but with five alternative residual components: homoscedastic normal, heteroscedastic normal, homoscedastic skewed, homoscedastic positive kurtosis, and homoscedastic negative kurtosis. The generated datasets were all analyzed by using both RS and MR methods. Table 1 reports, for each of the five scenarios, the average P values and the associated power at -value of 0.05, obtained by permutation as described in MATERIALS AND METHODS.

The relative merit of the RS and MR methods was evaluated by using the Wilcoxon matched pairs test as described in MATERIALS AND METHODS. As expected, multiple regression is superior to the rank-sum approach under the basic model of homoscedastic normal residual variance (P = 0.000014). The loss of power when using the rank-based method is estimated at 8% at -value of 0.05. The MR method proved also significantly superior to the RS method in the negative kurtosis model (model 5; P = 0.000001); the loss of power with the RS method was estimated at 14% at -value of 0.05. For the three remaining scenarios, however, the RS approach outperformed MR, the gain in power ranging from 8 to 20% at -value of 0.05 (Table 1).

The effect of the model on the power to detect the QTL was evaluated by using the Mann-Whitney U test (see MATERIALS AND METHODS), by using model 1 as reference. Comparisons were performed separately for the RS and MR approach. Interestingly, MR appears to be quite insensitive to the nonnormality of the residual variation, as the distribution of P values under the alternative models is never significantly different from those obtained under the basic model. This is likely due to the fact that significance levels are deduced from phenotype permutations rather than from the theoretical distribution of the test statistic. Using RS, on the contrary, significant increases in detection power are observed for models 2, 3, and 4 (respectively 9, 12, and 23% at -value of 0.05; Table 1), while the distribution of P values does not differ significantly between models 1 and 5.

Estimates of the precision in the estimation of QTL positions were also compared. Table 1 shows the standard deviation of the most likely QTL position for all simulations yielding a signal exceeding the 5% chromosome-wise significant threshold. Comparing the difference between real and estimated position by using the Mann-Whitney U test, we found no evidence for a significant effect either of the statistical method or of the model for the underlying residual variance. In essence, precision was as poor in all circumstances, standard deviations of the estimated position being 20 to 25 cM. While the actual position of the QTL was at 62 cM counting from the first marker, the estimates ranged from 0 to 118 cM, i.e., the entire chromosome length. A total of 95% of the estimates were within 43 cM (=1.9 ) from the actual position.

Real data:
Table 2 and Figure 2 show the most likely marker map as obtained from our genotypes. The map covers 125 cM (Kosambi) with average interval of 9 cM. The most likely order was in agreement with KAPPES et al. (1996). The same figure also compares information content when (1) exploiting marker allele frequency estimates to extract information from noninformative marker genotypes, and (2) when ignoring this information, i.e., when considering all microsatellite alleles to be equally frequent in the population. It can be seen that more than 80% of the maximal information is extracted for the central part of the chromosome; however, the information content drops at both extremities of the chromosome. Moreover, the figure shows that information content is improved only marginally by considering marker allele frequencies. This is especially true in the central, denser part of the marker map.

View larger version (27K): In this window In a new window Download PPT slide   Figure 2. Information content (in percentage of the theoretical maximum) map along the length of bovine chromosome 6 when using (+++) or ignoring (---) marker allele frequencies. Marker names and corresponding map position in centimorgans are shown along the x-axis.

Figure 3A and Figure 3B summarize the location score profiles obtained for the five different milk production traits by using both RS and MR approaches. Generally speaking, both methods clearly yield very similar curves for all traits along the entire chromosome length. For protein percentage, the location scores maximize at the same chromosome position (48 cM) using both approaches. The associated experiment-wise significance levels are P = 0.03 for RS and P = 0.01 for MR, therefore slightly superior for the latter.

View larger version (34K): In this window In a new window Download PPT slide   Figure 3. Location scores obtained along chromosome 6 for milk (), fat (), and protein () yield, as well as fat () and protein () percentage, using the RS (a) and MR methods (b). The y-axis corresponds to the log10 of the inverse of the corresponding chromosome-wise P value as determined by permutation. Horizontal bars on the graphs correspond to 15, 10, 5, and 1% experiment-wise thresholds, obtained by applying a Bonferroni correction to the chromosome-wise significance levels.

These results are in agreement with the report of a QTL affecting milk production on the centromeric half of chromosome 6, first identified by GEORGES et al. 1995 and later confirmed in independent studies in Holstein-Friesian by SPELMAN et al. 1996 and KUHN et al. 1996 , in Finnish Ayrshire by VILKKI et al. 1997 , and in Norwegian Red by GOMEZ-RAYA et al. 1996 . A detailed analysis of this chromosome region in the corresponding pedigree material is given in SPELMAN et al. 1996 .

	DISCUSSION

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

In this article, we have adapted a nonparametric QTL mapping method based on sum of ranks that was described previously for experimental crosses (KRUGLYAK and LANDER 1995A ) to outbred half-sib pedigrees. This is particularly relevant for mapping QTL in specific livestock and plant species where such pedigrees routinely are generated within the context of specific breeding designs. It extends the scope of QTL mapping in these pedigrees to a variety of not normally distributed traits, including counts generated by a Poisson process, truncated data, and probabilities and qualitative data (KRUGLYAK and LANDER 1995A ).

We confirm that this approach (the RS method) can be applied conveniently to normally distributed traits with minimal loss of power when compared to parametric methods. In the simulated example, we noticed a loss of power of 8% at -value of 5% when compared to the MR method. When simulating nonnormal or heteroscedastic residuals, however, the RS method outperformed the MR method in three out of four scenarios (models 2–4: heteroscedastic normal, homoscedastic skewed, and homoscedastic positively kurtosed). Interestingly, this was shown not to be due to a loss of power of the MR approach, which proved to be relatively robust in the scenarios that we simulated, but rather to a gain of power when applying the RS method. Our interpretation of this finding is that in the three scenarios where RS proved superior to MS, the phenotypic distribution is characterized by "outlyers" when compared to the normal distribution (see Figure 1). These outlyers contribute excessively to the residual variation, while the bulk of the observations actually are more centered around the mean (and therefore less variable) when compared to the normal distribution. When using ranks rather than the actual phenotypes, the contribution of the outlyers to the residual variation is tempered, therefore increasing the ratio QTL variance/residual variance and concommitantly increasing the power to detect the QTL.

A disadvantage of the rank-based methods is the fact that these do not provide convenient estimates of QTL effects. These methods therefore are suitable for the detection of QTLs but have to be complemented with alternative methods, such as least-squares or maximum likelihood techniques when quantifying the QTL effects.

Recently, a number of QTL mapping methods that account for multiple linked or unlinked QTL have been proposed. These include two QTL models (e.g., HALEY and KNOTT 1992 ), composite interval mapping (ZENG 1993 ), and multiple QTL mapping (JANSEN 1993 ). Rank-based approaches have been described to test three or more classes, including the Kruskal-Wallis test and the Jonckheere-Terpstra test, which would allow a two-QTL model to fit. Alternatively, it might be interesting to explore the possibility to use regression techniques directly on ranks, which, if applicable, would allow inclusion of additional markers as cofactors in the model.

Assuming paternal half-sib pedigrees, the proposed method allows for missing genotypes in the "dams." In such cases, estimates of marker allele frequencies can be used to improve inference about the identity of the transmitted paternal chromosome. However, it is shown that when performing multipoint analyses with dense marker maps, this contributes only a marginal improvement of the information content. The benefit of including marker allele frequency is therefore doubtful. Indeed, errors in the estimation of the marker allele frequencies may even cause an increase in type I errors or a loss of power if accounting for inaccuracies in the frequency estimates (CHARLIER et al. 1996 ).

As expected, the precision in the estimation of the QTL position using both proposed parametric and nonparametric approaches is mediocre. This illustrates the need to develop alternative strategies for fine-mapping QTL in outbred populations.

	ACKNOWLEDGMENTS

We acknowledge the financial support of Holland Genetics, Livestock Improvement Corporation, the Vlaamse Rundvee Vereniging, and the Ministère des Classes Moyennes et de l'Agriculture, Belgium. Continuous support from NANKE DEN DAAS, BRIAN WICKHAM, DENIS VOLCKAERT, and PASCAL LEROY is greatly appreciated. We thank JOHAN VAN ARENDONK, RICHARD SPELMAN, HENK BOVENHUIS, MARCO BINK, DAVE JOHNSON, and DORIAN GARRICK for fruitful discussions.

Manuscript received November 6, 1997; Accepted for publication March 30, 1998.

	LITERATURE CITED

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

CHARLIER, C., F. FARNIR, P. BERZI, P. VANMANSHOVEN, and B. BROUWERS et al., 1996  IBD mapping of recessive traits in livestock: application to map the bovine syndactyly locus to chromosome 15. Genome Res. 6:580-589[Abstract/Free Full Text].

CHURCHILL, G. A. and R. W. DOERGE, 1995  Empirical threshold values for quantitative trait mapping. Genetics 138:963-971[Abstract].

COWAN, C. W., M. R. DENTINE, R. L. AX, and L. A. SCHULER, 1990  Structural variation around prolactin gene linked to quantitative traits in elite Holstein sire family. Theor. Appl. Genet. 79:577-582.

FALCONER, D. S., and T. F. C. MACKAY, 1996 Introduction to Quantitative Genetics. 4th Edition. Longman, New York.

GEORGES, M., D. NIELSEN, M. MACKINNON, A. MISHRA, and R. OKIMOTO et al., 1995  Mapping quantitative trait loci controlling milk production by exploiting progeny testing. Genetics 139:907-920[Abstract].

GOMEZ-RAYA, L., D. I.VÅGE, I. OLSAKER, H. KLUNGLAND, G. KLEMETSDAL et al., 1996 Mapping QTL affecting traits of economical importance in Norwegian cattle. Proceedings of 47th Annual Meeting of the European Association for Animal Production, Lillehammer, Norway, August 25–29, 1996, p. 39.

HALEY, C. S. and S. A. KNOTT, 1992  A simple regression method for mapping quantitative trait loci in line crosses using flanking markers. Heredity 69:315-324[Medline].

HOLLANDER, M., and D. A. WOLFE, 1973 Nonparametric Statistical Methods. John Wiley & Sons, New York.

JANSEN, R. C., 1993  Interval mapping of multiple quantitative trait loci. Genetics 135:205-211[Abstract].

KAPPES, S. M., J. W. KEELE, R. T. STONE, R. A. MCGRAW, and T. S. SONSTEGARD et al., 1997  A second-generation linkage map of the bovine genome. Genome Res. 7:235-249[Abstract/Free Full Text].

KNOTT, S. A., J. M. ELSEN, and C. S. HALEY, 1996  Methods for multiple-marker mapping of quantitative trait loci in half-sib populations. Theor. Appl. Genet. 93:71-80.

KRUGLYAK, L. and E. S. LANDER, 1995a  A nonparametric approach for mapping quantitative trait loci. Genetics 139:1421-1428[Abstract].

KRUGLYAK, L. and E. S. LANDER, 1995b  Complete multipoint sib-pair analysis of qualitative and quantitative traits. Am. J. Hum. Genet. 57:439-454[Medline].

KÜHN, CH., R. WEIKARD, T. GOLDAMMER, S. GRUPE, and I. OLSAKER et al., 1996  Isolation and application of chromosome 6 specific microsattelite markers for detection of QTL for milk-production traits in cattle. J. Anim. Breed. Genet. 113:355-362.

LANDER, E. S. and P. GREEN, 1987  Construction of multilocus genetic linkage maps in humans. Proc. Natl. Acad. Sci. USA 84:2363-2367[Abstract/Free Full Text].

SPELMAN, R. L., W. COPPIETERS, L. KARIM, J. A. M. VAN ARENDONK, and H. BOVENHUIS, 1996  Quantitative trait loci analysis for five milk production traits on chromosome six in the dutch Holstein-Friesian population. Genetics 144:1799-1808[Abstract].

VAN RADEN, P. M. and G. R. WIGGANS, 1991  Derivation, calculation, and use of National Animal Model Information. J. Dairy Sci. 74:2737-2746[Abstract].

VILKKI, H. J., D.-J. DE KONING, K. ELO, R. VELMALA, and A. MAKI-TANILA, 1997  Multiple marker mapping of quantitative trait loci of Finnish dairy cattle by regression. J. Dairy Sci. 80:198-204[Abstract].

VOS, P., R. HOGERS, M. BLEEKER, M. REIJANS, and T. VAN DE LEE et al., 1995  AFLP: a new technique for DNA fingerprinting. Nucleic Acids Res. 23:4407-4414[Abstract/Free Full Text].

WEBER, J. L. and P. E. MAY, 1989  Abundant class of human DNA polymorphisms which can be typed using the polymerase chain reaction. Am. J. Hum. Genet. 44:388-396[Medline].

WELLER, J. I., Y. KASHI, and M. SOLLER, 1990  Power of daughter and granddaughter designs for determining linkage between marker loci and quantitatve trait loci in dairy cattle. J. Dairy Sci. 73:2525-2537[Abstract].

WILLIAMS, J. G., A. R. KUBELIK, K. J. LIVAK, J. A. RAFALSKI, and S. V. TINGEY, 1990  DNA polymorphisms amplified by arbitrary primers are useful as genetic markers. Nucleic Acids Res. 18:6531-6535[Abstract/Free Full Text].

ZENG, Z-B., 1993  Theoretical basis of separation of multiple linked gene effects on mapping quantitative trait loci. Proc. Natl. Acad. Sci. USA 90:10972-10976[Abstract/Free Full Text].