Accuracy of direct genomic values in Holstein bulls and cows using subsets of SNP markers
© Moser et al. 2010
Received: 26 May 2010
Accepted: 16 October 2010
Published: 16 October 2010
Skip to main content
© Moser et al. 2010
Received: 26 May 2010
Accepted: 16 October 2010
Published: 16 October 2010
At the current price, the use of high-density single nucleotide polymorphisms (SNP) genotyping assays in genomic selection of dairy cattle is limited to applications involving elite sires and dams. The objective of this study was to evaluate the use of low-density assays to predict direct genomic value (DGV) on five milk production traits, an overall conformation trait, a survival index, and two profit index traits (APR, ASI).
Dense SNP genotypes were available for 42,576 SNP for 2,114 Holstein bulls and 510 cows. A subset of 1,847 bulls born between 1955 and 2004 was used as a training set to fit models with various sets of pre-selected SNP. A group of 297 bulls born between 2001 and 2004 and all cows born between 1992 and 2004 were used to evaluate the accuracy of DGV prediction. Ridge regression (RR) and partial least squares regression (PLSR) were used to derive prediction equations and to rank SNP based on the absolute value of the regression coefficients. Four alternative strategies were applied to select subset of SNP, namely: subsets of the highest ranked SNP for each individual trait, or a single subset of evenly spaced SNP, where SNP were selected based on their rank for ASI, APR or minor allele frequency within intervals of approximately equal length.
RR and PLSR performed very similarly to predict DGV, with PLSR performing better for low-density assays and RR for higher-density SNP sets. When using all SNP, DGV predictions for production traits, which have a higher heritability, were more accurate (0.52-0.64) than for survival (0.19-0.20), which has a low heritability. The gain in accuracy using subsets that included the highest ranked SNP for each trait was marginal (5-6%) over a common set of evenly spaced SNP when at least 3,000 SNP were used. Subsets containing 3,000 SNP provided more than 90% of the accuracy that could be achieved with a high-density assay for cows, and 80% of the high-density assay for young bulls.
Accurate genomic evaluation of the broader bull and cow population can be achieved with a single genotyping assays containing ~ 3,000 to 5,000 evenly spaced SNP.
In genomic selection (GS), selection decisions are made on genomic breeding values predicted from high-density single nucleotide polymorphic (SNP) markers. In dairy cattle, GS has the potential to double the rate of genetic gain to that of traditional breeding schemes due to a substantial reduction in generation intervals and increased selection intensities [1, 2]. Significant additional gains in GS schemes could be made if cows to breed sires and cows to breed cows were selected on genomic breeding values . Another benefit of genotyping cows may be lower rates of inbreeding: according to Daetwyler et al. , the use of GS can be expected to decrease the rate of inbreeding relative to conventional selection using BLUP breeding values, this effect will be greatest when larger numbers of both cows and potential sires are genotyped .
At the current price, high-density SNP genotyping assays are limited to applications involving elite sires and dams. An alternative is to use a more cost-effective low-density assay for genotyping more animals from the population. As shown for a single trait by Weigel et al. , a low-density assay comprising selected SNP can deliver a substantial portion of the gain of a high-density assay, possibly for a fraction of the price. However, the use of such a low-density array may still be limited if multiple traits require so many SNP that their genotyping cost is similar to the cost of a high-density chip.
The utility of low-density arrays will depend in part on the genetic architecture of the target trait. In GS, prediction equations are derived from a training set, where animals are phenotyped and genotyped to predict breeding values based only on the genotype information of evaluation animals. This requires that the markers are in sufficient LD with the QTL and simulation studies have shown that accuracy of genomic predictions increases as LD increases [6–10]. In the ideal case where every QTL is in perfect LD with a single marker and where a limited number of QTL with large effects account for the genetic variation, the maximum accuracy could be obtained with very few markers. However, there is increasing evidence that most complex traits are affected by very many QTL with a small effect (e.g. height in humans, [11–14]). This would imply that the training population would need to be genotyped with a high-density SNP panel in order to capture the effects of all QTL. Selecting individual SNP from high-density genotype data is complicated because the multicollinearity between SNP, i.e. two or more SNP in high but not complete LD, makes it difficult to identify 'important' SNP, as each SNP masks a part of the effect of other SNP and a single marker might be in LD with several QTL.
Utility of SNP subsets will also be affected by the relationship of the selection candidates to the training set. Although genomic predictions rely on LD between SNP and QTL, this LD can operate or be interpreted at a number of levels. In addition to population level LD, simulation studies and empirical data have demonstrated that the accuracy of prediction depends on the relatedness between animals in the training and evaluation populations [10, 15, 16]. At the extreme, even in the absence of LD between markers and QTL, markers can predict family relationships between animals. If animals in the training and evaluation data share DNA segments from a small number of ancestors, relatively few markers are required to trace the segments shared between related animals separated by only a few generations. A low-density assay of evenly spaced SNP might then provide sufficient accuracies of prediction of evaluation animals, as long as the information content of the subset of SNP is sufficient to estimate effects of distinct haplotypes.
The objective of this study was to evaluate the use of low-density SNP genotyping assays to predict the direct genomic value (DGV) of bulls and cows for commercially important traits in Holstein-Friesian dairy cattle. The impact of two analysis methods, the number of SNP needed for accurate DGV prediction, as well as strategies for SNP selection were explored.
Phenotype and genotype data were available on 2,144 Holstein-Friesian bulls and 510 Holstein-Friesian cows. The traits analysed included milk production traits (milk yield, fat yield, protein yield, fat percentage and protein percentage), an overall confirmation trait (overall type), survival index, Australian Profit Ranking (APR) and Australian Selection index (ASI). The ASI is an index given by (3.8 × protein ABV) + (0.9 × fat ABV) - (0.048 × milk ABV), APR is given by (3.8 × protein ABV) + (0.9 × fat ABV) - (0.048 × milk ABV) + (1.2 × milking speed ABV) + (2.0 × temperament ABV) + (3.9 × survival ABV) + (0.34 × cell count ABV) - (0.26 × live weight ABV) + (3.0 × daughter fertility), whereas survival is given by (0.5 × likeability) + (1.8 × overall type) + (3.0 × udder depth) + (2.2 × pin set).
Phenotype information was provided by the Australian Dairy Herd Improvement Scheme (ADHIS, http://www.adhis.com.au). The phenotypes used were deregressed breeding values (DRBV) for protein percentage, fat percentage, ASI, APR and survival, and daughter trait deviations (DTD) for protein yield, fat yield, milk yield and overall type. The deregression procedure removed the contribution of relatives other than daughters to the breeding values, as detailed in . For cows, trait deviations (TD) were available for protein yield, fat yield, milk yield and overall type, but no DRBV information was available for the other traits.
SNP genotypes were derived from the Illumina BovineSNP50 BeadChip (Illumina Inc., San Diego, USA). After quality control  and omitting SNP located on the sex chromosomes a total of 42,576 markers remained for the analysis.
Number of animals in training and validation sets and median number of records contributing to the phenotype per animal
Protein, Fat, Milk
Protein%, Fat%, ASI
where DGV is the vector of direct genomic values estimated with the marker genotypes, X is an incidence matrix that relates genotypes to individuals, and is the vector of SNP effects which is estimated by either one of the two methods described below.
where N is the number of training animals, y is a vector of phenotypes, is an unknown constant, X is a (N × p) matrix of genotypes encoded as 0 (homozygote), 1 (heterozygote) or 2 (other homozygote), is a vector of SNP effects, and I is a p × p identity matrix. The penalty term λ, which is the same for all SNP, overcomes the problem of ill-conditioning when multicollinearity among columns in X causes X'X to be singular, or nearly so. The system of equations was solved iteratively by the preconditioned conjugate gradient method . The 10-fold cross-validation procedure described in Moser et al. , with golden segment search , was used to locate the optimal λ within a given range. RR is equivalent to the BLUP method of Meuwissen et al.  and Habier et al. , which assumes that regression coefficients are independent random draws from a common normal distribution. Under the BLUP model, λ = σ 2 e /σ 2 g , where σ 2 e is the residual variance and σ 2 g the genetic variance.
In RR, the contribution of each bull can be weighted according to the number of daughters contributing to the phenotype. However, reliabilities of the phenotypes expressed as 'equivalent daughter contributions' were uniformly high, with small differences between the majority of training bulls, and weighting the contributions of bulls had no impact on the accuracy of DGV for method RR (results not shown).
The main idea of PLSR is to build orthogonal components (called 'latent components') from the original genotype matrix X. A PLSR component t = X w is a linear combination of the SNP that have maximal covariance with the response vector, under the additional assumption that components are mutually orthogonal . Subsequently, y is regressed on the linear combinations of markers.
Different algorithms to extract the latent components and to obtain regression coefficients exist. We implemented PLSR using an algorithm described in . The optimal model complexity (i.e. number of latent components), was estimated by ten fold cross-validation . Note that the PLSR regression coefficients differ from the ordinary least squares regression coefficients and the RR regression coefficients. The magnitude of the PLSR regression coefficients can be used to determine the relative influence of each SNP on the model and to select relevant SNP .
The absolute magnitude of the regression coefficients was used to determine which SNP are most influential in the training data set. To select subsets of markers, all 42,576 SNP were ranked by their absolute value of . The ranking of SNP was derived using a backward elimination procedure. The process started with a model including the complete set of 42,576 SNP. Subsequently in each step, a fraction of SNP with the smallest absolute value of the regression coefficients was dropped from the SNP list and the regression coefficients were recomputed. This re-computation is important as the regression coefficient of an individual SNP can strongly depend on other SNP that are in LD with the SNP of interest. The optimal model complexity (i.e. number of latent components) for PLSR and the value of λ for RR was estimated at each step by cross-validation.
In detail, we first fitted models including all 42,576 SNP. In the first iteration 40,000 SNP with the highest absolute value of the regression coefficient were retained in the SNP list. The number of SNP subsequently dropped in each iteration was 2,000 for subsets of up to 10,000 SNP, 500 SNP for subsets of up to 1,000 SNP, 100 SNP for subsets of up to 300 SNP and 20 SNP for subsets of up to 100 SNP.
Four alternative strategies of SNP subset selection were compared. Under strategy 1, separate subsets including the highest ranked SNP for each individual trait were created. Strategies 2-4 used a single subset of evenly spaced SNP. To select a subset of n evenly spaced SNP, we divided the total length of the autosomes into n intervals flanked by two markers to give segments of approximately equal length. Chromosome lengths and SNP positions were based on the physical map of cattle genome assembly Btau 4.0. Subsequently, the highest ranked SNP for ASI (strategy 2), APR (strategy 3) or the SNP with the highest minor allele frequency (MAF, strategy 4) in each segment, was added to the subset. Using the same subset of SNP, a model was then fitted for each trait to derive the prediction equations. Subsets of evenly spaced SNP were generated for sets including between 100 and 5,000 SNP. The accuracy of DGV obtained using a subset of SNP was compared to the accuracy from the analysis of all 42,576 SNP.
Maximum accuracy of DGV of cows and bulls derived by RR and PLSR
The panels in Figures 2 are ordered from high to low heritable traits (left-right, top-bottom) based on reported heritability estimates [27, 28]. Heritability of APR and ASI was assumed to be intermediate between production traits and survival. Figure 2 shows a strong relationship between the accuracy of prediction of DGV and the heritability of the trait. Predictions of production traits with a higher heritability, such as protein percentage (h2 = 0.56), fat percentage (h2 = 0.52), and milk yield (h2 = 0.28), were more accurate than predictions of traits with a lower heritability, such as overall type (h2 = 0.18) and survival (h2 = 0.03).
Figure 2 shows a consistent trend in the accuracy of DGV when the SNP density decreased from 42,576 to approximately 1,000 SNP using trait-depended subsets of SNP. When SNP density exceeded 1,000 SNP the accuracy of DGV reached a plateau, and increases in accuracy with increasing number of SNP were marginal or fluctuated around the maximum accuracy (Table. 2). This plateau in accuracy of DGV was consistent in both bulls and cows (Figure 2). At densities below 1,000 SNP accuracies declined relatively rapidly, subsets of 100 SNP consistently showed the lowest accuracy within the range examined here (Figure 2).
When the number of SNP in the subset was 1,000 or larger, using trait-specific subsets gave higher accuracies than using a common subset of SNP in both validation sets, with the exception of overall type for both bulls and cows (Figure 3 and 4). In addition, the rate of decrease in accuracy, with respect to the size of the subset, was much more rapid for evenly spaced SNP than for trait-dependent SNP. The rate of decrease in accuracy tended to be lower for production traits, which have a higher heritability than traits related to fitness. Predictions based on at least 1,000 or 3,000 SNP appeared to be very robust to how SNP were selected, but were very sensitive when the subset included fewer SNP.
For the overall type trait, subsets including more than 1,000 of the highest ranked SNP for the trait gave lower accuracies than evenly spaced SNP selected for ASI and APR, which might be due the smaller number of training records available for this trait. All subsets containing less than 500 SNP performed poorly for survival, which has a low heritability (h2 = 0.03), particularly subsets of SNP selected for APR and ASI.
Summary of accuracy of DGV using low-density assays derived by PLSR
Number of SNP
Common assay of evenly spaced SNP
Common assay of evenly spaced SNP
The objective of the study was to evaluate the use of low-density SNP assays for genomic selection of dairy cattle. As also shown by Weigel et al.  for a single trait, the accuracy of DGV decreased with decreasing number of SNP in the subsets. However, a low-density assay comprising selected SNP can deliver a substantial portion of the gain of a high-density assay, even if a common set of SNP is used across traits. Our results show small differences between RR and PLSR when using high-density assays, but differences between the two methods become more evident for subsets containing fewer SNP.
Recently, a number of studies have reported on the accuracy of DGV for dairy traits [16–18, 20, 29–32]. These have shown that the accuracy of DGV depends on the size of the training data, SNP density, heritability and the genetic relationships between animals in the training and validation data. Although it is difficult to compare accuracies between studies, accuracies estimated in the current study are within the range of those reported previously.
There was a strong relationship between the accuracy of prediction and the heritability of the trait, with the prediction for production traits, which had with a higher heritability, being more accurate than that for traits with a low heritability. The generally low accuracies of DGV for survival are perhaps in part due to its low heritability (h2 = 0.03, ) and the low number of effective records contributing to the DRBV for young bulls (Table 1). For a trait with a low heritability, achieving an accuracy similar to that obtained for production traits requires more records [18, 33, 34]. Results for the overall type trait were less consistent across the various analyses, with larger differences between bulls and cows and between subset selection strategies compared to other traits. The differences between cows and bulls for overall type can be partly attributed to the fact that the cow's phenotype is derived from a single observation, and the smaller number of animals in the training and validation sets may be responsible for some of the variation between methods of SNP selection. In general, the estimated accuracies reported herein most likely underestimate the correlation between DGV and true breeding value, as the phenotypes (DRBV, DTD and TD) are not perfectly predicting the true breeding value.
Both, RR and PLSR performed very similar in predicting DGV and differences were generally small. However, the highest accuracy of prediction of PLSR was obtained with subsets including considerably fewer SNP than the high-density assay and fewer SNP than the best subset for RR. This might indicate that using PLSR is less appropriate when analysing very large numbers of SNP, although the differences between the maximum accuracy of DGV and the accuracy obtained with 42,576 SNP was small. A similar result has been found by Solberg et al.  who have compared PLSR and BayesB for different maker densities in simulated data and found that BayesB gives higher accuracies than PLSR and that the largest difference is obtained with high marker densities. In other simulation studies, Meuwissen et al.  and Habier et al.  have found higher accuracies for BayesB compared to RR. In all three simulation studies, a limited number of QTL with large effects accounts for most of the genetic variance. This situation is similar to the distribution of QTL effects for fat percentage, where a mutation in the gene DGAT1  is segregating which accounts for 30% of the genetic variance in our population. Of the 300 highest ranked SNP for fat percentage, 11 were located on BTA14 in the region of DGAT1, with the SNP with rank 1 closest to the known mutation. The highest accuracy for fat percentage was obtained with subsets including substantially less SNP than the high-density assay and this suggests that part of the advantage of BayesB over PLSR and RR in the simulations stems from the fact that it simultaneously performs shrinkage of marker coefficients and marker selection .
Comparisons of accuracies across traits between validation sets of cows and bulls were constrained by the fact that for cows the accuracy of DGV prediction, computed as the correlation between DGV and DRBV, could not be calculated for five out of the nine traits, as DRBV information was not available for cows. A possible remedy would be to use the correlation between DGV and estimated breeding value, r(DGV, EBV), as a measure of accuracy instead. When we computed r(DGV, EBV) in bulls and cows (results not shown) we found higher values of r(DGV, EBV) in cows than in bulls, and using r(DGV, EBV) would considerably overestimate the accuracy of genomic selection of cows. The higher correlation between DGV and EBV in cows is probably due to a larger contribution of shared information to the EBV through pedigree relationships.
Pedigree relationships between training and validation animals might also have contributed to the accuracy of DGV in cows, with 473 (92.7%) of cows having their sire in the training set compared to 80.8% of the young bulls, and the training set containing more sires of cows (N = 164) than sires of young bulls (N = 30). This could partly explain why for protein, fat and milk yield, the r(DGV, TD) of cows was similar in size as the r(DGV, DTD) for bulls, although the phenotypes of bulls are derived from many more daughter records than the cow phenotypes. Furthermore, the records of cows are included in the sire phenotype. However, the effect of cows contributing information to the sire phenotype is expected to be small, with bulls in the training set which sired a cow in the validation set having on average 5,259 phenotyped daughters. A larger variance of the phenotypes for cows compared to the pre-selected bull teams (Figure 1) has also positively contributed to the correlation between DGV and phenotypes for cows. In addition, the cows were selected from a repository of animals which have been well recorded, so the heritabilities in this subset are most likely higher than in the wider industry.
Habier et al.  have demonstrated that the maximum of the additive-genetic relationships between training and validation animals is a good indicator for the accuracy of DGV. The additive-genetic relationship did differ for bulls and cows whose sires were or were not represented in the training set, as shown in Figure 6. As part of the genetic relationship can be captured by SNP , one would expect higher accuracies for animals whose sire is included in the training set [5, 16]. Our results are less conclusive, as accuracies of DGV for bulls whose sire was included in the training set was not higher than those for bulls whose sire was not included in the training set for all traits. A partial explanation might be the relatively small number of bulls whose sire was not genotyped, resulting in a large sampling error of the correlation between DGV and phenotype. However, an important observation is that the differences in accuracies between the two groups appear to be independent of the number of SNP and the method of SNP selection [see Additional files 2 and 3], with the exception of survival and overall type and possible reasons for this are discussed above.
In practise, whilst an attractive application of genotyping using low-density assays is the selection of replacement heifers, for reasons given above, the accuracy is expected to be smaller than reported here for cows. However, the accuracy of genomic selection will be increased if DGV predictions are combined with information from pedigree [17, 18, 20].
Few SNP were in common between the trait-specific subsets (Figure 5) and, given that at least 1,000 SNP are required to obtain accurate DGV predictions for most traits, combining the highest ranked SNP for each trait onto a single chip or developing multiple low-density assays might not provide adequate reductions in genotyping costs. Irrespective of the method of SNP selection, subsets of 3,000 evenly spaced SNP provided more than 90% of the accuracy that can be achieved with a high-density assay in genomic selection of cows and 80% of the high-density assay in young bulls. Furthermore, the rate of increase in accuracy with increasing size of the subset was more rapid for evenly spaced SNP, so that the additional gain from using trait-specific assays or SNP related to a single index such as ASI or APR was small for subsets with a larger number of SNP.
Predictions using subsets including 3,000 of the highest ranked SNP were only 1.06 times more accurate in bulls and 1.03 times more accurate in cows than a common subset of evenly spaced SNP of the same size selected based on MAF. This suggests that the distribution of true effects is more or less spread among many loci across the genome and that 3,000 evenly spaced SNP largely capture the level of LD present in the population. While accuracies based on at least ~1,000 or 3,000 SNP appear to be very robust to the method used to select those SNP, accuracies are very sensitive when fewer than 1,000 SNP are used (Figure 3 and 4). However, without re-estimation of effects, accuracy when using evenly spaced low-density assays is expected to decrease steadily and faster over generations compared to accuracy from a high-density assay or from subsets including SNP with large effects, but this loss could be limited by genotyping the parents used for breeding with the high-density assay and retraining . Dense genotyping of some parents might not be necessary if the genotype information of the high-density assay can be imputed from a low-density SNP panel [36, 37].
Weigel et al.  have also assessed the ability to predict DGV using subsets of SNP with largest effects and subsets of evenly spaced SNP for the trait lifetime net merit in dairy cattle. Although it is difficult to compare accuracies between studies due to differences in the methods used to calculate DGV, the size of the training data and the number of SNP of the high-density assay, both studies agree well in that a trait-specific subset including 2,000 of the highest ranked SNP captured most of the gain achieved with a high-density assay. However, for subsets of evenly spaced SNP, the rate in loss of accuracy with decreasing SNP numbers was lower in our study compared to . Here, selection of SNP was performed by choosing the highest ranked SNP within segments of approximately equal size, whereas in  spacing of SNP was only informed by the position of the SNP. Using the latter approach, one would expect a larger number of low ranked SNP to be selected in the assay. Indeed, when compared to subsets including evenly spaced SNP selected for APR, we found that selecting SNP only on their position reduced the relative accuracy between 5% and 17% depending on the size of the subset and the trait analysed (results not shown).
Although ranking of SNP was based on the magnitude of the estimated SNP effects, SNP selected on their rank had a higher average minor allele frequency (p MAF) than all SNP on the high-density assay. For example, subsets containing 300 evenly spaced SNP selected for ASI by PLSR had a mean p MAF = 0.33 compared to p MAF = 0.30 for SNP selected by RR and p MAF = 0.27 for SNP selected on only their location. This suggests that selection of SNP should be based on their expected contributions to the genetic variance, which is a function of the allele frequency in the training set. This also implies that accuracy of prediction will be lower for a validation set where the distribution of allele frequencies does not resemble that of the training set and why prediction equations derived in one breed do not predict accurate DGV when applied to other breeds, as shown by Hayes et al. .
Arguably, the major advantage of a low cost SNP assay will be that training sets will become much larger, as relatively more animals are genotyped and hence accuracy of DGV will increase (e.g.). The current reference populations predominantly consist of elite progeny tested sires and to significantly increase the size of the training data will require the genotyping of cows.
Genomic selection has become a routine in dairy cattle breeding programs worldwide. The current cost of whole-genome selection based on dense SNP genotypes has limited the application to the selection of elite males and females that are likely to become parents of the next generation. Results of our study indicate that accurate genomic evaluation of the broader bull and cow population can be achieved with genotyping assays containing ~ 3,000 to 5,000 SNP. A chip containing 3,000 evenly spaced markers can provide approximately 90% of the accuracy achieved with a high-density SNP assay for genomic selection of bulls and cows combined across traits. Possible applications include the selection of replacement heifers and the pre-screening of young bulls and potential bull dams. Assays with evenly spaced markers are preferable as they can be used across traits and possibly across populations. It also allows for a high volume generic chip to be produced, which will lower assay cost per individual and will limit heterogeneity of genomic information compared to using multiple assays for different traits. Evenly spaced low-density assays might also permit the reconstruction of the genotype information of high-density assays through imputation, which is important in situations where, for example, high-density genotyping is limited to nucleus breeding herds. Increasing the proportion of animals genotyped will further increase the accuracy of genomic selection as the training data grows over time, particularly through genotyping of cows.
The authors wish to thank Genetics Australia for semen samples, the Australian Dairy Herd Improvement Scheme (ADHIS) for providing phenotype and pedigree data and Phillip Bowman for the extraction of the data. The study was supported by the Dairy Futures Cooperative Research Centre (CRC).
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.