Research Article| Volume 98, ISSUE 12, P9051-9059, December 2015
• PDF [182 KB]PDF [182 KB]
• Top

# Including different groups of genotyped females for genomic prediction in a Nordic Jersey population

Open AccessPublished:September 30, 2015

## Abstract

Including genotyped females in a reference population (RP) is an obvious way to increase the RP in genomic selection, especially for dairy breeds of limited population size. However, the incorporation of these females must be conducted cautiously because of the potential preferential treatment of the genotyped cows and lower reliabilities of phenotypes compared with the proven pseudo-phenotypes of bulls. Breeding organizations in Denmark, Finland, and Sweden have implemented a female-genotyping project with the possibility of genotyping entire herds using the low-density (LD) chip. In the present study, 5 scenarios for building an RP were investigated in the Nordic Jersey population: (1) bulls only, (2) bulls with females from the LD project, (3) bulls with females from the LD project plus non-LD project females genotyped before their first calving, (4) bulls with females from the LD project plus non-LD project females genotyped after their first calving, and (5) bulls with all genotyped females. The genomically enhanced breeding value (GEBV) was predicted for 8 traits in the Nordic total merit index through a genomic BLUP model using deregressed proof (DRP) as the response variable in all scenarios. In addition, (daughter) yield deviation and raw phenotypic data were studied as response variables for comparison with the DRP, using stature as a model trait. The validation population was formed using a cut-off birth year of 2005 based on the genotyped Nordic Jersey bulls with DRP. The average increment in reliability of the GEBV across the 8 traits investigated was 1.9 to 4.5 percentage points compared with using only bulls in the RP (scenario 1). The addition of all the genotyped females to the RP resulted in the highest gain in reliability (scenario 5), followed by scenario 3, scenario 2, and scenario 4. All scenarios led to inflated GEBV because the regression coefficients are less than 1. However, scenario 2 and scenario 3 led to less bias of genomic predictions than scenario 5, with regression coefficients showing less deviation from scenario 1. For the study on stature, the daughter yield deviation/daughter yield deviation performed slightly better than the DRP as the response variable in the genomic BLUP (GBLUP) model. Therefore, adding unselected females in the RP could significantly improve the reliabilities and tended to reduce the prediction bias compared with adding selectively genotyped females. Although the DRP has performed robustly so far, the use of raw data is recommended with a single-step model as an optimal solution for future genomic evaluations.

## Introduction

The size of the reference population (RP) is one of the important factors influencing the accuracy of genomic prediction (
• Goddard M.E.
• Hayes B.J.
Mapping genes for complex traits in domestic animals and their use in breeding programmes.
). To date, RP mainly has consisted of proven bulls in national or international dairy cattle genomic selection programs (
• Van Tassell C.P.
• Wiggans G.R.
• Sonstegard T.S.
• Schnabel R.D.
• Taylor J.F.
• Schenkel F.S.
Invited review: Reliability of genomic predictions for North American Holstein bulls.
;
• Harris B.L.
• Johnson D.L.
Genomic predictions for New Zealand dairy bulls and integration with national genetic evaluation.
;
• Jorjani H.
• Zumbach B.
• Dürr J.
• Santus E.
Joint Genomic Evaluation of BSW Populations.
;
• Muir B.
• Doormaal B.V.
• Kistemaker G.
International Genomic Cooperation—North American Perspective.
;
• Lund M.S.
• de Roos A.P.W.
• de Vries A.G.
• Druet T.
• Ducrocq V.
• Fritz S.
• Guillaume F.
• Guldbrandtsen B.
• Liu Z.T.
• Reents R.
• Schrooten C.
• Seefried F.
• Su G.S.
A common reference population from four European Holstein populations increases reliability of genomic predictions.
;
• Gao H.
• Su G.
• Janss L.
• Zhang Y.
• Lund M.S.
Model comparison on genomic predictions using high-density markers for different groups of bulls in the Nordic Holstein population.
). Due to the decreasing costs of genotyping and the increasing exchange of data on genotyped proven bulls, the prediction accuracies have been markedly enhanced due to the increased RP. This process has been beneficial for Holsteins, a widespread breed that is present in many countries. However, for the numerically smaller and geographically less widespread dairy breeds, such as Nordic Jersey or Nordic Red Cattle, the advantage of sharing reference data has been limited to international collaboration. An alternative solution to the problem of few proven bulls is to increase the RP by including genotyped females (heifers and cows) even though the information from females is much less reliable compared with information from bulls that have been progeny tested using a large daughter group.
Prediction accuracy is expected to be enhanced by increasing the size of the RP. However, because the number of progeny from bulls being tested is shrinking due to the use of genomic selection as a pre-selection tool for young bulls entering the progeny-testing scheme, the RP will increase less rapidly over time (
• Schaeffer L.R.
Strategy for applying genome-wide selection in dairy cattle.
;
• Lillehammer M.
• Meuwissen T.H.E.
• Sonesson A.K.
A comparison of dairy cattle breeding designs that use genomic selection.
). A simulation study by
• Thomasen J.R.
• Egger-Danner C.
• Willam A.
• Guldbrandtsen B.
• Lund M.S.
• Sorensen A.C.
Genomic selection strategies in a small dairy cattle population evaluated for genetic gain and profit.
showed that the inclusion of genotyped cows in the RP was an efficient way to increase the genetic gain and would be a profitable investment for the breeding schemes of small breeds. Furthermore, since 2010, genotyping of females with a low-density (LD) chip has been implemented at a large scale in Holsteins in the United States.
Therefore, an appealing and cost-effective approach could be to genotype females using an LD chip such as the Illumina BovineLD BeadChip (http://support.illumina.com/array/array_kits/bovineld_dna_analysis_kit.html), followed by the imputation to higher density (
• Browning B.L.
• Browning S.R.
A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals.
;
• Dassonneville R.
• Brondum R.F.
• Druet T.
• Fritz S.
• Guillaume F.
• Guldbrandtsen B.
• Lund M.S.
• Ducrocq V.
• Su G.
Effect of imputing markers from a low-density chip on the reliability of genomic breeding values in Holstein populations.
).
The advantage of genotyping females in dairy cattle breeding has been reported in some previous studies. The first empirical study of the inclusion of cows in the RP was reported by
• Wiggans G.R.
• Cooper T.A.
• Cole J.B.
Technical note: adjustment of traditional cow evaluations to improve accuracy of genomic predictions.
, where the records of the genotyped cows were pre-adjusted to be comparable with those of the genotyped bulls; these authors found an average gain in reliabilities of 3.5 and 0.9 percentage points in Holstein and Jersey populations, respectively.
• Pryce J.
• Hayes B.
• Goddard M.
Genotyping dairy females can improve the reliability of genomic selection for young bulls and heifers and provide farmers with new management tools.
demonstrated an improvement of 8 percentage points in the GEBV reliabilities by adding 10,000 genotyped cows to an RP consisting of approximately 3,000 bulls.
• Bapst B.
• Baes C.
• Seefried F.R.
• Bieber A.
• Simianer H.
• Gredler B.
Effect of cows in the reference population: First results in Swiss Brown Swiss.
added approximately 1,236 genotyped cows to an existing RP consisting of 4,085 bulls in a Brown Swiss population but did not achieve a significant improvement in the accuracy of genomic prediction. In the United States, 30,852 genotyped Holstein cows were incorporated into the RP of 21,883 Holstein bulls, and an extra 0.4 percentage points of genomic reliability was observed when averaged across all traits (
• Cooper T.
• Wiggans G.
Including cow information in genomic prediction of Holstein dairy cattle in the US.
). In general, the outcomes of adding genotyped females to the RP appeared to vary among different implementations. The value of adding cows to the RP mainly appears to be dependent on the proportion of added genotyped cows and the size of the original bull RP.
Deregressed proof (DRP), which is a back-calculation of phenotypes from EBV using the reliabilities obtained from the traditional genetic evaluation, has been widely adopted nationally and internationally as the pseudo-phenotype of choice in genomic evaluation procedures. This method is used due to the advantages of simplification over the daughter yield deviation (DYD) and nonregressed property compared with EBV (
• Gao H.
• Lund M.S.
• Zhang Y.
• Su G.
Accuracy of genomic prediction using different models and response variables in the Nordic Red cattle population.
). Therefore, DRP has worked fairly well when observations in the RP consisted of genotyped and progeny-tested bulls (
• Garrick D.J.
• Taylor J.F.
• Fernando R.L.
Deregressing estimated breeding values and weighting information for genomic regression analyses.
;
• Lund M.S.
• de Roos A.P.W.
• de Vries A.G.
• Druet T.
• Ducrocq V.
• Fritz S.
• Guillaume F.
• Guldbrandtsen B.
• Liu Z.T.
• Reents R.
• Schrooten C.
• Seefried F.
• Su G.S.
A common reference population from four European Holstein populations increases reliability of genomic predictions.
;
• Gao H.
• Lund M.S.
• Zhang Y.
• Su G.
Accuracy of genomic prediction using different models and response variables in the Nordic Red cattle population.
). However, for the genotyped females, the reliabilities of the DRP are much lower compared with the reliability of DRP for progeny-tested bulls. In such cases, the alternatives could be to use the yield deviation (YD) or raw phenotypic data in place of the DRP as the response variable for genotyped females because the YD is generated based only on the cows’ own records and avoids the deregression procedure. An untested hypothesis that we address in this study regards the importance of the choice of the response variable for the genotyped females as a factor influencing the accuracy of genomic prediction.
The breeding organizations in Denmark, Finland, and Sweden (Viking Genetics) have initiated a female genotyping project with the chance of genotyping all heifers in entire selected herds using a LD chip. To assess the effect of including genotyped females in the RP, the first batch of data from this project was used. The purposes of this study were to (1) examine the effect of adding different sources of genotyped females to the RP for Nordic Jersey, and (2) explore the effect of different response variables for the genotyped females on prediction accuracy.

## Materials and Methods

### Data

The animals used in this study consisted of 1,414 genotyped Nordic Jersey bulls born between 1981 and 2011 (with several individuals that were missing DRP). In addition, 1,154 proven Jersey bulls in the United States that were born between 1950 and 2009 were added to the RP through the collaboration for exchanging reference data to maximize the number of progeny-tested bulls in the RP (
• Su G.
• Nielsen U.
• Wiggans G.
• Aamand G.
• Guldbrandtsen B.
• Lund M.
Improving genomic prediction for Danish Jersey using a joint Danish-US reference population.
). A total of 4,251 genotyped females with DRP were classified into different subsets based on the genotyping strategy. Overall, 3,492 females born after 2010 were phenotyped and genotyped, along with the entire herd, using an LD chip through the LD project (hereafter referred as LD females). The remaining individuals, consisting of 759 genotyped and phenotyped females, were selected for genotyping by private breeders according to their individual breeding programs (hereafter referred as non-LD females). Depending on the genotyping date, among these 759 non-LD females, 240 heifers born between 2002 and 2011 were genotyped before their first calving (hereafter referred as non-LD heifers), whereas 143 cows born between 2002 and 2011 were genotyped after their first calving (hereafter referred as non-LD cows), the rest were 376 cows with unknown first calving date (hereafter referred as other cows). Therefore, non-LD heifers were genotyped before their own performance records were obtained, and therefore, they could be considered randomly selected individuals or selected based on parent average as opposed to non-LD cows, which could have been selected due to preferential treatment, thereby possibly resulting in a prediction bias.
Consequently, 5 scenarios for constructing an RP were considered in this study: (1) bulls only, (2) bulls with females from the LD project, (3) bulls with females from the LD project plus non-LD heifers, (4) bulls with females from the LD project plus non-LD cows, and (5) bulls with all genotyped females. For the scenarios that included genotyped females, the genotyped daughters of validation bulls were removed from the RP. The exact numbers of bulls and cows in the RP for the 5 scenarios are shown in Table 1.
Table 1Number of genotyped individuals in the reference population for each scenario
The numbers after removing the genotyped daughters of validation bulls from the RP.
ScenarioMalesFemales
Nordic

bulls
US

bulls
LD
Females that were genotyped with a low-density (LD) chip.

females
Non-LD

heifers
Non-LD females that were genotyped before their first calving.
Non-LD

cows
Non-LD females that were genotyped after their first calving.
Other

cows
Cows with an unknown first calving date.
11,0531,154
21,0531,1542,431
31,0531,1542,431214
41,0531,1542,431143
51,0531,1542,431214143332
1 The numbers after removing the genotyped daughters of validation bulls from the RP.
2 Females that were genotyped with a low-density (LD) chip.
3 Non-LD females that were genotyped before their first calving.
4 Non-LD females that were genotyped after their first calving.
5 Cows with an unknown first calving date.
The Nordic Jersey bulls and the non-LD females were genotyped using the Illumina Bovine SNP50 BeadChip (54,001 SNP for v1 and 54,609 SNP for v2; Illumina, San Diego, CA). Jersey bulls from the United States were genotyped with the Illumina Bovine SNP50 chip or the GeneSeek Genomic Profiler (76,999 SNP), and the extra SNP that are not in the Illumina Bovine SNP50 chip were removed for the prediction. The LD females were genotyped using the Illumina Bovine LD v1.1 BeadChip (6,909; Illumina), and the data were converted to regular 50k data. The Beagle package (
• Browning B.L.
• Browning S.R.
A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals.
) was applied for imputation, and the imputation from LD to 50k data in the Danish Jersey yielded an average allelic R2 (squared correlation between true and imputed genotype) value of 0.96 from Beagle (R. F. Brøndum, Aarhus University, Tjele, Denmark; personal communication). A total of 41,647 SNP remained after implementing the requirement for a minor allele frequency of at least 0.1%. Eight traits (sub-indices for milk, fat, and protein yields; mastitis; body conformation; udder conformation; milking speed; and yield index) in the Nordic total merit system were assessed. A cut-off year of 2006 was used to partition the Nordic Jersey bulls into the RP and validation population (VP) according to the birth date of the individuals, which means that the RP contained individuals born before January 1, 2006.

### Genomic Prediction Model

The following genomic BLUP (GBLUP) model (
Efficient methods to compute genomic predictions.
;
• Hayes B.J.
• Visscher P.M.
• Goddard M.E.
Increased accuracy of artificial selection by using the realized relationship matrix.
) was applied to predict GEBV in the current study:
$y = μ1 + Zg + e,$

where y is the vector of the DRP of the genotyped individuals in the RP; μ is a general mean; 1 is a vector of ones; Z is an incidence matrix allocating records to breeding values; g is a vector of GEBV with var(g) = $Gσg2,$ in which $σg2$ is the additive genetic variance (obtained from the routine genetic evaluation for the present study); and G is the realized genomic relationship matrix created using the method from (
Efficient methods to compute genomic predictions.
). That is, G = $(M − P)(M − P)′/∑i=1m2pi(1−pi),$ where M is an n × m matrix (number of individuals × number of loci) with SNP coded as 0, 1, and 2 for genotypes 11, 12, and 22, respectively. P is an n × m matrix containing twice the allele frequencies expressed as Pi = 2pi, where pi is the allele frequency of the homozygous genotype coded with 2 for all genotyped individuals at locus i. In principle, pi should be estimated from the unselected base population. However, in practice, this information is not available, and instead, pi was calculated based on the observed data in this study. Furthermore, (MP)(MP) was divided by $∑i=1m2pi(1−pi)$ to normalize the G matrix to be compatible with the traditional numerator relationship matrix. For random residuals, it is assumed that $e~N(0,Dσe2),$ where$σe2$ is the residual variance, and D is a diagonal matrix with the elements of reciprocals of weights. The weight factor used in the current data was $rDRP2/(1−rDRP2),$ as explained in
• Gao H.
• Christensen O.F.
• Nielsen U.S.
• Zhang Y.
• Lund M.S.
• Su G.
Comparison on genomic predictions using three GBLUP methods and two single-step blending methods in the Nordic Holstein population.
. The DRP were derived from the February 2014 evaluation by Nordic Cattle Genetic Evaluation (NAV; http://www.nordicebv.info/Routine+evaluation/). The deregression procedure was performed by the iterative method described in
• Jairath L.
• Dekkers J.C.
• Schaeffer L.R.
• Liu Z.
• Burnside E.B.
Genetic evaluation for herd life in Canada.
and
• Schaeffer L.R.
Multiple trait international bull comparisons.
using the MiX99 package (
• Strandén I.
• Mäntysaari E.M.
A recipe for multiple trait deregression.
). All analyses of the GBLUP models were conducted using the DMU package (
• Jensen J.
An User's Guide to DMU, Version 6, Release 5.1. Center for Quantitative Genetics and Genomics.
).

### Validation

To validate the genomic prediction and compare different scenarios, a maximum of 190 genotyped Nordic Jersey bulls born in 2006 and onward with DRP were used in the VP. A few of these bulls did not have an evaluation for all traits. Validation reliability of genomic prediction was estimated as the squared Pearson correlation coefficient between GEBV and DRP $(cor(GEBV,DRP)2)$ divided by the average reliability of the DRP of the bulls in the VP. The regression coefficient of the DRP on the GEBV was used as a measurement of the bias of the GEBV. Unbiased predictions are expected to have a regression coefficient of 1, whereas values smaller than 1 indicate the inflation in the variance of GEBV, and values larger than 1 indicate the deflation of the variance of GEBV (
• Reverter A.
• Golden B.L.
• Bourdon R.M.
• Brinks J.S.
Technical note: Detection of bias in genetic predictions.
). Confidence intervals of reliability and regression coefficients were estimated using a bootstrap procedure by resampling the entire original VP 2,000 times with replacement.

### Comparison of Response Variables

As shown in Table 2, the mean reliabilities of the DRP for cows (0.23–0.44) were less than half of the reliabilities of bulls (0.61–0.83) across all traits, indicating the uncertainty, which was also noted in terms of the weights applied in the prediction model when including females in the RP.
Table 2Mean and standard deviation of deregressed proof (DRP) reliabilities of bulls and cows based on scenario 2
DRP data with a reliability value less than 0.20 were removed from the reference population.
Traith
Heritabilities (h2) were determined by Nordic Cattle Genetic Evaluation (http://www.nordicebv.info/Routine+evaluation/).
n
Heritabilities (h2) were determined by Nordic Cattle Genetic Evaluation (http://www.nordicebv.info/Routine+evaluation/).
MeanSD
BullsCowsBullsCowsBullsCows
Milk0.392,2072,4230.830.440.100.06
Fat0.392,2072,3600.790.370.120.05
Protein0.392,2072,4230.830.440.100.06
Mastitis0.042,1781,2560.630.230.210.02
Body conformation0.302,0242,0010.750.260.090.03
Udder conformation0.252,0232,0010.720.260.090.03
Milking speed0.261,0781,4640.610.300.170.02
Yield index0.392,2072,4230.830.440.100.06
1 DRP data with a reliability value less than 0.20 were removed from the reference population.
2 Heritabilities (h
Heritabilities (h2) were determined by Nordic Cattle Genetic Evaluation (http://www.nordicebv.info/Routine+evaluation/).
) were determined by Nordic Cattle Genetic Evaluation (http://www.nordicebv.info/Routine+evaluation/).
To investigate the effect of response variables on the prediction accuracies when including the genotyped females in the RP, a study was conducted using raw data, the DYD/YD and the DRP for a selected model trait. Considering the simplicity of a single recording per individual, stature, which is one of the sub-traits in the combined index of body conformation, was used as the model trait in this comparison. Performance data collected for stature were used for the evaluation run in February 2015, with 184,905 first lactation records of females remaining after removing the daughters of the validation bulls. To compute the DYD/YD, a program was used to compute trait deviations and progeny trait deviations in the Gibbs sampler module (RJMC) in the DMU package (
• Jensen J.
An User's Guide to DMU, Version 6, Release 5.1. Center for Quantitative Genetics and Genomics.
). This program computes trait deviations and progeny trait deviations for each stored sample. Traditional BLUP keeps the dispersion parameter constant and only samples the location parameter. Consequently, the DYD/YD were obtained from the posterior means of progeny trait deviations and trait deviations along with their posterior standard deviations (i.e., a measure of accuracy without any approximation). For the implementation of the Markov chain Monte Carlo method, a single chain was run for 4,000 iterations, with a burn-in of 1,000 iterations and a thinning interval of 10. Variance components from the official NAV routine were used as the fixed dispersion parameters. The computation needs to be conducted separately for bulls and cows because the information on the genotyped daughters has to be removed when calculating the DYD for the bulls to avoid double counting. After the cut-off, the RP consisted of 1,049 genotyped Nordic Jersey bulls and 4,024 genotyped cows with the DYD/YD and DRP, with 230 Nordic Jersey bulls remaining with the DYD and DRP in the VP. For raw phenotypic data, the GBLUP model was replaced by a single-step model (
• Legarra A.
• Aguilar I.
• Misztal I.
A relationship matrix including full pedigree and genomic information.
;
• Christensen O.
• Lund M.
Genomic prediction when some animals are not genotyped.
) to compute the GEBV; a pedigree with 338,097 individuals was used to construct the numerator relationship matrix, and combined with the genomic relationship matrix. A detailed description of the model can be found in
• Gao H.
• Christensen O.F.
• Nielsen U.S.
• Zhang Y.
• Lund M.S.
• Su G.
Comparison on genomic predictions using three GBLUP methods and two single-step blending methods in the Nordic Holstein population.
.

## Results

Table 3 presents the reliabilities for the 5 different scenarios of the RP. Generally, a benefit of including genotyped females in the RP was observed in terms of an increase in reliability (not for body conformation). For mastitis, scenario 5 achieved a gain of 5.1 percentage points in reliability, and for udder conformation, scenario 2 had a higher reliability by 2.3 percentage points. Gains of 1.9 to 4.5 percentage points in reliability have been achieved, when using averages across the 8 traits, by including genotyped females in the RP. Averaged over the 8 traits, the addition of all genotyped females in the RP achieved the maximum gain in reliability (scenario 5), followed by scenario 3, scenario 2, and scenario 4, especially for protein, mastitis, and milking speed. The mean reliabilities of scenarios 2 and 3 were 1.7 and 1.1 percentage points lower, respectively, than for scenario 5. For scenario 3, the reliabilities of production traits were further improved with additional genotyped cows included in the RP in contrast to scenario 2. With respect to scenario 4, which has a smaller RP compared with scenario 5, the mean reliability was 2.6 percentage points lower than that in scenario 5. Furthermore, a decrease of 0.9 percentage points in reliability was observed when including the extra 143 genotyped cows in scenario 4 compared with scenario 2.
Table 3Reliabilities of the genomically enhanced breeding value (GEBV) for animals in the validation population (VP)
TraitVP

(no.)
Scenario 1
Bulls only.

(90% CI)
Confidence intervals of 90% of the reliability using the bootstrap analysis.
Scenario 2
Bulls + females genotyped with a low-density (LD) chip (LD females).

(90% CI)
Scenario 3
Bulls + LD females + non-LD females with genotyping before first calving.

(90% CI)
Scenario 4
Bulls + LD females + non-LD females with genotyping after first calving.

(90% CI)
Scenario 5
Bulls + all genotyped females.

(90% CI)
Milk1880.443 (0.344–0.543)0.481 (0.378–0.583)0.497 (0.399–0.598)0.478 (0.377–0.581)0.496 (0.399–0.595)
Fat1880.233 (0.128–0.346)0.260 (0.143–0.384)0.269 (0.150–0.390)0.251 (0.137–0.374)0.268 (0.153–0.391)
Protein1880.317 (0.213–0.427)0.327 (0.210–0.450)0.348 (0.224–0.477)0.330 (0.214–0.451)0.356 (0.230–0.483)
Mastitis1900.362 (0.247–0.484)0.366 (0.250–0.479)0.362 (0.251–0.477)0.337 (0.227–0.454)0.413 (0.300–0.521)
Body conformation1690.374 (0.270–0.483)0.370 (0.259–0.493)0.370 (0.262–0.486)0.362 (0.244–0.486)0.367 (0.254–0.482)
Udder conformation1800.372 (0.248–0.504)0.395 (0.275–0.518)0.377 (0.260–0.499)0.372 (0.258–0.488)0.356 (0.247–0.471)
Milking speed1810.241 (0.133–0.355)0.364 (0.255–0.483)0.363 (0.250–0.483)0.352 (0.238–0.473)0.423 (0.313–0.538)
Yield index1880.258 (0.145–0.380)0.264 (0.144–0.392)0.283 (0.157–0.423)0.268 (0.141–0.401)0.277 (0.154–0.408)
Mean0.3250.3530.3590.3440.370
1 Bulls only.
2 Confidence intervals of 90% of the reliability using the bootstrap analysis.
3 Bulls + females genotyped with a low-density (LD) chip (LD females).
4 Bulls + LD females + non-LD females with genotyping before first calving.
5 Bulls + LD females + non-LD females with genotyping after first calving.
6 Bulls + all genotyped females.
Table 4 shows the regression coefficients of DRP on GEBV, and unbiasedness was assessed by the degree of deviation from the expected value of unity (i.e., smaller deviation indicates less bias of GEBV). In the present study, the variance of GEBV from all scenarios was inflated to some extent because most of regression coefficients were significantly below 1. The mean deviations ranged from 0.265 to 0.345 for all scenarios. For the RP with genotyped females included (scenarios 2–5), regression coefficients deviated by more from 1 for all traits except for milking speed, indicating that the GEBV predicted by including genotyped females in the RP yielded more bias than predictions from the RP with proven bulls only (scenario 1).
Table 4Regression coefficients of deregressed proof (DRP) on the genomically enhanced breeding value (GEBV) for animals in the validation population (VP)
TraitVP

(no.)
Scenario 1
Bulls only.

(90% CI)
Confidence intervals of 90% of the regression coefficient using bootstrap analysis.
Scenario 2
Bulls + females genotyped with a low-density (LD) chip (LD females).

(90% CI)
Scenario 3
Bulls + LD females + non-LD females with the genotyping date before the first calving date.

(90% CI)
Scenario 4
Bulls + LD females + non-LD females with the genotyping date after the first calving date.

(90% CI)
Scenario 5
Bulls + all genotyped females.

(90% CI)
Milk1880.783 (0.666–0.916)0.689 (0.578–0.802)0.698 (0.588–0.811)0.685 (0.580–0.800)0.676 (0.581–0.777)
Fat1880.640 (0.477–0.801)0.591 (0.441–0.741)0.594 (0.439–0.740)0.581 (0.426–0.731)0.545 (0.420–0.667)
Protein1880.675 (0.537–0.818)0.595 (0.468–0.720)0.604 (0.481–0.724)0.595 (0.467–0.714)0.583 (0.469–0.684)
Mastitis1900.822 (0.660–1.000)0.776 (0.606–0.943)0.774 (0.613–0.936)0.730 (0.567–0.891)0.823 (0.667–0.982)
Body conformation1690.747 (0.607–0.890)0.695 (0.567–0.824)0.683 (0.566–0.805)0.675 (0.552–0.805)0.653 (0.547–0.763)
Udder conformation1800.938 (0.761–1.118)0.828 (0.683–0.967)0.780 (0.644–0.913)0.779 (0.638–0.921)0.750 (0.619–0.893)
Milking speed1810.626 (0.455–0.786)0.669 (0.527–0.821)0.646 (0.499–0.795)0.646 (0.499–0.795)0.687 (0.557–0.818)
Yield index1880.648 (0.481–0.810)0.567 (0.411–0.714)0.577 (0.426–0.721)0.567 (0.415–0.710)0.527 (0.399–0.646)
Mean deviation
Mean of absolute deviation from 1.
0.2650.3240.3310.3430.345
1 Bulls only.
2 Confidence intervals of 90% of the regression coefficient using bootstrap analysis.
3 Bulls + females genotyped with a low-density (LD) chip (LD females).
4 Bulls + LD females + non-LD females with the genotyping date before the first calving date.
5 Bulls + LD females + non-LD females with the genotyping date after the first calving date.
6 Bulls + all genotyped females.
7 Mean of absolute deviation from 1.
Among the 4 scenarios with genotyped females included in the RP, scenario 2 and scenario 3 had slightly lower mean deviations from 1 compared with scenario 4 and scenario 5. For scenario 5, regression coefficients deviated the most from 1. However, for mastitis and milking speed, the smallest mean deviations from 1 were observed in scenario 5. In general, the differences in regression coefficients between scenario 2 and scenario 3 were negligible except for udder conformation and milking speed.
To test the effects of different response variables on prediction accuracy, the DRP (both bulls and cows), DYD/YD (DYD for bulls and YD for cows), YD (only cows), and raw data (only cows) were compared. Table 5 presents the reliability and regression coefficient using different response variables. The DRP and DYD/YD were compared in the GBLUP model, whereas the YD and raw data were compared based on a single-step model. Validations were carried out between different combinations of GEBV and response variables when comparing DRP and DYD/YD, as shown in the last column of Table 5. Generally, DYD/YD generated slightly better prediction results than the DRP as the response variable, although the difference is not significant. The correlation between the GEBV predicted using the DRP and the GEBV predicted using the DYD/YD was 0.961. The same pattern was observed between the YD and raw data; the correlation between the GEBV from YD and GEBV from raw data was 0.999.
Table 5Reliability (R
Confidence intervals of 90% on the reliability using bootstrap analysis.
) and regression coefficients (b) using different response variables in the prediction model for stature
VP=validation population; DRP=deregressed proof; DYD=daughter yield deviation; YD=yield deviation; DYD=daughter yield deviation; GEBV=genomically enhanced breeding value.
Response variableNo. of

VP
R
Confidence intervals of 90% on the reliability using bootstrap analysis.

(90% CI)
Confidence intervals of 90% on the reliability using bootstrap analysis.
b

(90% CI)
Confidence intervals of 90% on the regression coefficient using bootstrap analysis.
ModelValidation variables
DRP_bulls + DRP_cows2300.567 (0.478–0.655)0.843 (0.758–0.931)GBLUPDYD, GEBVDRP
DYD_bulls + YD_cows2300.610 (0.523–0.695)0.893 (0.821–0.965)GBLUPDYD, GEBVDYD
DRP_bulls + DRP_cows2300.604 (0.507–0.700)0.873 (0.777–0.968)GBLUPDRP, GEBVDRP
DYD_bulls + YD_cows2300.642 (0.552–0.732)0.914 (0.841–0.985)GBLUPDRP, GEBVDYD
YD_cows2300.658 (0.577–0.736)0.976 (0.900–1.045)Single stepDYD, GEBVYD
RAW_cows2300.658 (0.572–0.733)0.966 (0.892–1.039)Single stepDYD, GEBVRAW
1 VP = validation population; DRP = deregressed proof; DYD = daughter yield deviation; YD = yield deviation; DYD = daughter yield deviation; GEBV = genomically enhanced breeding value.
2 Confidence intervals of 90% on the reliability using bootstrap analysis.
3 Confidence intervals of 90% on the regression coefficient using bootstrap analysis.

## Discussion

In this study, Nordic Jersey data were used to compare different strategies to include genotyped females in the RP, which resulted in higher reliabilities. However, prediction bias was considerably increased with the addition of genotyped females. Scenarios 2 and 3, with LD females and randomly genotyped females added to the RP, led to better results when GEBV reliability and prediction bias were considered simultaneously, followed by scenario 5 with all the genotyped females pooled together. Scenario 4, with LD females and selectively genotyped females, resulted in the smallest improvement in reliabilities.
Generally, improvements were more profound for production traits and milking speed. The main reason for this result is that the heritabilities for production traits were higher than the heritabilities of other traits in the present study (0.39 for milk, fat, and protein), indicating that the DRP that were calculated backward from EBV had relatively higher reliabilities for genotyped cows (Table 2). Thus, the contributions from the information on the genotyped females could be more reliable for production traits. For milking speed, only half of the bulls (1,078) have DRP records compared with other traits; therefore, significant gains in reliabilities were achieved by increasing the size of the RP by approximately 1.5-fold. This result demonstrates that the extra information obtained from genotyped females has more capability to improve the prediction accuracy for traits with a smaller RP.
The results in the current study are in line with experience from the Australian Holstein population (
• Pryce J.
• Hayes B.
• Goddard M.
Genotyping dairy females can improve the reliability of genomic selection for young bulls and heifers and provide farmers with new management tools.
), from which 10,000 genotyped cows had been added to the approximately 3,000 Holstein bull RP, yielding an increase of 4 to 8 percentage points in the reliability of breeding values in the 13 traits investigated. Unlike the significant gains experienced in Australia, in which the added genotyped cows numbered more than 3 times the number of bulls in the RP, the addition of only 1,236 genotyped Brown Swiss cows to an RP of 4,085 bulls had a very small improvement on the accuracy of genomic prediction (
• Bapst B.
• Baes C.
• Seefried F.R.
• Bieber A.
• Simianer H.
• Gredler B.
Effect of cows in the reference population: First results in Swiss Brown Swiss.
). Recently, the addition of 30,852 Holstein cows in the United States to an RP with 21,883 Holstein bulls resulted in limited gains in the reliabilities due to diminishing returns when adding new information to a large RP (
• Cooper T.
• Wiggans G.
Including cow information in genomic prediction of Holstein dairy cattle in the US.
).
The deviations from the expected regression coefficient of DRP on GEBV (i.e., unity) illustrate how much prediction bias exists in each scenario. As shown in Table 4, including genotyped females resulted in more prediction bias (scenarios 2–5) because the regression coefficients deviated more from 1 compared with scenario 1. On average, the mean deviations were 5.9 to 8.0 percentage points higher compared with the bulls-only RP (scenario 1). This negative effect resulting from the inclusion of genotyped females could be explained by the fact that information from the genotyped females was biased. Among these scenarios, scenarios 2 and 3 had less prediction bias than scenarios 4 and 5 because the females in scenario 2 are from the LD project where all females in the herds were genotyped, without the selection of individuals. Furthermore, the additional females in scenario 3 were genotyped before their first calving date, which indicated that no phenotypes were observed at that time; thus, they could be close to free of selection or selected only based on parent average, which is the same mechanism employed for the genotyped bulls. Note that compared with scenario 3, the time point for genotyping the extra cows for inclusion in the RP (non-LD cows) in scenario 4 was different; that is, the cows were genotyped after their first calving, meaning that these cows might be highly selected based on their first lactation by the breeders. Therefore, it is likely that these cows were preferentially treated when chosen to be genotyped. In this case, their breeding values may be inflated. However, the differences among these 4 scenarios were very small based on the current validation data set.
Until now, no standard solution to the problem of preferential treatment has been available, which explains why most countries still leave their female populations out of the RP for genomic prediction even though a large proportion of females have already been genotyped. However, for traits that are less commonly affected by preferential treatment, such as mastitis and milking speed, the addition of all the genotyped females in the RP is not expected to induce extra prediction bias for the GEBV, which was confirmed in the present study (Table 4). Furthermore, the DRP for bulls is less reliable for those 2 traits (Table 2). A similar result was reported by
• Dassonneville R.
• Baur A.
• Fritz S.
• Boichard D.
• Ducrocq V.
Inclusion of cow records in genomic evaluations and impact on bias due to preferential treatment.
, where 2 different genotyped cow groups were classified based on whether the cows were randomly selected or not (the latter denoted as elite cows) before adding them into the RP. The validation results showed that no difference was observed for SCC, but a difference was observed for milk yield in the Normande and Holstein breeds.
Response variable is another factor worthy of concern when including females in the RP. Currently, DRP has been widely used as the pseudo-phenotype of choice in dairy cattle breeding schemes in genomic prediction by most countries. The DRP is essentially calculated backward from the traditional genetic evaluation system and is actually a simplified version of the DYD (
• Garrick D.J.
• Taylor J.F.
• Fernando R.L.
Deregressing estimated breeding values and weighting information for genomic regression analyses.
). Moreover, an important aspect of DRP is that data can be exchanged across countries because foreign EBV are available from the Multiple Across Country Evaluation system of Interbull (http://www.interbull.org/index) and can be easily integrated into the national database for deregression, which is not possible for DYD. However, considering that the EBV reliabilities of cows were far away from those of proven bulls, the deregression procedure used for bulls may not be optimal for cows.
• Wiggans G.R.
• Cooper T.A.
• Cole J.B.
Technical note: adjustment of traditional cow evaluations to improve accuracy of genomic predictions.
noted that the traditional EBV for cows is higher than the EBV and GEBV for bulls. Consequently, a strategy was proposed to adjust the mean and Mendelian sampling variance of cows to ensure that the EBV of cows were consistent with bull EBV. Using this approach, prediction accuracies were increased compared with those using unadjusted data.
To supplement the method of pre-correcting the EBV of cows, our analysis used a novel test of alternative response variables for cows such as the YD or raw phenotypes. For the pseudo-phenotypes used in the GBLUP, the DYD of bulls, combined with the YD for cows, is the best choice, where YD is derived by removing all the nongenetic effects from the cow’s own record and DYD is a weighted mean of the YD of daughters corrected by the mate’s breeding value (
• Wiggans G.
Derivation, calculation, and use of national animal model information.
). In this study, stature was used as the model trait to compare the DYD/YD and DRP. Gibbs sampler was implemented for computing the DYD/YD, the posterior SD of DYD/YD, and the posterior mean, which was obtained with the same process, thereby ensuring that optimal weights were used for each observation in the prediction. Table 5 shows that small advantages in both reliability and prediction bias have been obtained using the DYD/YD instead of the DRP as the response variable in the GBLUP model, indicating that DYD/YD performed slightly better than DRP for stature, although the difference was not significant based on 230 validation bulls in the present study. The first reason for the small improvement could be due to the weights; that is, the posterior SD used to account for the heterogeneous residual variances of the DYD/YD is considered more accurate than the effective daughter contribution used for the DRP. The second possible reason is that in the different models used for calculation of the DRP and DYD/YD, the DRP were generated from the NAV routine deregression procedure in which a sire-maternal grandsire model was applied, whereas the DYD/YD were calculated based on an animal model. The limited increase could be mainly due to the high heritability of stature, which is 0.6 in the Nordic Jersey population. The improvement might be larger for traits with medium and low heritability because the advantage of the YD over the DRP for genotyped cows would be more significant. Moreover, the comparison between the YD and raw data showed that a good alternative is the straightforward use of the females’ information through the single-step model. The same reliabilities were observed, and the difference in prediction bias was negligible when comparing YD and raw data in the single-step model. Because extra information was contributed from dams, the reliabilities were higher and the regression coefficients were closer to 1 in a single-step model compared with the DYD/YD and DRP comparison scenarios in the GBLUP. When adding female data into the breeding program, the extra step of computing the DYD/YD using raw data should be avoided to achieve the highest reliability and minimize the prediction bias via a single-step model.

## Conclusions

The inclusion of genotyped females in the RP improves the reliability of genomic prediction by 1.9 to 4.5 percentage points. The benefit is larger for production traits than for conformation traits. For numerically small dairy breeds, including females could be an efficient way to increase the benefits of genomic selection in the breeding program. The comparison among different scenarios shows that the addition of unselected females into the RP tends to reduce the prediction bias compared with adding selectively genotyped females. The analysis based on stature confirmed that the DRP works quite well as a pseudo-phenotype in genomic prediction models. However, using raw data as the response variable with a single-step model leads to a clearly higher reliability and lower prediction bias over the use of pseudo-phenotypes with the GBLUP model.

## Acknowledgments

Seges (Aarhus, Denmark), Faba Co-op (Hollola, Finland), Växa Sverige (Uppsala, Sweden), and Nordic Cattle Genetic Evaluation (Aarhus, Denmark) are thanked for supplying the data. The first author acknowledges the funding from Innovation Fund Denmark (Copenhagen), Nordic Cattle Genetic Evaluation (Aarhus, Denmark), and Aarhus University.

## References

• Bapst B.
• Baes C.
• Seefried F.R.
• Bieber A.
• Simianer H.
• Gredler B.
Effect of cows in the reference population: First results in Swiss Brown Swiss.
in: Proc. Interbull Bull, Nantes, France Interbull, Uppsala, Sweden2013: 187-191
• Browning B.L.
• Browning S.R.
A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals.
Am. J. Hum. Genet. 2009; 84: 210-223
• Christensen O.
• Lund M.
Genomic prediction when some animals are not genotyped.
Genet. Sel. Evol. 2010; 42: 2
• Cooper T.
• Wiggans G.
Including cow information in genomic prediction of Holstein dairy cattle in the US.
in: Proc. 10th World Congr. Genet. Appl. Livest. Prod, Vancouver, BC, Canada2014: 803
• Dassonneville R.
• Baur A.
• Fritz S.
• Boichard D.
• Ducrocq V.
Inclusion of cow records in genomic evaluations and impact on bias due to preferential treatment.
Genet. Sel. Evol. 2012; 44: 40
• Dassonneville R.
• Brondum R.F.
• Druet T.
• Fritz S.
• Guillaume F.
• Guldbrandtsen B.
• Lund M.S.
• Ducrocq V.
• Su G.
Effect of imputing markers from a low-density chip on the reliability of genomic breeding values in Holstein populations.
J. Dairy Sci. 2011; 94: 3679-3686
• Gao H.
• Christensen O.F.
• Nielsen U.S.
• Zhang Y.
• Lund M.S.
• Su G.
Comparison on genomic predictions using three GBLUP methods and two single-step blending methods in the Nordic Holstein population.
Genet. Sel. Evol. 2012; 44: 8
• Gao H.
• Lund M.S.
• Zhang Y.
• Su G.
Accuracy of genomic prediction using different models and response variables in the Nordic Red cattle population.
J. Anim. Breed. Genet. 2013; 130 (a): 333-340
• Gao H.
• Su G.
• Janss L.
• Zhang Y.
• Lund M.S.
Model comparison on genomic predictions using high-density markers for different groups of bulls in the Nordic Holstein population.
J. Dairy Sci. 2013; 96 (b): 4678-4687
• Garrick D.J.
• Taylor J.F.
• Fernando R.L.
Deregressing estimated breeding values and weighting information for genomic regression analyses.
Genet. Sel. Evol. 2009; 41: 55
• Goddard M.E.
• Hayes B.J.
Mapping genes for complex traits in domestic animals and their use in breeding programmes.
Nat. Rev. Genet. 2009; 10: 381-391
• Harris B.L.
• Johnson D.L.
Genomic predictions for New Zealand dairy bulls and integration with national genetic evaluation.
J. Dairy Sci. 2010; 93: 1243-1252
• Hayes B.J.
• Visscher P.M.
• Goddard M.E.
Increased accuracy of artificial selection by using the realized relationship matrix.
Genet. Res. (Camb.). 2009; 91: 47-60
• Jairath L.
• Dekkers J.C.
• Schaeffer L.R.
• Liu Z.
• Burnside E.B.
Genetic evaluation for herd life in Canada.
J. Dairy Sci. 1998; 81: 550-562
• Jorjani H.
• Zumbach B.
• Dürr J.
• Santus E.
Joint Genomic Evaluation of BSW Populations.
in: Proc. Interbull Bull., Paris, France Interbull, Uppsala, Sweden2010: 8-14
• Legarra A.
• Aguilar I.
• Misztal I.
A relationship matrix including full pedigree and genomic information.
J. Dairy Sci. 2009; 92: 4656-4663
• Lillehammer M.
• Meuwissen T.H.E.
• Sonesson A.K.
A comparison of dairy cattle breeding designs that use genomic selection.
J. Dairy Sci. 2011; 94: 493-500
• Lund M.S.
• de Roos A.P.W.
• de Vries A.G.
• Druet T.
• Ducrocq V.
• Fritz S.
• Guillaume F.
• Guldbrandtsen B.
• Liu Z.T.
• Reents R.
• Schrooten C.
• Seefried F.
• Su G.S.
A common reference population from four European Holstein populations increases reliability of genomic predictions.
Genet. Sel. Evol. 2011; 43: 43
• Jensen J.
An User's Guide to DMU, Version 6, Release 5.1. Center for Quantitative Genetics and Genomics.
Dept. of Molecular Biology and Genetics, University of Aarhus. Research Centre Foulum, Tjele, Denmark2013
• Muir B.
• Doormaal B.V.
• Kistemaker G.
International Genomic Cooperation—North American Perspective.
in: Proc. Interbull Bull, Paris, France Interbull, Uppsala, Sweden2010: 71-76
• Pryce J.
• Hayes B.
• Goddard M.
Genotyping dairy females can improve the reliability of genomic selection for young bulls and heifers and provide farmers with new management tools.
in: Proc. 38th ICAR Conf., Cork, Ireland ICAR, Rome, Italy2012: 28
• Reverter A.
• Golden B.L.
• Bourdon R.M.
• Brinks J.S.
Technical note: Detection of bias in genetic predictions.
J. Anim. Sci. 1994; 72: 34-37
• Schaeffer L.R.
Multiple trait international bull comparisons.
Livest. Prod. Sci. 2001; 69: 145-153
• Schaeffer L.R.
Strategy for applying genome-wide selection in dairy cattle.
J. Anim. Breed. Genet. 2006; 123: 218-223
• Strandén I.
• Mäntysaari E.M.
A recipe for multiple trait deregression.
in: Proc. Interbull Bull., Riga, Latvia Interbull, Uppsala, Sweden2010: 21-24
• Su G.
• Nielsen U.
• Wiggans G.
• Aamand G.
• Guldbrandtsen B.
• Lund M.
Improving genomic prediction for Danish Jersey using a joint Danish-US reference population.
in: Page 60 in Proc. 10th World Congr. Genet. Appl. Livest. Prod, Vancouver, BC, Canada2014
• Thomasen J.R.
• Egger-Danner C.
• Willam A.
• Guldbrandtsen B.
• Lund M.S.
• Sorensen A.C.
Genomic selection strategies in a small dairy cattle population evaluated for genetic gain and profit.
J. Dairy Sci. 2014; 97: 458-470
• Wiggans G.
Derivation, calculation, and use of national animal model information.
J. Dairy Sci. 1991; 74: 2737-2746
Efficient methods to compute genomic predictions.
J. Dairy Sci. 2008; 91: 4414-4423
• Van Tassell C.P.
• Wiggans G.R.
• Sonstegard T.S.
• Schnabel R.D.
• Taylor J.F.
• Schenkel F.S.
Invited review: Reliability of genomic predictions for North American Holstein bulls.
J. Dairy Sci. 2009; 92: 16-24
• Wiggans G.R.
• Cooper T.A.