Impacts of both reference population size and inclusion of a residual polygenic effect on the accuracy of genomic prediction

Background The purpose of this work was to study the impact of both the size of genomic reference populations and the inclusion of a residual polygenic effect on dairy cattle genetic evaluations enhanced with genomic information. Methods Direct genomic values were estimated for German Holstein cattle with a genomic BLUP model including a residual polygenic effect. A total of 17,429 genotyped Holstein bulls were evaluated using the phenotypes of 44 traits. The Interbull genomic validation test was implemented to investigate how the inclusion of a residual polygenic effect impacted genomic estimated breeding values. Results As the number of reference bulls increased, both the variance of the estimates of single nucleotide polymorphism effects and the reliability of the direct genomic values of selection candidates increased. Fitting a residual polygenic effect in the model resulted in less biased genome-enhanced breeding values and decreased the correlation between direct genomic values and estimated breeding values of sires in the reference population. Conclusions Genetic evaluation of dairy cattle enhanced with genomic information is highly effective in increasing reliability, as well as using large genomic reference populations. We found that fitting a residual polygenic effect reduced the bias in genome-enhanced breeding values, decreased the correlation between direct genomic values and sire's estimated breeding values and made genome-enhanced breeding values more consistent in mean and variance as is the case for pedigree-based estimated breeding values.

Results: As the number of reference bulls increased, both the variance of the estimates of single nucleotide polymorphism effects and the reliability of the direct genomic values of selection candidates increased. Fitting a residual polygenic effect in the model resulted in less biased genome-enhanced breeding values and decreased the correlation between direct genomic values and estimated breeding values of sires in the reference population.
Conclusions: Genetic evaluation of dairy cattle enhanced with genomic information is highly effective in increasing reliability, as well as using large genomic reference populations. We found that fitting a residual polygenic effect reduced the bias in genome-enhanced breeding values, decreased the correlation between direct genomic values and sire's estimated breeding values and made genome-enhanced breeding values more consistent in mean and variance as is the case for pedigree-based estimated breeding values.

Background
With the availability of the bovine genome sequence and the development of high-density arrays of single nucleotide polymorphism (SNP) markers, the accuracy of genetic predictions has improved compared to conventional breeding value estimations based on phenotypic data and pedigree [1][2][3][4][5][6][7][8][9]. In order to model genetic variation for quantitative traits, Meuwissen et al. [10] have proposed a genetic evaluation model that includes a large number of SNP markers simultaneously. This genomic model assumes that, all the loci that affect the trait are in linkage disequilibrium (LD) with at least one SNP marker and thus marker genotypes can be used as predictors for breeding values. A main advantage of the availability of genome-enhanced breeding values (GEBV) in dairy cattle comes from the improved accuracy in pre-selecting animals for breeding. Therefore, more and more countries have been implementing genomic evaluations in dairy cattle breeding. The genomic BLUP model, which has been used to include high-density SNP data in most of the dairy cattle applications [11][12][13][14][15][16][17], assumes that all SNP contribute equally to the genetic variance, because field data results support the infinitesimal model [11,15,18].
The reliability of genomic predictions strongly depends on the number of genotyped bulls in the reference population that is used to estimate SNP effects [15,18]. The increase in genomic reliability appears to be approximately linearly correlated with the number of reference bulls [15]. However, little is known on how the size of reference populations impacts the estimation of SNP effects. A German national genomic dataset has been used to study this question. Genomic models [10,[15][16][17]19] usually assume that a given SNP marker chip, such as the Illumina Bovine54K (Illumina Inc., San Diego, CA), explains all the genetic variation of a trait, and as a consequence no residual polygenic effect (RPG) is typically fitted in genomic prediction [10,[15][16][17]19]. Fitting the RPG effect can account for the fact that SNP markers may not explain all the genetic variance [13,20,21]. Including the RPG effect in the genomic model can also render the estimates of SNP effect less biased and more persistent over generations [22]. To investigate the impact of including an RPG effect on genomic prediction, a larger dataset from the EuroGenomics reference population [18] was used. The objectives of this study were to investigate (1) the impact of the size of a genomic reference population using German reference bulls on the estimation of SNP effects and on direct genomic values (DGV) and (2) the impact of including an RPG effect on the accuracy of genomic prediction using EuroGenomics reference bulls.

German national genomic and phenotypic data
Holstein bulls from the German national genomic reference population originating partially from the national genome project GenoTrack and partially from routinely genotyped populations, were genotyped using the Illumina Bovine50k (Illumina Inc., San Diego, CA). The genotyping was conducted after ethnical review and approval by the project committee. Only SNP with a minor allele frequency greater than 1% and a call rate threshold greater than 95% i.e. 45,181 SNP were used for the analysis. Since male animals have only one allele for the 533 markers on chromosome X, the procedure to estimate marker effects developed for markers with two alleles was modified for these SNP. A genotyped animal was excluded if less than 95% of all SNP markers were called. Deregressed EBV (DRP) and effective daughter contributions (EDC) were obtained from the January 2010 German national conventional evaluation for all bulls. Forty-four traits from seven trait groups were analysed: milk production (three traits), udder health (one trait), functional longevity (one trait), calving (four traits), female fertility (six traits), workability (four traits) and conformation (25 traits). Table 1 shows the number of genotyped bulls per year of birth in the analyzed reference and validation sets. A total of 10,487 animals were genotyped. The reference bull population for milk yield comprised 5,025 German Holstein bulls. To validate the genomic evaluation system, genotyped bulls born between September 2003 and December 2004 were used for validation, and 3,676 genotyped bulls born before September 2003 were used to estimate SNP effects. To compute DGV of validation bulls, the estimated SNP effects multiplied by genotype were summed, which were then combined with the conventional pedigree index from the reference population using the pseudo-record BLUP method [14,23] to derive GEBV. Subsequently, the combined GEBV of the validation bulls were compared with their actual deregressed EBV to validate the genomic model and to check the consistency of the genetic trend and variance based on GEBV versus EBV according to the Interbull genomic validation test procedure [24]. Realised reliabilities for the pedigree-based EBV and the combined GEBV of the validation bulls were computed as the square of observed correlations with deregressed EBV, adjusted for the average reliability of the conventional EBV of their daughters [18]. The gain in reliability from genomic information was calculated as the difference between the realised reliability of the pedigree-based EBV and the combined GEBV of the validation bulls.

Scenarios to study the impact of the residual polygenic effect
To investigate the impact of including an RPG effect on GEBV, another dataset was used, which originated from the EuroGenomics collaboration [18]. This dataset comprised 17,429 genotyped Holstein bulls, representing 21.4 million daughters from the EuroGenomics countries i.e. France, Germany, Nordic countries and The Netherlands [18]. The total number of genotyped animals in the German Holstein population, including domestic candidates, was 26,191. Deregressed Multiple Across Country Evaluation (MACE) EBV from the April 2010 Interbull evaluation were used as dependent variables. In order to apply the Interbull genomic validation test [24], the genotyped bulls were divided into two groups: 14,494 reference bulls born before September 2003 and 1,377 German national validation bulls born between September 2003 and December 2004. The GEBV and parental average of pedigree-based EBV of the validation bulls were compared to their actual deregressed MACE EBV to evaluate the predictive ability of the genomic model. To investigate the impact of including an RPG effect on genomic predictions, three different percentages of residual polygenic variance to total genetic variance were considered, 5%, 10% and 15%. These three scenarios were compared to a scenario with a very small residual polygenic variance by setting the heritability of the RPG effect to 0.0001 [14], which was equivalent to 0.02% of the total genetic variance for milk yield. In order to determine the optimal residual polygenic variance for each trait in the German Holstein breed, a genomic validation study was conducted according to the Interbull genomic validation test [24], in which SNP effects were estimated using genotypic and phenotypic information of older bulls and the resulting GEBV of younger validation bulls were compared to their daughters' actual performance, i.e. deregressed EBV of the validation bulls. Observed regression coefficients of validation bulls' DRP on GEBV were compared to their expected value of 1. The scenario with observed regression coefficients close or equal to the expectation of 1 was chosen as the one with the most optimal residual polygenic variance.
In the literature [25,26], some concern has been raised that, under the BLUP genomic model, estimated SNP marker effects may model mainly family relationships. Solberg et al. [22] have suggested fitting an RPG effect to reduce this problem. In order to investigate whether incorporation of an RPG effect into the genomic model would reduce the correlation of animal DGV with EBV of sires in reference population, milk yield was analysed for the scenarios of residual polygenic variance of 0.02%, 5%, 10% and 20%.

A genomic model for German Holstein cattle
The following BLUP SNP model was applied to the DRP of reference bulls: Where q i is the DRP of bull i, μ is a general mean, ν i is the RPG effect of bull i, p is the number of fitted SNP, z ij is a genotype indicator (-1 or 1 for the two homozygotes and 0 for the heterozygote) of marker j of bull i, u j is the random regression coefficient for marker j, and e i is the residual effect of bull i. The total additive genetic variance, σ 2 a , was obtained from a conventional pedigree-based analysis, e.g. for milk production traits [6] and for female fertility traits [7], and was partitioned into two components: the residual polygenic variance σ 2 RPG = wσ 2 a , where w is the proportion of additive genetic variance explained by the RPG effect, and additive genetic variance explained by the p markers We assumed that all markers contribute equal genetic variance. The proportion of residual polygenic variance w was assumed to vary across traits. The optimal w value was determined by applying the Interbull genomic validation test [24]. Residual variance associated with the deregressed EBV q i was var(e i ) = σ 2 e /ϕ i , where σ 2 e is the error variance obtained from the pedigree-based evaluation and i is the EDC for bull i. The RPG was fitted in the same way as in conventional genetic evaluations, i.e. using full pedigree and the same grouping procedures of phantom parents [14].
Since the BLUP SNP model (1) has a large number of parameters, i.e. SNP effects that need to be estimated simultaneously, a Gauss-Seidel iteration with residual updating [27] was applied to estimate all the effects of model (1). To further improve convergence, the SNP were processed in descending order of heterozygosity.

Results and discussion
Genomic validation using German national data Table 2 shows the results of genomic validation based on the national genomic and phenotypic data of German Holstein cattle. Gains in reliability were high in general, due to the large reference population, except for fertility and calving traits. For the three milk production traits, the gain in reliability was about 30%, with the highest gain found for fat yield. Low heritability traits, such as fertility traits and stillbirth, had the lowest gain in reliability, which can be partially explained by the fact that reliabilities of conventional EBV of the reference bulls were much lower than for other traits. The realised gains in reliability of conformation traits ranged between 10% and 28%. When the genomic reference population for German Holstein cattle was switched from the German national to the EuroGenomics reference population, the number of reference bulls increased from 5,025 to 17,429. Additionally, the dependent variable DRP was derived from MACE EBV, which included phenotypic information from foreign countries, in contrast to German national EBV. In comparison to the validation results from the German national reference population in Table 2, when the larger EuroGenomics reference population was used the gain in reliability over pedigree-based EBV was 12% greater on average across four of the analyzed traits, protein yield, somatic cell score, udder depth and nonreturn rate. A significant gain in genomic reliability has also been reported in another genomic validation study using the EuroGenomics reference population [18].

Effect of the genomic reference population size
During the development of the German genomic evaluation system, a number of test runs were conducted over time, which enabled a comparison of the estimates of SNP effects across different reference populations. Table  3 shows the comparison among estimates of SNP effects for milk yield from eight genomic test runs, differing in the number of reference bulls. Because only a few young reference bulls added some daughter information over the time period of the test runs, the difference in phenotypic information on bulls already genotyped was neglected when interpreting the results in Table 3. As the number of reference bulls increased from 735 to 5,025, the observed variance of the SNP effect estimates increased more than five times. The estimate for the SNP with the largest effect increased continuously, up to 4.13 fold, as the size of the reference population increased. As expected, the correlation of SNP effect estimates was higher between any two runs, when the numbers of genotyped bulls were similar. Note that the correlation of SNP effect estimates is much lower than the correlation of DGV which was close to 1 for the reference bulls (unpublished data). It can be seen that even under the BLUP genomic model assuming equal variance for all markers, effect estimates can vary greatly between markers, and even more when new genotyped animals are added to the reference population. Table 4 shows the correlations between DGV estimates from the most recent genomic evaluations (February 2010) with the largest reference population of 5,025 bulls and DGV from each of the previous test runs. For all selection candidates, born between 2006 and 2009 and for which no phenotypic information was available, correlations between DGV increased from 0.824 to 0.993 as the number of reference bulls increased from 1,939 to 4,896. Candidates with sires included in both reference populations had somewhat higher DGV correlations than those without a genotyped sire in the reference population; however this difference in DGV correlations almost disappeared when the number of reference bulls reached 4,896. When bulls changed from candidate to reference individuals from one run to the next, the correlations between their DGV were much lower, ranging from 0.72 to 0.875, as expected. The increase in DGV correlations due to the inclusion of more reference bulls clearly shows that the genomic prediction for candidates becomes more consistent with an increasingly larger reference population.

Impact of the residual polygenic effect
Estimated SNP effects from three scenarios using the EuroGenomics reference population were compared to the scenario with the lowest residual polygenic variance for milk yield ( Table 5). The correlation of SNP effect estimates decreased only marginally with an increasing difference in residual polygenic variance assumed in the genomic model. Correlations were greater than 0.9, except for the correlation between the two most different scenarios with 0.02% and 20% residual polygenic variance (i.e. 0.86). As the residual polygenic variance increased, the variance of SNP effect estimates and the value of the estimate for the SNP with the largest effect decreased. Similar results were also obtained for all the other traits (data not shown). Table 6 shows the observed variance of estimated DGV defined as the sum of SNP marker effects and the variance of DGVt, which was defined as the sum of DGV and the estimate of the residual polygenic effect, and their correlations with conventional EBV for the reference bulls. It can be seen that the correlation between DGV and EBV decreased and the correlation between DGVt and EBV increased slightly with increasing residual polygenic variance. The variance of DGV estimates was also significantly lower for the scenarios with the higher residual polygenic variance. However, the observed variance of DGVt remained constant, indicating that the information lost from the DGV was captured by the residual polygenic effect for the reference bulls. For all scenarios, regressions of conventional EBV or DRP on DGV or RPG were unity for the reference bulls, and the regression intercepts were very close to zero (results not shown). The estimates of RPG effects and DGV were positively correlated for milk yield, with somewhat higher correlations for the scenarios with a higher percentage of residual polygenic variance, e.g. 0.42 and 0.47 for 5% and 20% residual polygenic variance respectively.
Following the Interbull genomic validation test procedure [24], conventional deregressed EBV of the validation bulls were compared to their DGV or combined GEBV estimates, which were calculated based on the reduced subset of the reference population. Table 7 shows the correlations observed between deregressed EBV, without adjusting for the reliability contributed by the daughters' performance, and DGV or GEBV estimates for the validation bulls. These correlations were high, indicating a high reliability of the genomic evaluation with 14,494 reference bulls. The correlations between DGV and deregressed EBV decreased as the polygenic variance increased, especially for milk yield. In contrast, the correlations between GEBV and deregressed EBV decreased less when the polygenic variance increased or remained constant, e.g. around 0.72 for somatic cell score. Based on the relatively small decrease in correlations between DRP and DGV or GEBV, we can conclude that the impact of the assumed percentage of residual polygenic variance on accuracy is limited.  from previous generations can identify some possible biases of a genomic evaluation model [24]. The intercept of the linear regression model was not significantly different from zero for all traits. The estimate of the regression slope was nearly unity for the validation population according to the validation procedure [24]. A regression slope estimate that is lower (higher) than its expected value indicates that the variance of the GEBV is too high (too low). According to the regression slope estimates in Table 8, the optimal percentage of residual polygenic variance seems to vary across traits. For traits with a high heritability or reliability, e.g. production traits, somatic cell score, stature and rump angle, the optimal residual polygenic variance appeared to be less than 5%. For the conformation traits, rump width and body conditional score, 10% or higher residual polygenic variances gave the least biased GEBV estimates. Genomic validation results have revealed that either fitting a residual polygenic effect in the BLUP SNP model or blending the G matrix with the pedigree relationship matrix A in the G-matrix BLUP model [13,20,21] was necessary to avoid over-prediction of candidates' GEBV. The optimal proportion of genetic variance assigned to the RPG effect or the optimal weight on matrix A varies across traits. As a result, a trait-specific residual polygenic variance was assumed in routine genomic evaluations for German Holstein cattle. The magnitude of the assumed polygenic variance had a minor effect on the correlation between GEBV and deregressed EBV for selection candidates (Table 7); however, the variance of GEBV decreased significantly with increasing residual polygenic variance. Including the RPG effect in the genomic model (1) provided a similar scale of variances for GEBV and EBV, making them more comparable and consequently resulting in a more accurate joint ranking of genomic selection candidates and proven bulls. However, the problem of optimal partitioning of the additive genetic variance between the residual polygenic and SNP-based components is not resolved. More appropriate statistical methods, such as REML or Bayesian methods [28], should be used to estimate the residual polygenic variance, preferably also including non-genotyped animals.

Influence of the sires' EBV on direct genomic values
A concern that under the genomic BLUP model, animals' DGV are highly correlated with the sires' EBV [25,26] was addressed in this study by fitting an RPG effect with varying residual polygenic variances: 0.02%, 5%, 10% and 20% of the total genetic variance for milk yield. The DGV or the sum of DGV and RPG of 11,978 reference bulls that had genotyped sires in the reference population, were regressed on the conventional EBV of their 580 sires that were also included in the genomic reference population. The corresponding R 2 values indicate the fraction of the sons' genetic variation that is explained by their sires and are shown in Figure 1 for the genomic models with different residual polygenic variances for milk yield. When the RPG effect was given a nearly zero variance i.e. 0.02%, the R 2 value was 0.42  for both DGV and the sum. As the residual polygenic effect increased to 20% of the total genetic variance, the R 2 value between the DGV of the son and the EBV of the sire dropped below 0.20. In contrast to DGV, corresponding R 2 values for the sum of DGV and RPG remained constant, regardless of the level of residual polygenic variance. Figure 2 shows the influence of the sires' EBV on the DGV of validation bulls. The R 2 values of the regression of DGV on the sires' EBV dropped from 0.29 for the scenario with a 0.02% residual polygenic variance to about 0.10 for the scenario with a 20% residual polygenic variance, suggesting a decreasing impact of the sires' EBV on the DGV of validation bulls. With increasing residual polygenic variances, R 2 values decreased much less for combined GEBV of the validation bulls than for DGV alone, because in the combined GEBV the influence of sires was added back via the pedigree index. By fitting an RPG effect in the genomic model, the estimated DGV were less dependent on the sire's EBV, which was indicated by the lower R 2 value of the DGV regression on sire's EBV. The two figures showed that fitting an RPG effect in a genomic model can reduce the correlation between sires' EBV and animals' DGV.

Estimation of SNP effects
Convergence of the BLUP SNP model was improved when the SNP markers were processed in descending order of heterozygosity. The processing order was particularly important when some reference bulls with extremely high or low EBV happened to have extremely high EDC, because those extreme phenotypic values could lead to extreme regression estimates of SNP markers with a low heterozygosity and thus could cause a convergence problem in the estimation of SNP effects. For the currently and most widely used 54 K Illumina Bead-Chip (Illumina Inc., San Diego, CA), we observed that SNP effects did not converge as well as their sum, i.e. DGV. Due to higher LD, convergence of SNP effects could become even lower for a higher density chip, although the convergence of DGV should remain unchanged. An alternative modelling of marker information from high-density chips should be explored.

Conclusions
The tremendous advances in conventional genetic evaluations during the last decades have formed a solid basis for genomic evaluation and selection in dairy cattle. Genomic validation studies worldwide have demonstrated that the genomic model proposed by Meuwissen et al. [10] is highly effective to increase the reliability of evaluations in dairy cattle breeding. In this study, we have shown that the size of the genomic reference population is an important factor affecting the reliability of genomic prediction. Fitting a residual polygenic effect in the genomic model is necessary to avoid the variance of DGV being too high, to make the GEBV of candidates less biased, and to reduce the correlation between reference sires' EBV and animals' DGV. The optimal residual polygenic variance appears to differ between traits. Our validation study has clearly shown that genomic evaluation is efficient.