Molecular Dissection of a Quantitative Trait Locus: A Phenylalanine-to-Tyrosine Substitution in the Transmembrane Domain of the Bovine Growth Hormone Receptor Is Associated With a Major Effect on Milk Yield and Composition

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

Communicating editor: C. HALEY

	ABSTRACT

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

We herein report on our efforts to improve the mapping resolution of a QTL with major effect on milk yield and composition that was previously mapped to bovine chromosome 20. By using a denser chromosome 20 marker map and by exploiting linkage disequilibrium using two distinct approaches, we provide strong evidence that a chromosome segment including the gene coding for the growth hormone receptor accounts for at least part of the chromosome 20 QTL effect. By sequencing individuals with known QTL genotype, we identify an F to Y substitution in the transmembrane domain of the growth hormone receptor gene that is associated with a strong effect on milk yield and composition in the general population.

WITH the development of genome-wide marker maps for several species it is becoming possible to map quantitative trait loci (QTL) underlying the genetic variation for continuously distributed phenotypes of medical and agronomical importance (e.g., ANDERSSON 2001 ; FLINT and MOTT 2001 ; MACKAY 2001 ; MAURICIO 2001 ). The mapping resolution that is achieved with most experimental designs, however, is in the tens of centimorgans at best, therefore precluding efficient marker-assisted selection let alone positional cloning of the corresponding genes.

When working with model organisms or plants, strategies to improve the mapping resolution most often require breeding of a large number of progeny to increase the density of crossovers in the chromosome regions of interest (e.g., DARVASI 1998 ). When working with humans or farm animals, this approach is not practical and alternative strategies need to be identified. One approach that has recently received considerable attention is linkage disequilibrium (LD) mapping, which aims at exploiting historical recombinants. However, as useful LD is expected to extend only over limited distances (e.g., 60 kb in the human; REICH et al. 2001 ) this approach requires a commensurate increase in marker density.

In some livestock populations, including dairy cattle, LD has been shown to extend over very long chromosome segments when compared to human populations (FARNIR et al. 2000 ). This has raised hope that LD would be readily exploitable in these species using the presently available medium density maps (e.g., KAPPES et al. 1997 ). These hopes have been supported by a series of initial, promising results dealing with both simple (e.g., CHARLIER et al. 1996 ) and complex inherited traits (e.g., GRISART et al. 2002 ). Potential downsides associated with this long-range LD are (i) a possible limited mapping resolution and (ii) the occurrence of association in the absence of linkage due to gametic association between nonsyntenic loci.

We herein use LD to improve our molecular understanding of a QTL influencing milk yield and composition that was previously mapped to bovine chromosome 20 (GEORGES et al. 1995 ; ARRANZ et al. 1998 ).

	MATERIALS AND METHODS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Pedigree material:
The pedigree material used in this study was composed of the following:

• Data set I: a previously described Black-and-White Holstein-Friesian granddaughter design (GDD) sampled in the Netherlands and composed of 22 paternal half-sib families for a total of 987 bulls (SPELMAN et al. 1996 ; COPPIETERS et al. 1998A ).

• Data set II: 276 progeny-tested Holstein-Friesian sires sampled in the Netherlands.

• Data set III: 1550 progeny-tested Holstein-Friesian sires sampled in New Zealand.

• Data set IV: 959 progeny-tested Jersey sires sampled in New Zealand.

• Data set V: 485 Holstein-Friesian cows sampled in New Zealand.

• Data set VI: 387 Jersey cows sampled in New Zealand.

Phenotypes:
Phenotypes were, respectively, daughter yield deviations (DYD) for bulls, lactation values (LV, the unregressed first lactation yield deviations) for cows, average parental predicted transmitting abilities (PTA) for bulls and cows for milk, protein, and fat yield, as well as protein and fat percentage (VAN RADEN and WIGGANS 1991 ). DYDs, lactation values, and PTA were obtained directly from CR-DELTA (Netherlands; data sets I and II) or Livestock Improvement Corporation (LIC, New Zealand; data sets III–VI), respectively.

Map construction:
Microsatellite genotyping, map construction, and information content mapping were performed as previously described (COPPIETERS et al. 1998A ). Sequence information for the primers used for PCR amplification of anonymous type II microsatellite markers can be obtained from ArkDB (http://www.thearkdb.org/species.html). The following primers were designed on the basis of HEAP et al. 1995 to amplify a microsatellite in the promotor region of the growth hormone receptor gene: GHRJA.UP, 5'-TGCTCTAATCTTTTCTGGTACCAGG-3', and GHRJA.DN, 5'-TCCTCCCCAAATCAATTACATTTTCTC-3'.

Conventional QTL mapping:
QTL mapping was performed by multimarker regression (KNOTT et al. 1996 ) using the previously described HSQM software (COPPIETERS et al. 1998B ). Chromosome-wide significance thresholds were determined by permutation as previously described (CHURCHILL and DOERGE 1995 ; COPPIETERS et al. 1998B ). Segregating sire families were identified on the basis of the results of within-family analyses as previously described (COPPIETERS et al. 1998A ).

Haplotype-based test for association:
Assumptions: Assume a QTL that is characterized by two additively acting alleles, Q and q, that segregate in the population of interest with respective allelic frequencies of f_Q and (1 - f_Q). Assume that the Q allele appeared in the population by mutation or migration on a chromosome with haplotype "H" for a series of flanking markers. All other haplotypes are pooled and referred to as "O." At the present generation the H haplotype may still be in LD with the Q allele by an amount D. The H to O haplotype substitution effect can then be shown to equal

where a corresponds to one-half the difference between the phenotypic values of QQ vs. qq individuals, and f_H corresponds to the population frequency of the H haplotype (FALCONER and MACKAY 1996 ).

Test for association: Knowing that in our GDD, phased-marker genotypes are available for all sons and their sires but not their dams as these are not marker genotyped, and defining T_i as [DYD_i - PA_i], where DYD_i is the daughter yield deviation of son i and PA_i is the average predicted transmitting ability (VAN RADEN and WIGGANS 1991 ) of the sire and dam of son i, one can express the expected value of T_i as a function of the marker genotype of the sire's chromosomes (SC) and the marker genotypes of the paternal (PC) and maternal gametes (MC) inherited by son i, as shown in Table 1.

Expected values of T_i can be seen to be linear functions of the unknown haplotype substitution effect, . A least-squares estimator of can therefore easily be obtained by linear regression, while the ratio

which is distributed as an F-statistic with 1 and n - 2 d.f., can be used to measure the evidence in favor of a statistically significant haplotype substitution effect. In this n corresponds to the number of sons available in the GDD, SSE to the residual sum of squares, and SSR to the regression sum of squares.

By using T_i as phenotype, one essentially performs a transmission disequilibrium test (TDT; SPIELMAN et al. 1993 ), which simultaneously tests for association and linkage. As the dams are not genotyped, however, our TDT reduces in part to a conventional association test.

Choice of markers and haplotypes: So far, we have not defined which of the m markers available on the chromosome have to be considered when defining a haplotype. As the exact location of the QTL and the size of the haplotype that will maximize are both unknown, all possible windows comprising between one and m adjacent markers are tested separately. We thus examine m windows of one marker, (m - 1) windows of two markers, (m - 2) windows of three markers, ... , and one window of m markers.

Having selected the markers comprising the haplotype, one still has to choose the H haplotype among all haplotypes encountered in the population. In the proposed approach, the haplotypes that were successively considered as H haplotypes correspond to the chromosomes of the s sires in the GDD that were known to be heterozygous Qq for the QTL on the basis of the results of a marker-assisted segregation analysis performed on their sons (see above). As one does not a priori know which of the sire's homologs carries the Q allele, the haplotypes corresponding to both chromosomes are examined, for a total of 2s homologs.

When estimating the substitution effect of the haplotypes of a given sire, his sons were eliminated from the data set, to avoid extracting information that would be redundant with the linkage analysis.

Significance thresholds: The F-ratio defined above does not account for the multiple tests that are performed, i.e., the (m² + m)/2 marker windows tested for each of the 2s homologs. We accounted for multiple testing by applying a permutation test. The phenotypes and marker genotypes were shuffled 1000 times and the 2s(m² + m)/2 tests were performed on each permutated data set. The highest F-ratios obtained with the real data were then compared with the highest F-ratios obtained across the 1000 permutations.

Simultaneous mining of linkage and linkage disequilibrium:
QTL fine-mapping exploiting both linkage and LD: The utilized mapping method is implemented in the LD variance component mapping (LDVCM) programs and can be summarized as follows. Testing for the presence of a QTL at map position p of the studied chromosome was performed as follows:

1. For all markers on the studied chromosome, we determine the marker linkage phase of the sires and sons as described (FARNIR et al. 2002 ). As a consequence, the marker data then consist of 2s sire chromosomes (SC), n paternally inherited chromosomes of the sons (PC), and n maternally inherited chromosomes of the sons (MC), where s and n correspond, respectively, to the number of sire families and the number of sons in the GDD. From the genotypes of the PC, we can easily compute the probability that son i inherited the "left" (_p) or "right" (_p = 1 - _p) SC from its sire at map position p as described (COPPIETERS et al. 1998B ).

2. We compute identity-by-descent (IBD) probabilities (_p) for all pairwise combinations of SC and MC using the method described by MEUWISSEN and GODDARD 2001 . This method approximates the probability that two chromosomes are IBD at a given map position conditional on the identity-by-state (IBS) status of flanking markers, on the basis of coalescent theory (HUDSON 1985 ). Windows of 16 markers were considered to compute _p.

3. Using (1 - _p) as a distance measure, we apply the UPGMA hierarchical clustering algorithm (e.g., MOUNT 2001 ) to generate a rooted dendrogram representing the genetic relationship—at position p—between all SC and MC haplotypes encountered in the population.

4. We use the logical framework provided by this dendrogram to group the SC and MC in functionally distinct clusters. A cluster is defined as a group of haplotypes that coalesce into a common node. A useful feature of UPGMA trees in this regard is that the distance (1 - _p) between all the haplotypes that coalesce into a given node is less than or equal to two times the distance between the node and any of these haplotypes. As a consequence, the tree is scanned downward from the root and branches are cut until nodes are reached such that all coalescing haplotypes (i.e., all haplotypes within the cluster) have a distance measure (1 - _p) < T (KIM and GEORGES 2002 ).

5. We model the sons' phenotypes (DYDs) using the following linear model:
   

   

   y is the vector of phenotype records of all sons. b is a vector of fixed effects, which in this study reduces to the overall mean. X is the incidence matrix relating fixed effects to individual sons, which in this study reduces to a vector of ones. h is the vector of random QTL effects corresponding to the defined haplotype clusters. Z_h is an incidence matrix relating haplotype clusters to individual sons. In Z_h, a maximum of three elements per line can have nonzero value: "1" in the column corresponding to the cluster to which the MC haplotype belongs and _p and _p in the columns corresponding, respectively, to the haplotype clusters of the right and left SC. If the SC and/or MC belong to the same cluster, the corresponding coefficients are added. u is the vector of random individual polygenic effects ("animal model", LYNCH and WALSH 1997 ). Z_u is a diagonal incidence matrix relating individual polygenic effects to individual sons. e is the vector of individual error terms.

   Haplotype cluster effects with corresponding variance, ²_H, individual polygenic effects with corresponding variance, ²_A, and individual error terms with corresponding variance, ²_E, were estimated using average information restricted maximum likelihood (AIREML; JOHNSON and THOMPSON 1995 ), by maximizing the restricted log-likelihood function L:
   

   

   In this,
   

   

   Because we assume that the covariance between the QTL effects of the different haplotype clusters is zero, H reduces to an identity matrix. This differentiates our approach from that of MEUWISSEN and GODDARD 2000 , in which H is the matrix of between-haplotype IBD probabilities. A is the additive genetic relationship matrix (LYNCH and WALSH 1997 ).

6. Steps 4 and 5 are repeated for all possible values of T (from 0 to 1), to identify a restricted maximum likelihood (REML) solution for map position p. By analogy with FARNIR et al. 2002 , we denote the hypothesis corresponding to this REML solution as H₂.

QTL mapping exploiting linkage only: Note that the previous model can be extended with minor modifications to map QTL by exploiting linkage information only. This is simply achieved by ignoring all MCs and considering that all SCs belong to distinct haplotype clusters, irrespective of their marker genotype. REML solutions for the different parameters can be found as described in the previous section. Again by analogy with FARNIR et al. 2002 , we refer to the corresponding hypothesis as H₁.

Hypothesis testing and significance thresholds: The log likelihoods of the data under the H₂ and H₁ hypotheses are compared with that under the null hypothesis, H₀, of no QTL at map position p. The latter is computed as described above but using the reduced model,

Evidence in favor of a QTL at map position p can then be expressed as a lod score:

As customary when performing interval mapping, we are sliding the hypothetical position of the QTL throughout the chromosome map and compute lod scores at each map position as described to generate chromosome-wide lod score profiles.

KIM and GEORGES 2002 have shown by simulation that when analyzing a chromosome of 100 cM with a marker density of one marker every 5 cM, 2 x (L_H1/2 - L_H0) has (under the null hypothesis) an approximate chi-square distribution with 2 d.f. corrected (Bonferroni correction) for two and six independent tests when considering, respectively, H₁ and H₂. Chromosome-wide significance levels were computed from these distributions in this study.

Sequencing the coding portion of the growth hormone receptor from genomic DNA:
To develop primers that would allow us to conveniently amplify and sequence the entire growth hormone receptor (GHR) coding sequence from bovine genomic DNA, we screened a bovine bacterial artificial chromosome (BAC) library (WARREN et al. 2000 ) using standard procedures with an oligonucleotide probe complementary to exon 10 and isolated eight GHR-containing clones. DNA from one of these clones was used as template for sequencing the intron-exon boundaries using exonic primers designed on the basis of the bovine cDNA sequence (e.g., HAUSER et al. 1990 ) and predicted to flank exon-intron boundaries assuming conservation of intron position between human and cattle (e.g., GODOWSKI et al. 1989 ). On the basis of the obtained intronic information, we then designed primers (Table 2) to amplify and sequence most of the GHR coding sequence from genomic DNA using standard procedures. Sequence traces were analyzed with the POLYPHRED software (NICKERSON et al. 1997 ).

Oligonucleotide ligation assay:
An oligonucleotide ligation assay (OLA) test to genotype the GHR F279Y, Nt864-33(T-G), Nt933+21(A-G), Nt1095(T-C), N528T, and Nt1922(C-T) single-nucleotide polymorphisms (SNPs) in multiplex was developed as previously described (KARIM et al. 2000 ). The primers used for the PCR amplification step and the ligation reaction are reported in Table 3.

Estimating the effect on milk yield and composition associated with the F279Y polymorphism in the general dairy cattle population:
The effect of the F279Y genotype on milk yield and composition was estimated using the model

where y_i were DYDs when studying bulls or lactation values when studying cows, g_i is a fixed effect corresponding to the F279Y genotype (FF, FY, or YY), a_i is a random polygenic component accounting for all known pedigree relationships ["animal model" (LYNCH and WALSH 1997 ) including ungenotyped individuals whose phenotypes were ignored], and e_i is a random residual. Maximum-likelihood solutions for g_i, a_i, and e_i were obtained using the MTDFREML program (BOLDMAN 1993 ), setting ²_a/(²_a + ²_e) for yield (percentage) traits at 70% (75%) and 35% (50%) for DYDs and LVs, respectively.

The statistical significance of the F279Y genotype effect was estimated from

where SSM_F, SSM_R, and SSE_F are the sum of squares due to the full model, reduced model, and error (full model), respectively, which is distributed as an F-statistic with 3 and (n - 3) d.f.

The proportion of the trait variability due to the QTL genotype (r²_QTL) was approximated as

where f_FF, f_FY, and f_YY are the frequencies of the FF, FY, and YY genotypes; , , and are the genotype means estimated using the mixed model, SST is the total sum of squares, and is the overall trait mean estimated as

	RESULTS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Construction of a high-density microsatellite map of bovine chromosome 20:
To refine the map position of the chromosome 20 QTL, we first increased the marker density on this chromosome. We genotyped data set I for 22 additional, publicly available microsatellites known to map to bovine chromosome 20 as well as for a microsatellite in the promotor region of the bovine growth hormone receptor gene (GHRJA). A male linkage map composed of 29 markers covering 85 cM [Kosambi (K)] was constructed with average marker interval of 3 cM (K; Fig 1). The information content of the corresponding map was computed as previously described (COPPIETERS et al. 1998A ). It is >80% for most of the chromosome length. This is considerably larger than the average marker heterozygosity among the sires of data set I (45%), which illustrates the benefit of multipoint mapping methods. The map, shown in Fig 1, also reports the position of the prolactin receptor (PRLR) gene deduced from segregation data of prolactin receptor SNPs in the same pedigree material (S. MOISIO, unpublished observations). Note that in the human, the GHR gene is located in band 5p13.3 at map position 37.4 Mb on the "golden path" human sequence (Ensembl Human Genome Server: http://www.ensembl.org), while the PRLR gene is located in band 5p13.1 at map position 50.9 Mb, i.e., at 15 Mb from the former. The genetic distance separating the bovine GHR and PRLR genes is therefore compatible with the human data.

View larger version (48K): In this window In a new window Download PPT slide   Figure 1. Chromosome 20 microsatellite map, conventional QTL mapping, and haplotype-based test for association. Chromosome 20 microsatellite map: the name of the corresponding markers is given at the top and their respective position in centimorgans (Kosambi) at the bottom. GHRJA corresponds to a microsatellite marker in the promotor of the growth hormone receptor gene. The most likely position of the prolactin receptor (PRLR) inferred from the segregation of SNP markers (S. MOISIO, unpublished results) is given. Markers that could not be ordered with odds >1000 are in parentheses. The thin black curve corresponds to the information content (IC; expressed as a percentage—right y-axis) obtained in the GDD. The open diamonds correspond to the percentage of founder sires in the GDD that are heterozygous for the corresponding microsatellite marker. Conventional QTL mapping: the thick black and gray curves correspond to the location scores obtained, respectively, for milk protein percentage and milk fat percentage. Location scores are expressed as log(1/p) (left y-axis), where p corresponds to the chromosome-wide probability of obtaining the corresponding signal under the null hypothesis of no QTL determined by phenotype permutation. Most likely QTL positions obtained across 1000 bootstrap samples (left y-axis) are given as black vertical bars. The resulting 95% confidence interval is shown as a thick black horizontal bar at the top (C.I.). Haplotype-based test for association: marker windows showing significant effects in the haplotype-based association test are shown as gray cylinders. Their position with respect to the left y-axis corresponds approximately to their significance level determined as described in MATERIALS AND METHODS.

Conventional QTL mapping using a dense marker map:
We then used these novel microsatellite genotypes to repeat a QTL-mapping analysis in data set I. Fig 1 reports the location scores that were obtained by multimarker regression in the across-family analysis along the newly generated chromosome 20 marker map. As expected, these results confirm the presence of a QTL with strong effect on protein percentage at most likely position 49 cM. The QTL affects fat percentage to a lesser extent and has only very modest influence on the yield traits (data not shown).

Bootstrap analyses were performed for protein percentage according to VISSCHER et al. 1996 to estimate the 95% confidence interval (C.I.) for the position of the QTL. Fig 1 illustrates the distribution of the most likely position of the QTL across 1000 bootstrap samples as well as the deduced 95% C.I. It can be seen that the C.I. covers 50 cM, which in essence corresponds to the distal half of chromosome 20 and therefore to a very poor location of the QTL.

We then performed within-family regression analyses on protein percentage as described (ARRANZ et al. 1998 ) to identify sire families that were segregating for this QTL. Two such families were identified in data set I: families 1 and 18 (data not shown).

Refining the map position of a QTL—use of a haplotype-based test for association:
The previously described within-family analyses indicate that sires 1 and 18 are heterozygous for QTL alleles with large substitution effects (Q) on chromosome 20. Previous work within the same population revealed extensive genome-wide linkage disequilibrium due to random drift (FARNIR et al. 2000 ). We therefore hypothesized that the marker haplotypes flanking the Q alleles in the segregating sires might well be in linkage disequilibrium with the same Q alleles in the general population as well. To test this hypothesis, we measured the effect on protein percentage of the sire haplotypes in the general population using the haplotype-based test for association described in MATERIALS AND METHODS.

By doing so, we identified five haplotype windows that were yielding significant F-ratios (P < 0.01 after correction for multiple testing), corresponding to substitution effects of 0.03% milk protein. The corresponding haplotypes were all derived from a chromosome segment that was shared identical by state by sires 1 and 18. The sons of both sires were eliminated from the data set prior to performing the test for association. Fig 1 shows the position and statistical significance of the corresponding marker windows. It can be seen that their position centers around the TGLA153-GHRJ marker pair, corresponding to a minor peak for protein percentage, but the most likely QTL position when analyzing fat percentage. This result strongly suggests that a gene in the vicinity of these markers indeed contributes to the QTL effect observed on bovine chromosome 20.

Refining the map position of a QTL—combined linkage and LD analysis:
To confirm the findings obtained with the haplotype-based test for association, we analyzed data set I using the LDVCM program for combined linkage and LD mapping. Fig 2 shows the location scores that were obtained with this approach for protein percentage. The profile obtained when considering only linkage information essentially parallels that obtained by multimarker regression (cf. Fig 1), although the lod scores are slightly less significant (z _max = 1.8; chromosome-wide P value = 0.016). When including linkage disequilibrium information, however, a very significant lod score of 8.5 corresponding to a chomosome-wise P value of 1.5E-8 was obtained at map position 43 cM, i.e., very close to the chromosome region identified by the haplotype-based association test. Using the same approach, highly significant lod scores were obtained in the same chromosome region for fat percentage (position 43 cM, lod score 5.9, P value 7.5E-6), milk yield (position 43 cM, lod score 4.5, P value 0.00018), fat yield (position 46 cM, lod score 3.2, P value 0.0047), and protein yield (position 43 cM, lod score 5.2, P value 3.7E-5; data not shown). These results therefore supported the existence of a QTL influencing milk yield and composition in the vicinity of the GHR gene.

View larger version (41K): In this window In a new window Download PPT slide   Figure 2. Lod score profiles obtained for protein percentage along the chromosome 20 map using the LDVCM programs. The names of the markers comprising the map are given at the top and their respective positions in centimorgans (Kosambi) at the bottom. The gray curve is obtained by considering linkage information only, while all other curves are obtained by considering both linkage and LD. Black, basic chromosome 20 microsatellite marker map; red, chromosome 20 microsatellite marker map plus six GHR SNPs [F279Y, Nt864-33(T-G), Nt933+21(A-G), Nt1095(T-C), N528T, and Nt1922(C-T)]; light green, chromosome 20 microsatellite marker map plus five GHR SNPs (F279Y dropped); Blue, chromosome 20 microsatellite marker map plus four PRLR SNPs; dark green, chromosome 20 microsatellite map plus F279Y treated as a fixed effect. The yellow diamonds correspond to the lod scores obtained by single-point analysis with the individual GHR SNPs. The names of the corresponding SNPs are given in the adjacent boxes.

Scanning the bovine growth GHR gene for DNA sequence polymorphisms:
Knowing the key role played by the growth hormone axis in the initiation and maintenance of lactation, the GHR obviously stood out as a strong positional candidate for the corresponding QTL effect. Assuming that the GHR gene would indeed account for at least part of the QTL effect, we predicted on the basis of the haplotype-based test for association that sires 1 and 18 would both be heterozygous for a mutation causing the GHR to be functionally different. We therefore decided to scan the coding portion of the GHR gene for DNA sequence polymorphisms in these animals. Intronic primers allowing for the convenient amplification and sequencing of exons 3–10 of the GHR were developed as described in MATERIALS AND METHODS. Analysis of the sequence traces obtained from five Holstein-Friesian individuals including sires 1 and 18 revealed 10 SNPs in the GHR gene. Fig 3 reports the position and nature of the corresponding SNPs.

View larger version (28K): In this window In a new window Download PPT slide   Figure 3. Schematic representation of the bovine GHR gene. The 10 exons are shown as large cylinders and labeled by exon number. Coding sequences are shown in dark gray, sequenced 3' and 5' UTR sequences in gray, and unsequenced UTR regions in light gray. Introns are shown as interrupted thin cylinders. SNPs are marked as black lines connected with a box detailing the corresponding DNA sequences. The SNPs for which sires 1 and 18 were found to be hetereozygous are marked by stars.

Four of these are SNPs located in introns [Nt71-85(del1), Nt71-12(T-C), Nt864-33(T-G), and Nt933+21(A-G)], one is an SNP located in the 3' UTR of the GHR gene [Nt1922(C-T)], and three are synonymous mutations in third codon positions [Nt1095(C-T), Nt1635(C-T), and Nt1809(C-T)]. None of these are a priori likely to affect the function of the GHR gene.

The two remaining SNPs, however, modify the amino acid sequence of the GHR receptor. A T to A substitution in exon VIII results in the nonconservative replacement of a neutral phenylalanine with an uncharged but polar tyrosine residue (F279Y). The corresponding phenylalanine residue is located within the transmembrane domain of the GHR and is conserved among all analyzed mammals (human, baboon, rabbit, mouse, rat, dog, pig, sheep, and opossum) except guinea pig, where it is nevertheless replaced by a neutral leucine residue. In chicken and pigeon, the corresponding residue is also a neutral isoleucine.

An A to C substitution in exon X results in the replacement of an asparagine with a threonine (N528T), both amino acids being polar uncharged residues. This residue is less conserved during evolution, being either an asparagine (human, rabbit, pig, and chicken) or a serine residue (ovine, mouse, and rat). It is, however, only one residue removed from one of the highly conserved, potentially phosphorylated intracellular tyrosine residues.

Sires 1 and 18, which were both heterozygous for the GHR-containing marker haplotype associated with a highly significant substitution effect on protein percentage in the association test, were heterozygous for SNPs Nt71-85(del1), Nt864-33(T-G), Nt933+21(A-G), and most importantly F279Y. Given the effect of this SNP on the primary sequence of the GHR and therefore possibly on its function, F279Y stood out as a prime candidate for the mutation causing the observed QTL effect.

Inclusion of SNPs in the combined linkage and LD analysis dramatically increases the lod score at the GHR locus:
We constructed an OLA as described (KARIM et al. 2000 ) for multiplex genotyping of the F279Y, Nt864-33(T-G), Nt933+21(A-G), Nt1095(T-C), N528T, and Nt1922 (C-T) SNPs and applied it to data set I. The linkage phase was determined as described (FARNIR et al. 2002 ). Fig 4 shows the frequency distribution of the GHR haplotypes as measured in the MCs (see above). It shows that at least 13 distinct haplotypes occur in the Dutch Holstein-Friesian population; however, 3 of these account for 85% of the chromosomes in this population.

View larger version (20K): In this window In a new window Download PPT slide   Figure 4. Frequency distribution of the GHR SNP haplotypes in the Dutch Holstein-Friesian population.

We placed the GHR SNP haplotype by linkage analysis on the chromosome 20 marker map at position 42.7 cM, coinciding with the GHRJ microsatellite as expected.

A combined linkage and LD analysis was then performed using the LDVCM software, including the new GHR SNP genotypes. As shown in Fig 2 for protein percentage, inclusion of the GHR SNPs increased the maximum lod score by 3.8 units, yielding a maximum lod score of 12.3 at position 43 cM, i.e., just distal of the GHR gene. Table 4 reports the corresponding variance component estimates.

Including the GHR SNPs in the LDVCM analysis had a comparable effect when analyzing fat percentage. The lod score increased from 5.9 to 7.8, maximizing exactly at the GHR gene (Table 4). The effect was more modest for milk yield and fat yield, increasing the lod scores by, respectively, 0.4 and 0.1 units but maximizing in both instances on the GHR gene (Table 4). For protein yield, however, inclusion of the GHR SNPs resulted in a marked decrease of the lod scores, dropping from 5.2 to 1.7 or less in the region of the GHR gene (Table 4).

For comparison, performing a combined linkage and LD analysis after inclusion of a haplotype composed of four PRLR SNPs resulted in a local decrease in the lod score values for all traits (see Fig 2 for protein percentage).

These results strongly suggested that either the GHR itself or a gene in strong LD with GHR underlies the observed QTL effect.

Unique status of the F279Y polymorphism with regard to the chromosome 20 QTL effect:
We then performed three tests to determine the relative contribution of the different SNPs to the increase in signal noted for protein percentage. First, we reran LDVCM analyses by sequentially dropping one of the six GHR SNPs comprising the GHR SNP haplotype. While dropping the Nt864-33(T-G), Nt933+21(A-G), Nt1095(T-C), N528T, and Nt1922 (C-T) SNPs did not significantly alter the lod score profiles (data not shown), dropping the F279Y SNP virtually annihilated the entire gain obtained by considering the complete GHR SNP haplotype (Fig 2).

Second, we ran LDVCM analyses considering only the microstallite markers to compute pairwise IBD probabilities (_p) between SCs and MCs, but sequentially including the GHR and PRLR SNP genotypes as fixed effects in our mixed model. Inclusion of the F279Y SNP as a fixed effect completely erased the lod score peak in the GHR region (Fig 2). With none of the other SNPs, whether in the GHR or PRLR genes, did the lod scores drop below 5.7 (data not shown). Note that including the PRLR SNPs also did not markedly affect the minor lod scores peak around position 65 cM, i.e., the approximate location of the PRL gene.

Finally, we used LDVCM to estimate the effects of the different GHR SNPs individually (i.e., without considering flanking marker data): F279Y would yield a lod score of 11.5 while the other SNPs would respectively yield lod scores of only 0.75 [Nt864-33(T-G)], 0.22 [Nt933+21 (A-G)], 0 [Nt1095(T-C)], 2.18 (N528T), and 1.77 [Nt1922 (C-T)] (Fig 2).

Altogether, these results clearly point toward a unique status of the F279Y polymorphism with regard to the chromosome 20 QTL effect, indicating that this SNP is either directly responsible for the QTL effect or tightly associated with the actual causal mutation.

Fig 5 shows the segregation of the F and Y alleles within the haplotype clusters maximizing the LDVCM lod scores when analyzing protein percentage including all six GHR SNPs. The REML solution is associated with a grouping in 22 haplotype clusters of which 17 are perfectly homogeneous with regard to the F279Y polymorphism.

View larger version (29K): In this window In a new window Download PPT slide   Figure 5. UPGMA dendrogram representing the genetic relationship between the SC and MC haplotypes at position 43.4 cM (interval GHR-TGLA53), yielding the most likely REML solution when analyzing protein percentage. The dotted vertical line corresponds to the optimal threshold value, T. The arrows point toward the nodes in which all haplotypes belonging to a given haplotype cluster coalesce. The arrows are solid when the corresponding cluster contains either only the F or only the Y allele for the F279Y SNP. They are open when the cluster contains both alleles. The bars at the right of the tree indicate the F279Y genotype of the corresponding haplotype: black equals Y and gray equals F. The bar is broken to mark the separation between clusters.

Effect of the F279Y GHR polymorphism on milk yield and composition in the general dairy cattle population:
To more accurately estimate the effect of the F279 GHR polymorphism on milk yield and composition, we genotyped data sets II–VI—corresponding to an additional 2772 bulls and 872 cows—for this SNP. Effects of the F279Y genotype on DYDs and LVs for milk yield (kilograms), protein yield (kilograms), fat yield (kilograms), protein percentage, and fat percentage were estimated using a mixed model including a fixed genotype effect and a random animal model to account for the polygenic background.

It can be seen from Table 5 that the F279Y polymorphism behaved in a very similar fashion in all analyzed populations, whether Dutch or New Zealander, Holstein-Friesian or Jersey. As expected, the effect of the F279Y genotype was—in all five data sets—most pronounced on protein percentage, accounting for 5.3–9.4% of the trait variance. The effect of the F279Y mutation was also clearly detectable in all these populations on fat percentage and to a lesser extent on milk yield. It accounts for between 2.3 and 5.5% of the variance in fat percentage and between 0.7% and 2.9% of the variance in milk yield. For milk yield, inheriting one Y allele increases the DYD for milk yield by an estimated ±67 kg to ±112 kg and the LV for milk yield by ±87 kg to ±162 kg. Effects of the F279Y genotype on fat and protein yield were in essence nonsignificant although a tendency toward a decrease in fat yield of 1.5–2.5 kg for every dose of Y allele was notable.

The fact that the F279Y mutation showed very comparable effects in all five analyzed populations strongly supports their bona fide nature and the causality of either F279Y itself or a tightly linked mutation.

The effects on LVs were expected to be approximately twice as large as their equivalent on DYDs (as LVs and DYDs are by definition measured on a different scale; FALCONER and MACKAY 1996 ), but to explain a smaller proportion of the trait variability as the heritability of LVs is lower than that of DYDs. For reasons that are as yet not explained, this tendency was not systematically apparent from our data (Table 5). Note also that the absolute values of the effects on the percentage traits cannot be directly compared between data sets I and II (Netherlands) vs. data sets III–VI (New Zealand) as the percentage traits are computed from the yield traits using different formulas in both countries.

	DISCUSSION

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

We herein illustrate how LD can be exploited to dissect and refine the map position of QTL in outbred dairy cattle populations. In this work, this led to the demonstration that a chromosomal region harboring the GHR gene contributes to the QTL effect on milk yield and composition that was previously identified on bovine chromosome 20 (e.g., GEORGES et al. 1995 ; ARRANZ et al. 1998 ). No evidence was found, on the contrary, for an effect of the chromosomal region harboring the PRLR gene, an equally intriguing positional candidate in the initial stages of this analysis (ARRANZ et al. 1998 ).

Two distinct approaches aimed at capturing LD information were utilized in this study and led to similar conclusions: a haplotype-based association test and a variant component-mapping approach. The latter differs from the method that was recently described by MEUWISSEN and GODDARD 2000 , MEUWISSEN and GODDARD 2001 , by including a hierarchical haplotype-clustering step aimed at solving numerical issues and increasing the computation speed. This clustering step has the additional advantage that it allows for the identification of haplotypes that are likely to be functionally different, which could then become the focus for further molecular work. This proposition is presently being tested in the analysis of other QTL segregating in dairy cattle populations.

By sequencing a very small number of individuals that were carefully selected on the basis of their inferred QTL genotype, we subsequently identified a nonconservative substitution of a highly conserved F residue in the transmembrane domain of the GHR. We present a series of systematic tests that strongly suggest that the corresponding F279Y polymorphism may be the direct cause (quantitative trait nucleotide; MACKAY 2001 ) of the consistently associated effects on milk yield and composition. Preliminary evidence (R. SPELMAN, unpublished observations) suggests that the F279Y polymorphism affects live weight, which would also be compatible with a direct effect of the GHR.

However, these observations do not formally prove the causality of the GHR and F279Y mutation. The fact that the lod score is highest in the multipoint analysis and maximizes just distally from the GHR gene might indicate that the causal mutation is different, affecting either the GHR or another gene. Examination of the genes located in the vicinity of the GHR on the human genomic sequence (http://www.ensembl.org/homo_sapiens/), however, does not reveal any obvious alternative candidate. It is also noteworthy that administering growth hormone to lactating cows has a major effect on protein yield (BAUMAN et al. 1999 ), while this trait is in essence not affected by this QTL.

The effects of the F279Y genotype on the indexes that are used as the basis for selection in the Netherlands and New Zealand [net economic merit index (INET) and breeding worth (BW), respectively] are highly significant (R. SPELMAN, unpublished observations). As a matter of fact, a retrospective survey of the genotypes of the New Zealand sires clearly indicates that the frequency of the F allele has increased in recent years and that the FF genotype increases the likelihood for a sire to be selected for breeding (R. SPELMAN, unpublished observations). As a consequence, we anticipate that this marker has the potential to be very useful for marker-assisted selection and to more effectively increase the frequency of the favorable F allele.

Data sets V and VI (composed of cows) allowed for the analysis of potential dominance effects between the F and Y alleles. Modest evidence in favor of dominance of the Y over the F allele was found for protein percentage (P < 0.05; data not shown). However, as the number of YY individuals was small, the power to detect significant dominance interactions was very limited. Preliminary analyses in these data sets also suggest that the F279Y mutation and the previously described K232A mutation in the bovine DGAT gene (GRISART et al. 2002 ) act in an additive manner.

We believe it unlikely that the F279Y or tightly associated polymorphism accounts for the entire chromosome 20 QTL effect. Indeed, examination of the location scores (e.g., Fig 1), suggests that additional more distally located genes might contribute to the QTL effect on BTA20 as well. We have also identified two sires that would clearly be heterozygous for a QTL on BTA20 despite being homozygous for the F279Y polymorphism (data not shown). For one of these (sampled in the Netherlands), the QTL would primarily affect fat percentage and milk yield (rather than protein percentage) and the most likely QTL location would be at map position 48 cM. This could indicate either that the F279Y is not the causative mutation or that another closely linked QTL is segregating in this population as well. For the other sire (sampled in New Zealand), the QTL would primarily affect protein and fat percentage but the most likely QTL position would be at map position 85 cM. For the latter sire, the GHR gene falls outside of the 95% confidence interval for the QTL. Having identified a first polymorphism with major effect on milk yield and composition on BTA20 should facilitate the identification of these other putative QTL.

	ACKNOWLEDGMENTS

This work was funded by grants from Holland Genetics, Livestock Improvement Corporation, Via Lactia Biosciences, the Vlaamse Rundvee Vereniging, the Ministère des Classes Moyennes et de l'Agriculture (Belgium), and European Union grants B104-CT95-0073 and PL970471. Sarah Blott was supported by a Wellcome International research fellowship.

Manuscript received July 1, 2002; Accepted for publication October 9, 2002.

	LITERATURE CITED

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

ANDERSSON, L., 2001  Genetic dissection of phenotypic diversity in farm animals. Nat. Rev. Genet. 2:130-138.[Medline]

ARRANZ, J.-J., W. COPPIETERS, P. BERZI, N. CAMBISANO, and B. GRISART et al., 1998  A QTL affecting milk yield and composition maps to bovine chromosome 20: a confirmation. Anim. Genet. 29:107-115.[Medline]

BAUMAN, D. E., R. W. EVERETT, W. H. WEILAND, and R. J. COLLIER, 1999  Production responses to bovine somatotropin in Northeast dairy herds. J. Dairy Sci. 82:2564-2573.[Abstract]

BOLDMAN, K. G., 1993 A Manual for Use of MTDFREML: A Set of Programs to Obtain Estimates of Variances and Covariances. U.S. Department of Agriculture, Agriculture Research Service.

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.

COPPIETERS, W., J. RIQUET, J.-J. ARRANZ, P. BERZI, and N. CAMBISANO et al., 1998a  A QTL with major effect on milk yield and composition maps to bovine chromosome 14. Mamm. Genome 9:540-544.[Medline]

COPPIETERS, W., A. KVASZ, J.-J. ARRANZ, B. GRISART, and F. FARNIR et al., 1998b  A rank-based nonparametric method to map QTL in outbred half-sib pedigrees: application to milk production in a granddaughter design. Genetics 149:1547-1555.[Abstract/Free Full Text]

DARVASI, A., 1998  Experimental strategies for the genetic dissection of complex traits in animal models. Nat. Genet. 18:19-24.[Medline]

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

FARNIR, F., W. COPPIETERS, J.-J. ARRANZ, P. BERZI, and N. CAMBISANO et al., 2000  Extensive genome-wide linkage disequilibrium in cattle. Genome Res. 10:220-227.[Abstract/Free Full Text]

FARNIR, F., B. GRISART, W. COPPIETERS, J. RIQUET, and P. BERZI et al., 2002  Simultaneous mining of linkage and linkage disequilibrium to fine map quantitative trait loci in outbred half-sib pedigrees: revisiting the location of a quantitative trait locus with major effect on milk production on bovine chromosome 14. Genetics 161:275-287.[Abstract/Free Full Text]

FLINT, J. and R. MOTT, 2001  Finding the molecular basis of quantitative traits: successes and pitfalls. Nat. Rev. Genet. 2:437-445.[Medline]

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]

GODOWSKI, P. J., D. W. LEUNG, L. R. MEACHAM, J. P. GALGANI, and R. HELLMISS et al., 1989  Characterization of the human growth hormone receptor gene and demonstration of a partial gene deletion in two patients with Laron-type dwarfism. Proc. Natl. Acad. Sci. USA 86:8083-8087.[Abstract/Free Full Text]

GRISART, B., W. COPPIETERS, F. FARNIR, L. KARIM, and C. FORD et al., 2002  Positional candidate cloning of a QTL in dairy cattle: identification of a missense mutation in the bovine DGAT gene with major effect on milk yield and composition. Genome Res. 12:222-231.[Abstract/Free Full Text]

HAUSER, S. D., M. F. MCGRATH, R. J. COLLIER, and R. J. KRIVI, 1990  Cloning and in vivo expression of bovine growth hormone receptor mRNA. Mol. Cell. Endocrinol. 72:187-200.[Medline]

HEAP, D., M. C. LUCY, R. J. COLLIER, C. K. BOYD, and W. C. WARREN, 1995  Nucleotide sequence of the promoter and first exon of the somatotropin receptor gene in cattle. J. Anim. Sci. 73:1529.[Medline]

HUDSON, R. R., 1985  The sampling distribution of linkage disequilibrium under an infinite alleles model without selection. Genetics 109:611-631.[Abstract/Free Full Text]

JOHNSON, D. L. and R. THOMPSON, 1995  Restricted maximum likelihood estimation of variance components for univariate animal models using sparse matrix techniques and average information. J. Dairy Sci. 78:449-456.[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]

KARIM, L., W. COPPIETERS, L. GROBET, A. VALENTINI, and M. GEORGES, 2000  Convenient genotyping of six myostatin mutations causing double-muscling in cattle using a multiplex oligonucleotide ligation assay. Anim. Genet. 31:396-399.[Medline]

KIM, J. J. and M. GEORGES, 2002  Evaluation of a new fine-mapping method exploiting linkage disequilibrium: a case study analysing a QTL with major effect on milk composition on bovine chromosome 14.. Asian-Aust. J. Anim. Sci. 15:1250-1256.

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

LYNCH, M., and B. WALSH, 1997 Genetics and Analysis of Quantitative Traits. Sinauer Associates, Sunderland, MA.

MACKAY, T. F. C., 2001  Quantitative trait loci in Drosophila. Nat. Rev. Genet. 2:11-20.[Medline]

MAURICIO, R., 2001  Mapping quantitative trait loci in plants: uses and caveats for evolutionary biology. Nat. Rev. Genet. 2:370-381.[Medline]

MEUWISSEN, T. H. and M. E. GODDARD, 2000  Fine mapping of quantitative trait loci using linkage disequilibria with closely linked marker loci. Genetics 155:421-430.[Abstract/Free Full Text]

MEUWISSEN, T. H. and M. E. GODDARD, 2001  Prediction of identity by descent probabilities from marker-haplotypes. Genet. Sel. Evol. 33:605-634.[Medline]

MOUNT, D. W., 2001 Bioinformatics: Sequence and Genome Analysis. Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY.

NICKERSON, D. A., V. O. TOBE, and S. A. TAYLOR, 1997  PolyPhred: automating the detection and genotyping of single nucleotide substitutions using fluorescence-based resequencing.. Nucleic Acids Res. 25:2745-2751.[Abstract/Free Full Text]

REICH, D. E., M. CARGILL, S. BOLK, J. IRELAND, and P. C. SABETI et al., 2001  Linkage disequilibrium in the human genome. Nature 411:199-204.[Medline]

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]

SPIELMAN, R. S., R. E. MCGINNIS, and W. J. EWENS, 1993  Transmission test for linkage disequilibrium: the insulin gene region and insulin-dependent diabetes mellitus (IDDM). Am. J. Hum. Genet. 52:506-516.[Medline]

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]

VISSCHER, P. M., R. THOMPSON, and C. S. HALEY, 1996  Confidence intervals in QTL mapping by bootstrapping. Genetics 143:1013-1020.[Abstract]

WARREN, W., T. P. SMITH, C. E. REXROAD, III, S. C. FAHRENKRUG, and T. ALLISON et al., 2000  Construction and characterization of a new bovine bacterial artificial chromosome library with 10 genome-equivalent coverage. Mamm. Genome 11:662-663.[Medline]

This article has been cited by other articles:

				C. Sandor and M. Georges On the Detection of Imprinted Quantitative Trait Loci in Line Crosses: Effect of Linkage Disequilibrium Genetics, October 1, 2008; 180(2): 1167 - 1175. [Abstract] [Full Text] [PDF]

				G. Banos, J. A. Woolliams, B. W. Woodward, A. B. Forbes, and M. P. Coffey Impact of Single Nucleotide Polymorphisms in Leptin, Leptin Receptor, Growth Hormone Receptor, and Diacylglycerol Acyltransferase (DGAT1) Gene Loci on Milk Production, Feed, and Body Energy Traits of UK Dairy Cows J Dairy Sci, August 1, 2008; 91(8): 3190 - 3200. [Abstract] [Full Text] [PDF]

				C. Zhu, M. Gore, E. S. Buckler, and J. Yu Status and Prospects of Association Mapping in Plants The Plant Genome, July 1, 2008; 1(1): 5 - 20. [Abstract] [Full Text] [PDF]

				X. Wang, C. Maltecca, R. Tal-Stein, E. Lipkin, and H. Khatib Association of Bovine Fibroblast Growth Factor 2 (FGF2) Gene with Milk Fat and Productive Life: An Example of the Ability of the Candidate Pathway Strategy to Identify Quantitative Trait Genes J Dairy Sci, June 1, 2008; 91(6): 2475 - 2480. [Abstract] [Full Text] [PDF]

				F. Guillaume, S. Fritz, D. Boichard, and T. Druet Short Communication: Correlations of Marker-Assisted Breeding Values with Progeny-Test Breeding Values for Eight Hundred Ninety-Nine French Holstein Bulls J Dairy Sci, June 1, 2008; 91(6): 2520 - 2522. [Abstract] [Full Text] [PDF]