Genomic prediction from observed and imputed high-density ovine genotypes

Background Genomic prediction using high-density (HD) marker genotypes is expected to lead to higher prediction accuracy, particularly for more heterogeneous multi-breed and crossbred populations such as those in sheep and beef cattle, due to providing stronger linkage disequilibrium between single nucleotide polymorphisms and quantitative trait loci controlling a trait. The objective of this study was to evaluate a possible improvement in genomic prediction accuracy of production traits in Australian sheep breeds based on HD genotypes (600k, both observed and imputed) compared to prediction based on 50k marker genotypes. In particular, we compared improvement in prediction accuracy of animals that are more distantly related to the reference population and across sheep breeds. Methods Genomic best linear unbiased prediction (GBLUP) and a Bayesian approach (BayesR) were used as prediction methods using whole or subsets of a large multi-breed/crossbred sheep reference set. Empirical prediction accuracy was evaluated for purebred Merino, Border Leicester, Poll Dorset and White Suffolk sire breeds according to the Pearson correlation coefficient between genomic estimated breeding values and breeding values estimated based on a progeny test in a separate dataset. Results Results showed a small absolute improvement (0.0 to 8.0% and on average 2.2% across all traits) in prediction accuracy of purebred animals from HD genotypes when prediction was based on the whole dataset. Greater improvement in prediction accuracy (1.0 to 12.0% and on average 5.2%) was observed for animals that were genetically lowly related to the reference set while it ranged from 0.0 to 5.0% for across-breed prediction. On average, no significant advantage was observed with BayesR compared to GBLUP.


Background
The development of high-throughput genotyping based on single nucleotide polymorphisms (SNPs) in livestock species has made the implementation of genomic evaluation more practical. In genomic prediction, the breeding values of selection candidates are evaluated according to their genotypes and a prediction equation derived from a reference population with both phenotypes and genotypes [1]. The accuracy of genomic prediction relies on several factors including linkage disequilibrium (LD) between genome-wide SNPs and quantitative trait loci (QTL) that are responsible for the phenotypic variation of traits of interest [1]. High-density (HD) SNP genotypes can result in stronger LD between SNPs and QTL which can improve the accuracy of genomic prediction in livestock, e.g. [2][3][4][5].
Results of simulation studies in livestock show various degrees of improvement in genomic prediction when using HD genotypes compared to genotypes from moderate-density SNP panels such as 50k. For example, based on simulation studies, Meuwissen and Goddard [6] reported a large gain (>40%) in prediction accuracy from HD genotypes, while VanRaden et al. [7] and Harris and Johnson [8] found zero to only small gains in prediction accuracy. Such differences can be attributed to the assumption made about the distribution of QTL effects in the simulated models. Meuwissen and Goddard [6] and Clark et al. [9] showed that both the number and distribution of QTL effects that control a polygenic trait have a significant impact on the advantage of using HD genotypes in genomic prediction, with only small benefits for the "infinitesimal' model for which most of the variation of a trait is due to a large number of QTL each with a relatively small effect.
Analyses of real data are available from dairy cattle and show zero to relatively small increases in prediction accuracy from HD genotypes. Solberg et al. [10] reported between 0.0 and 9.0% improvement in prediction accuracy across seven production and functional traits in Norwegian Red bulls. VanRaden et al. [7] found up to 6.5% (on average 0.4%) extra accuracy across 28 production traits using HD genotypes in Holstein dairy cattle.
Initially, the first factor that was suggested to affect the accuracy of genomic prediction was the LD between genome-wide SNPs [1,2]. However, it was later shown that genomic prediction accuracy depends both on cosegregation of SNP alleles in related individuals and information from SNP alleles being in LD with QTL alleles e.g. [11]. Prediction accuracy based on LD is more persistent over distant relationships and the expectation is that higher density SNP arrays are better at capturing effects of QTL that are in LD with SNPs. Therefore, the advantage of using HD genotypes is expected to be greater for animals that are less genetically related to the reference set, and this could apply to both within-breed and across-breed genomic prediction. Thus, denser SNP genotypes may have a favorable effect on the accuracy of genomic prediction in multi-breed and crossbred populations, which are common in the sheep and beef cattle industries. Harris et al. [12] and Erbe et al. [13] showed that there was very limited improvement from using HD genotypes in across-breed prediction in Holstein and Jersey dairy cattle, but differences may be larger in sheep where breeds are genetically more related to each other and have a larger effective population size.
The objective of this study was to compare the accuracy of genomic prediction for weight, ultra-sound scanned fat and muscle traits, and wool quality and quantity traits in Australian sheep breeds based on both observed and imputed HD genotypes (600k Illumina Ovine SNP) to accuracies based on moderate-density SNP genotypes (Illumina ovine SNP50k). Using a reference set comprised of purebred, crossbred or mixed crossbred and purebred animals, prediction accuracies were compared for purebred industry sires for which very accurate estimated breeding values based on a progeny test were available. Furthermore, we contrasted accuracy of genomic prediction within a breed between animals with low and high genetic relatedness to the reference set as well as prediction within and across breeds.

Reference set, phenotypes and validation population
The genomic prediction reference set consisted of about 20,000 animals that were recorded for a large number of production traits measured in the "Sheep Cooperative Research Centre Information Nucleus Flock" (INF) and "Sheep Genomics Flock" (SGF). The INF consisted of eight flocks that are located across different regions of Australia and are linked to each other because artificial insemination with common sires was used between 2007 and 2011 [14]. The SGF was a single research flock located in southern New South Wales, Australia, for which data were collected between 2005 and 2006 [15]. All animals in the reference set were from multiple breeds or crossbreds with the sires comprising approximately 40% animals from Terminal breeds [Poll Dorset (PD) and White Suffolk (WS)], 20% from a Maternal breed [Border Leicester (BL)] and 40% from Merino and the dams comprising 80% Merino and 20% BL × Merino crossbreds. The dominant purebred animals were Merinos which included three sheep strains that have different wool qualities, i.e. strong wool, fine wool and ultra-fine wool types. The traits analyzed were live body weights from birth to adult age, ultra-sound scanned muscle and fat depth measured at post-weaning age and wool quantity and quality measured at yearling and adult age. The data used in this study was collected according to the guidelines of the "University of New England Animal Ethics committee" reference number AEC 09/115. The number of records and basic statistics per trait are summarized in Table 1.
A validation population was used to find the empirical accuracy of genomic prediction. The validation population was a group of industry purebred sires with accurate estimated breeding values (EBV) (accuracy ranging from 0.70 to 0.99 and on average 0.92), which were calculated based on progeny records. The phenotypes of INF and SGF animals (genomic prediction reference set) were not used in the calculation of EBV of the validation sires.

Genotypes
The reference and validation populations were genotyped using a 50k SNP panel (Illumina Inc., San Diego, CA, USA). This 50k SNP panel provided 48,559 SNP genotypes after applying quality control based on the following criteria: individual SNP genotypes were removed if their call rates were lower than 90%, or if the GenCal (GC) scores were <0.6, if the heterozygosity rate for a given SNP deviated more than 3 SD from the population mean, if the SNP minor allele frequency was lower than 0.01, and for SNPs located on chromosomes X and Y or SNPs that deviated from Hardy-Weinberg equilibrium (P < 1 × 10 −15 ). Furthermore, an individual sample was removed if the correlation of its genotypes (coded 0, 1 or 2 per locus) with those of another sample was equal or greater than 0.98.
Most of the sires and 1735 progeny from the four main breeds including Merino, BL, PD and WS were genotyped using the HD (Illumina Inc., San Diego, CA, USA) ovine SNP panel. This SNP panel provided 510,174 SNPs after applying the same quality controls as above. Using all HD genotyped animals as imputation reference set, the un-typed genotypes of the rest of the population were imputed to HD genotypes using the software program FImpute [16]. The accuracy of imputation, which was tested within subsets of animals with observed HD genotypes, was high (on average 0.98).

Statistical methods
For the analysis based on pedigree relationships, the following mixed model was fitted using ASReml 3.0 [17]: where y is a vector of phenotypes, b is a vector of fixed effects, a is a vector of random additive polygenic effects, w is a vector of random maternal effects, q is a vector of random breed effects, s is a vector with random sire y = Xb + Z 1 a + Ww + Z 1 Qq + Z 2 s + e, by flock interaction effects and e is a vector of random residual effects. X, Z 1 and W and Z 2 are incidence matrices relating fixed effect, additive genetic, maternal effects and sire by flock interaction effects to phenotypes. Q is a matrix with breed proportions for each animal derived from pedigree data. Up to 28 breed effects, including those of the three Merino strains, were estimated via the Q matrix, however the major breeds were Merino, BL, PD and WS. All random effects are identically and independently distributed except for a which is distributed as: a ∼ N 0, Aσ 2 a , where A is a numerator relationship and σ 2 a is the additive genetic variance. The fixed effects in the model were birth type, rearing type, gender, age at measurement, weight at measurement and contemporary group which was defined as a cohort of site × birth year × management group. The model used for the estimation of variance components and prediction of genomic breeding values (GBV) was the same except that A was replaced by G, where G is a genomic relationship matrix calculated based on 50k or HD SNP genotypes using VanRaden's [18]

equation as below:
where M is a matrix of the size n × m (i.e. number of individual by number of SNPs) with coefficients equal to (2 − 2p j ), (1 − 2p j ) and (−2p j ) for genotype (A 1 A 1 ), (A 1 A 2 ) and (A 2 A 2 ) of the jth SNP genotype respectively, p j is the frequency of allele A 1 for the jth SNP genotype. σ 2 g is the additive genetic variance estimated from SNPs. Variance components were estimated according to the restricted maximum likelihood (REML) method using either pedigree information or genomic information from 50k or HD genotypes. Genomic EBV (GEBV) were also calculated based on a Bayesian method (BayesR [13]) in which BESSiE [19] was used for prediction of GBV based on the following model: where m refers to the random effects of SNPs, M 1 is an incidence matrix relating SNP effects to phenotypes and the other terms are the same as described above. A mixture of four normal distributions for SNP effects with variances σ 2 1 = 0, σ 2 2 = 0.0001σ 2 g , σ 2 3 = 0.001σ 2 g , and σ 2 4 = 0.01σ 2 g was considered in BayesR where σ 2 g is the assumed total genetic variance. The starting values for σ 2 g were taken from GREML analysis and the prior distribution of the proportion of SNPs in each distribution was the Dirichlet distribution. A total of 50,000 iterations (with 10,000 burn-in) were run for analysis.
The accuracy of GBV was assessed in a separate population of purebred industry rams including Merino, Maternal and Terminal sires (validation set), as the   Pearson correlation coefficient between GBV and an accurate EBV estimated from progeny test. Correlations were estimated for each breed separately, while an effect due to the Merino strain was fitted to avoid GBV accuracy to be biased upward for merinos by evaluating accuracy across strains. The size of the validation set for different traits was 341 to 389 sires for Merino, 79 to 88 for BL, 161 to 188 for PD and 189 to 204 for WS. We also contrasted the accuracy of GBV for animals with high or low genomic relationships with the reference set. Animals with high genomic relatedness were those for which the average value of their 30 highest genomic relationships to the reference population was at least 0.20. Animals with low genetic relatedness were those for which the genomic relationship with any of the individuals in the reference set was not higher than 0.10. Table 2 shows the genetic and residual variance components of the studied traits as well as the estimated heritability based on the genetic covariance matrix among animals that was estimated from pedigree or marker genotypes (50k or HD). Additive genetic variances and heritability estimates based on 50k SNP genotypes tended to be lower than those based on pedigree data (heritability was on average 4.9% lower across different traits). Other variance components including the maternal effect and the sire by site (genotype by environment) interaction effects varied little between different models and are not reported in Table 2. In most cases, estimated residual variances were slightly larger from a model based on 50k genotypes compared with those based on pedigree relationships. Variance components estimated by using HD genotypes resulted in larger additive genetic variance, smaller residual variance and hence higher heritability across all studied traits, when compared to 50k genotypes. However, the increase in additive variance and heritability was small (up to 4% of the absolute value for heritability). Variance components and heritability estimates were similar between models that used HD genotypes and pedigree. Less than 1% differences were found between heritability estimates based on HD genotypes and pedigree when averaged across all weight, carcass scan and wool traits.

Table 2 Additive (V A ) and residual (V R ) variance components and heritability estimate based on pedigree (PBLUP) and 50k (GBLUP-50k) or HD SNP genotypes (GBLUP-HD)
a Standard error of heritability was between 0.02 and 0.09; for trait abbreviations see Table 1 Trait PBLUP GBLUP-50k GBLUP-HD  Tables 3, 4 and 5 show the accuracy of genomic prediction for weight and scanned carcass traits for Merino, BL, PD and WS sires, based on GBLUP (both for 50k and HD SNP genotypes) and BayesR and using the complete multi-breed reference set. Compared to 50k SNP genotypes, the HD SNP genotypes provided higher prediction accuracy but the extra accuracy was on average small. The maximum improvement in prediction accuracy as absolute value was 7.7% and was on average equal to 1.6, 1.2, 4.3 and 3.1% for Merino, BL, PD and WS sires, respectively. Terminal breeds showed a higher increase in prediction accuracy (3.7%) compared to Merino and Maternal breeds (1.4%), which suggests a tendency for greater improvement in accuracy from HD genotypes for breeds with a lower overall accuracy. When using HD genotypes, the accuracy of genomic prediction was very similar between GBLUP and BayesR across all traits, with an average absolute value of the difference in genomic prediction accuracy between GBLUP-HD and BayesR of −0.008, −0.006 and 0.03 for Merino, Maternal and Terminal breeds, respectively. Table 6 shows the accuracy of genomic prediction of breeding value for wool traits in Merino sires based on GBLUP-with 50k and HD SNP density, and BayesR using HD SNP density with only Merinos in the reference set. The extra accuracy resulting from HD genotypes ranged from 0.0 to 8.0% with an average of 5.0%. No considerable difference in accuracy was observed between GBLUP and BayesR. Table 7 shows the accuracy of genomic prediction within and across breeds for three weight traits and two Table 3  where r is the correlation coefficient and n is the number of paired observations; for trait abbreviations see Table 1 Trait

Genomic prediction within and across breeds from purebred or crossbred reference sets
Size GBLUP-50k where r is the correlation coefficient and n is the number of paired observations; for trait abbreviations see Table 1 Trait

Merino
Size

Table 5 Accuracy of genomic prediction of weight and scanned traits for Merino, Border Leicester (BL), Poll Dorset (PD) and White Suffolk (WS) sires based on the multibreed reference set and BayesR based on HD genotypes
a Standard error (SE) calculated according to: 1−r 2 n−2

0.5
where r is the correlation coefficient and n is the number of paired observations; for trait abbreviations see Table 1 Trait Size BayesR scanned carcass traits. Using HD genotypes and a purebred Merino reference set resulted in a small increase in GBV accuracy (0.0 to 2.5%) for Merino sires, which was similar to the increase in genomic prediction accuracy in Tables 3, 4 and 5. A larger increase (0.3 to 9.6%) was observed for Merino sires based on prediction from crossbred Merinos. However, it should be noted that the magnitude of the prediction accuracy for Merino sires from crossbred Merinos is still much lower than the prediction from purebred Merinos. The data in Table 7 can be used to infer the accuracy of genomic prediction across breeds. The increase in genomic prediction accuracy for BL, PD or WS sires from a purebred Merino reference set, which is genetically distant to the target breeds, was low and showed a small non-significant improvement in prediction accuracy when moving from 50k to HD prediction. However, genomic prediction of PD and WS sires based on a combined crossbred reference set (PD × M + WS × M) where r is the correlation coefficient and n is the number of paired observations; for trait abbreviations see Table 1 Trait

Merino
Size GBLUP-50k GBLUP-HD BayesR  where r is the correlation coefficient and n is the number of paired observations; for trait abbreviations see Table 1 Trait Reference set Size GBV accuracy (50k) GBV accuracy (HD) showed a greater improvement in prediction accuracy (up to 8.0%). It should be noted that this accuracy was still low, even when using HD genotypes. Figures 1 and 2 compare the accuracy of genomic prediction for two groups of Merino sires used as validation animals, one with a high and one with a low genomic relationship to the purebred Merino reference set. For highly related animals, the gain in accuracy from using HD genotypes was very low (on average 0.8%) but it was significantly higher for lowly related animals (up to 12% and on average 5.2%).

Discussion
This study investigates the possible improvement in accuracy of genomic prediction of breeding values for weight, scanned carcass and wool quantity and quality traits in Australian sheep when using high-density SNP genotypes. First, we compared the variance components that were estimated based on relationships derived from 50k and HD genotypes to those based on pedigree relationships. Estimated additive genetic variances based on HD genotypes were larger than those based on the 50k SNP panel, which suggests that the HD panel captures more genetic variation; this is likely due to higher LD between SNPs and QTL. Estimated genetic variances based on the HD panel were similar while the estimates based on the 50k panel were lower than those based on pedigree data. However, the A and G matrices are not necessarily on the same scale (e.g. the G matrix is derived as a genomic relationship) so these estimates cannot be directly compared. Haile-Mariam et al. [20] also reported that the additive genetic variances and heritabilities estimated from Bovine50k genotypes were lower than those based on pedigree BLUP for 29 production traits in Australian dairy cattle. Legarra [25] argued that the relationship matrices used to estimate genetic variances should be comparable, i.e. the same average relationship and the same average inbreeding. In any case, the difference between 50k and HD panels is the most relevant comparison and this is not affected by scaling.  HD SNP panels provided higher prediction accuracies but the increase had only practical significance for individuals that were not closely related to the reference population. The average improvement in prediction accuracy was small, ~2.2% which is likely due to the effect of closer relationships providing information that is not much improved by higher marker density. SNPs can capture cosegregation of alleles (family relationships) as well as the LD between SNPs and QTL [5,11,21,22]. Co-segregation is based on linkage between SNPs and QTL which exists over much larger chromosomal regions, therefore not requiring a very high SNP density for adequate prediction. Van der Werf et al. [22] pointed out that prediction from closer relatives is similar to prediction in populations with a lower effective size in which fewer effective chromosome segments are segregating. This observation leads to the same conclusion, i.e. that higher SNP density will have a limited effect on the prediction accuracy when the relationship between reference and target set is stronger.
Previous reports based on real data in dairy cattle also showed a very limited improvement in prediction accuracy when using HD genotypes [7,8], which confirm results from some simulation studies [12,23]. However, Meuwissen and Goddard [6] showed a larger gain in prediction accuracy, using a simulation model that included more QTL with large effects, e.g. all the genetic variation of a polygenic trait was due to three to 30 QTL segregating on one chromosome. Meuwissen and Goddard [6] and Clark et al. [9] showed that the use of denser SNP panels was more beneficial if traits are controlled by fewer QTL with larger effects. Our results show limited extra accuracy from HD genotypes, which could indicate that the distribution of QTL effects is closer to the infinitesimal model assumption.
Genomic prediction in a multi-breed reference set could potentially benefit from across-breed prediction when using HD genotypes, as has been suggested in various studies [24,26,27]. However, we observed only a small (from 0 to a slightly positive value) increase in accuracy when using information from other breeds. Acrossbreed prediction could be lower due to differences in both QTL and SNP allele frequencies, incomplete LD between SNPs and QTL across breeds and different allele substitution effects at QTL in different breeds, e.g. due to epistatic interactions [28]. Using higher density SNPs would address only the incomplete LD aspect but not the other two factors. In this study, a slightly greater improvement in GBV accuracy from using HD genotypes was observed for purebred Merinos (5%) based on a Merino reference set compared to a larger multi-breed reference set. Very limited prediction accuracy from HD genotypes was found for PD and WS breeds based on the Merino sheep reference set, which is likely due to the large genetic distance between Merino and PD or WS as terminal breeds. These results are in line with those of other across-breed prediction studies, e.g. [12,27] who reported small to no across-breed prediction accuracy from a combined Holstein and Jersey dairy cattle reference set. Interestingly, our results showed a notable (on average 5%) improvement in genomic prediction of PD or WS sheep based on a combined crossbred PD or WS reference set. This suggests that HD SNP panels could be useful to improve LD between SNPs and QTL within diverse breeds or between closely-related breeds, in which case it is also more likely that QTL effects are similar. However, predictions across more distant breeds will not benefit from HD genotypes due to lower levels of LD and possibly larger differences in QTL effects.
Some studies have shown that using moderate-density SNP panels (~50k) provide a more marked improvement in genomic prediction accuracy over low-density SNP panels in different livestock species. Moghaddar et al. [29] compared prediction based on panels of 5k, 10k, 20k and 50k SNPs and showed on average a 11 to 13% gain in prediction accuracy for different production traits in Merino sheep. In dairy cattle, Moser et al. [30] reported on average 10% extra accuracy by switching from very low-density SNP genotypes (3000 to 5000) to moderatedensity SNP genotypes (50k). Other studies have also reported relatively large improvements in prediction accuracy from using moderate-density SNP panels compared to low-density SNP sets [3,5,31]. However, this study showed improvements in prediction accuracy from using ovine HD genotypes compared to moderate-density genotypes (ovine 50k) seems generally much smaller, but significant improvements were still observed for individuals distantly related to the reference population. This is consistent with the theory about genomic prediction accuracy [32].
The regression coefficient of EBV on GBV was on average higher (less biased) based on HD SNPs than on 50k SNPs. This could be related to the larger additive genetic variances that were estimated when using HD genotypes and are more similar to the estimates of additive genetic variance based on pedigree data. Bias could also occur if selected SNPs were used for genomic prediction. To some extent, the BayesR method uses selected SNPs, in the sense that it uses some priors to emphasize a larger effect for some SNPs by giving them more weight. However, regression coefficients did not differ between GBV based on GBLUP using HD genotypes and GBV based on the BayesR method, which suggests that this explanation is less likely.
Regression coefficients of EBV on GBV were generally lower than 1.00 (0.74 to 0.94). This may be due to the G-matrix not being expressed at the same scale as the numeric relationship matrix (A) used in the genetic evaluation that produces the EBV, or because of differences in the method for accounting for genetic groups in the reference and validation populations. The A-matrix is based on pedigree relationships whereas GBV are calculated with a G-matrix that uses relationships across various subpopulations within the population. Since this study was mainly aimed at evaluating genomic prediction accuracy, we did not attempt to rescale the G-matrix, since accuracy is calculated as a correlation which is independent of scale. Furthermore, the averages of diagonal and off-diagonal elements of A and G were similar (1.01 and 0.00 for A, 1.00 and 0.00 for G based on 50k SNP density and 1.03 and 0.00 for G based on HD density) as was suggested by Legarra [25] as a requirement to obtain unbiased estimation of breeding values.

Conclusions
Our results show that the use of high-density (600k) SNP genotypes for the genomic prediction of weight and wool production traits in a multi-breed sheep population resulted in a small improvement in accuracy compared to a moderate SNP density (50k). Improvement in accuracy was greater for individuals that were distantly related to the reference set. Prediction accuracy based on a reference set from other breeds was low and showed limited improvement with HD genotypes. Results of GBLUP and BayesR were not significantly different.