Selection of haplotype variables from a high-density marker map for genomic prediction

Using haplotype blocks as predictors rather than individual single nucleotide polymorphisms (SNPs) may improve genomic predictions, since haplotypes are in stronger linkage disequilibrium with the quantitative trait loci than are individual SNPs. It has also been hypothesized that an appropriate selection of a subset of haplotype blocks can result in similar or better predictive ability than when using the whole set of haplotype blocks. This study investigated genomic prediction using a set of haplotype blocks that contained the SNPs with large effects estimated from an individual SNP prediction model. We analyzed protein yield, fertility and mastitis of Nordic Holstein cattle, and used high-density markers (about 770k SNPs). To reach an optimum number of haplotype variables for genomic prediction, predictions were performed using subsets of haplotype blocks that contained a range of 1000 to 50 000 main SNPs. The use of haplotype blocks improved the prediction reliabilities, even when selection focused on only a group of haplotype blocks. In this case, the use of haplotype blocks that contained the 20 000 to 50 000 SNPs with the highest effect was sufficient to outperform the model that used all individual SNPs as predictors (up to 1.3 % improvement in prediction reliability for mastitis, compared to individual SNP approach), and the achieved reliabilities were similar to those using all haplotype blocks available in the genome data (from 0.6 % lower to 0.8 % higher reliability). Haplotype blocks used as predictors can improve the reliability of genomic prediction compared to the individual SNP model. Furthermore, the use of a subset of haplotype blocks that contains the main SNP effects from genomic data could be a feasible approach to genomic prediction in dairy cattle, given an increase in density of genotype data available. The predictive ability of the models that use a subset of haplotype blocks was similar to that obtained using either all haplotype blocks or all individual SNPs, with the benefit of having a much lower computational demand.


Background
Since genomic selection methods were introduced [1], single nucleotide polymorphisms (SNPs) are usually used to perform genomic prediction, each as an individual explanatory variable independent from each other. An alternative is to use haplotypes as explanatory variables in genomic prediction. The main benefit of using haplotypes for genomic prediction is that haplotypes are expected to be in higher linkage disequilibrium (LD) with considerations, it is reasonable to assume that haploblocks based on LD should be good explanatory variables for genomic prediction. Assuming a correct marker map, one of the advantages of using haploblocks based on LD is the non-fixed number of SNPs in a haploblock which allowed the grouping of tightly linked adjacent SNPs. Therefore, it efficiently limits the number of "alleles" per haploblock. However, it is clear that the construction of haploblocks based on LD requires an appropriate choice of LD measure and a minimum threshold of LD between markers. The choices that were applied in our work are described in the Methods section. One interesting alternative to LD-based haploblocks are genomic prediction methods based on identity-by-descent (IBD) relationships [5,6]; this alternative benefits from linkage analysis over the genomic data. However, prediction accuracies of breeding values depend strongly on the number of phenotyped and genotyped relatives within the population [6]. Besides, this method using IBD relationships aims at decreasing marker density to reduce genotyping cost, whereas the method based on haploblocks in the current study aims at reducing prediction variables from HD marker data. Thus, the LD-based haploblocks were preferred in this study.
Genomic prediction using a set of appropriately selected haploblocks is expected to achieve higher prediction accuracy while reducing the amount of predictor variables in prediction models. A recent study showed that better predictions in dairy cattle traits can be obtained by using a set of haploblocks with a fixed size (number of SNPs) [7]. Our hypothesis is that using haploblocks that contain the main SNP effects (i.e. the SNPs with the highest absolute effect estimated using the models which estimated effects of all individual SNPs simultaneously, and the haploblocks containing these SNPs are referred to as QTL-haploblocks hereafter) can improve genomic prediction. By using QTL-haploblocks, it is possible to identify the parts of the genome that strongly influence the predictions of livestock traits. In addition, a large proportion of haploblocks may have no effect on a trait, and removing them as covariates may reduce the noise in genomic prediction models. The idea of a QTL-haploblock approach is similar to marker-assisted selection (MAS). However, MAS performs predictions using only a few genetic markers with a significant effect, that was previously estimated by a model including all individual SNPs [8,9], while the proposed QTL-haploblock approach uses genome-wide dense markers and allows a large number of markers into the model, not necessarily all with a significant effect.
Using pre-selected haploblocks for genomic prediction is an important area of research, especially when considering the availability of increasingly denser SNP chips. Reliability of genomic prediction for a trait is expected to be improved by identifying the most influential haploblocks for this trait to be included in the prediction model. In addition, genomic prediction models including a selected group of haploblocks will reduce computing time considerably, compared to models using all haploblocks. This is not necessarily relevant when dealing with moderate-density marker data but plays an important role when dealing with high-density marker data, and is more important when using whole-genome sequence data.
Therefore, this study compared genomic predictions using an individual SNP approach, a haploblock approach with all available haploblocks, and a haploblock approach using a set of haploblocks that contained the main SNPs. The analyses were performed using data from the Nordic Holstein population. The key objective of this work was to assess prediction reliability obtained by using QTLhaploblocks as covariates, and to compare them to those achieved when using all individual SNPs or all haploblocks from a high-density marker chip.

Marker and phenotypic data
The genomic prediction models performed in this study were based on marker and phenotypic data from the Nordic Holstein population. The original marker data was obtained from a 54k SNP chip and then imputed to highdensity (HD) data of 777k SNPs (Ilumina BovineHD Bead-Chip [10]), by applying the Beagle package [11], using 557 HD genotyped reference bulls from the EuroGenomics project [12]. The imputed data was then edited by removing SNPs with a minor allele frequency (MAF) less than 0.01 and also SNPs that were in complete LD with adjacent ones [13]. After editing, the final marker data set included 492 057 SNPs for 5214 bulls.
The phenotypic values to perform genomic prediction for protein yield, fertility and mastitis were pseudoobservations in the form of deregressed proofs (DRP), obtained from the estimated breeding values (EBV) and effective daughter contributions [14][15][16]. All three traits are index traits. The EBV for protein yield and mastitis were composed of EBV for each parity. The EBV for fertility comprised the EBV for interval from calving to first insemination, interval from first to last inseminations, and number of inseminations (heifer and cow separately, pooled over parities for cow). The DRP of all the animals (both training and validation) were derived from the EBV of the official evaluation in August 2010 by the Nordic Cattle Genetic Evaluation (NAV).
To validate the predictive ability of the prediction models, the marker and phenotypic data sets were divided into training and test subsets by the cut-off birth date of bulls on October 1, 2001. The size of training and test data sets are in Table 1, as well as trait reliabilities for these populations.

Animal ethics
The phenotypic data were obtained from routine records of dairy cattle farms. Genotyped animals used in this work were progeny-tested bulls, and the semen samples for genotyping were obtained from routine bull semen collection. Therefore, no ethical approval was necessary.

Genomic prediction covariates
Genomic predictions were performed using individual SNPs and haploblocks. Haploblocks were built based on LD and then selected according to specific criteria. In this section, first we briefly describe the construction of haploblocks and then their selection for genomic prediction.
There are three common pairwise LD measures, D, r 2 and D [17,18]. In this study, D was chosen to define haploblocks according to a previous study [19], and due to the fact that it depends less on individual allele frequencies than D. In addition, a pilot study was performed to compare predictions using haploblocks that were defined using r 2 and D , and no difference on predictive ability was observed. Because the use of r 2 led to many more haploblocks, D was finally chosen as the most adequate LD measure to build haploblocks.
Following our previous study [4], a haploblock was defined as a group of adjacent SNPs such that the LD between any pair of SNPs in this group satisfies |D | ≥ 0.45. This threshold of 0.45 was considered as optimal, considering the prediction reliability, to predict genomic breeding values for the three traits of interest using all the haploblocks built from the HD marker data [4].
Using this LD criterion to define the haploblocks resulted in a total of 76 062 haploblocks. Because haploblocks are "multi-allelic" it summed up a total of 368 709 haploblock variables. The number of SNPs in a haploblock ranged from 1 to 78, with a mean of 6. The number of variants within a haploblock ranged from 1 to 16, with a mean number of 5 [4]. Haploblocks that had only one variant were excluded.
Selection of haploblocks was based on the estimated SNP effects obtained from prediction models using either Bayesian best linear unbiased prediction (Bayesian BLUP) or a Bayesian mixture model, based on the training dataset. Detailed description of the models is provided in the next section, entitled Genomic prediction models. For each trait, the absolute values of the estimated SNP effects were ranked. Then, a determined number k of SNPs with the highest effects was defined. Finally, the haploblocks containing those SNPs were selected to perform genomic prediction. The number k of SNP effects used to select the haploblocks varied from 1000 to 50 000. In the following, the haploblocks selected according to the SNP effects estimated from the training dataset will be referred to as QTL-haploblocks.
Haploblocks that were selected by SNP effects estimated by the Bayesian BLUP model were used for genomic prediction using the Bayesian BLUP model. Analogously, haploblocks that were selected by SNP effects estimated by the Bayesian mixture model were used for genomic prediction using the Bayesian mixture model. Because the estimates of SNP effects differed according to trait and model, the ranking of SNP effects differed as well, thus the number of main SNPs within a haploblock varied. This resulted in different selected haploblocks for each trait, and the number of haploblocks, used to perform genomic prediction.
In order to confirm that genomic prediction using QTLhaploblocks obtains more accurate results than selecting the haploblocks randomly, protein yield was analysed using haploblocks containing 1000 to 50 000 randomly selected SNPs. This procedure was repeated 10 times, and the reliabilities of the predictions were compared to the reliabilities of predictions obtained using QTLhaploblocks.

Genomic prediction models
For the three traits mentioned previously, genomic predictions were performed using a Bayesian BLUP or a Bayesian mixture model, both including the QTL-haploblocks effect and a polygenic effect. The two models used a Bayesian algorithm and were performed using the BayZ package [20], running a single Markov chain Monte Carlo (MCMC) with a length of 50 000, of which the first 20 000 cycles were taken as the burn-in of the chain. Estimates were assessed by the posterior means of the non-discarded 30 000 cycles. Convergence and length of MCMC were monitored by graphical inspection of the dispersion parameter in the models and the correlation between the genomic estimated breeding values (GEBV) from two separate rounds in a pilot study.

Bayesian BLUP model
The Bayesian BLUP model is defined by the equation where y represents the vector containing the DRP of training bulls, μ a general mean, M the haploblock matrix, g the vector of additive haploblock effects, Z the incidence matrix linking a to y, a the vector of residual polygenic additive genetic effects and the vector of random errors of the model. It is assumed that the distributions are as follows, where A is the genetic relationship matrix constructed according to pedigree, D is a diagonal matrix with [16,21]. Furthermore, w i is a weighting factor accounting for heterogeneous residual variances due to differences in r 2 DRPi , the i − th DRP's reliability [22]. The prior uniform distributions were always improper, care was taken to ensure that the overall mean was within the real values and the variances were positive real values.
Taking into account that each haploblock may have more than two variants, matrix M may have more than one column for each haploblock and had dimension n × q (n = number of animals, q = total number of haploblock variables).

Bayesian mixture model
The Bayesian mixture model is defined by the same equation and variables as the Bayesian BLUP model but differs in the assumed distribution of g, the additive haploblock effects, given by This Bayesian mixture model [23] is an extended version of simpler ones [24,25], and intends to facilitate the mixing of the MCMC on the high-density marker data. The mixing proportions π k were fixed at π 1 = 0.889, π 2 = 0.1, π 3 = 0.01 as π 4 = 0.001, and the variances were estimated under the constraint σ 2 π 1 < σ 2 π 2 < σ 2 π 3 < σ 2 π 4 assuming a non-informative prior uniform distribution.

Evaluation of prediction models
GEBV obtained from the prediction models were calculated as GEBV i = j m ijĝj +â i , the performance of each model was assessed by the estimated reliability of GEBV, r 2 and the bias of GEBV. The bias was assessed as b − 1, where b is the regression coefficient b of DRP on the GEBV [22].
The reliability of the prediction for breeding values was obtained as the squared correlation between DRP and GEBV of individuals in the test population corrected for the average reliability of DRP of the test animals r 2 DRP [16]. Thus, the average reliability of GEBV in the test population was calculated as, One of the objectives of this study was to test if fitting a group of selected haploblocks performed as well or better than fitting all haploblocks from the marker data. Thus, reliabilities of models with selected haploblocks were compared to the reliability of the model using all haploblocks with the Hotelling-Williams' t-test [26,27]. Testing whether where r * * refers to the observed values of the correlations ρ * * , as described above, n the number of observations and |R| is the determinant of the correlation matrix R for DRP and GEBV for models i and j. If the null hypothesis is true, then then H 0 is rejected and it is considered that ρ drp,i = ρ drp,j at a significance level α. T statistics and their associated p-values were calculated using R [28]. Table 2 presents the number of QTL-haploblocks selected for each trait and for both statistical models used for genomic prediction. Because the selection of haploblocks was based on the SNP effects obtained from two models that included all individual SNPs, the haploblocks selected differed by trait and model. It can be observed that when using up to 10 000 main SNPs to select QTL-haploblocks, the number of haploblocks did not differ much from the number of main SNPs. Between 20 000 and 50 000 this difference was more pronounced, which means that the first 10 000 SNPs with the highest effects were located in different haploblocks, while thereafter more than one main SNP fell in the same haploblock. Table 3 shows the total number of haploblock variables. Since haploblocks are "multi-allelic", the numbers in Table 3 represent the sum of these alleles, for the selected haploblocks. These were the total number of covariates used in the genomic prediction models. Table 4 presents the prediction reliabilities and bias for the three traits using each prediction model. The row with 492 057 main SNPs corresponds in fact to the haploblock approach using all haploblocks (full haploblocks model) and the last row is the SNP approach. These results were the basis for the comparison of predictions using  cases. An exception was observed for the prediction of fertility using 40 000 main SNPs, for which case the prediction reliability was 0.1 % lower than that achieved by the individual SNP approach. In addition, the differences between the bias obtained by the QTL-haploblock model and the full haploblock model were very small for the three traits. The observed biases, measured as the deviation of the regression coefficients of DRP on GEBV, to 1, was between 0.002 and 0.181 among the three traits. Figures 1, 2 and 3 present the prediction reliabilities in graphs, for protein yield, fertility and mastitis, respectively. These figures show the fast increase of prediction reliabilities when using up to 10 000 main SNP effects to select QTL-haploblocks. Thereafter, the curves stabilize around the reliabilities obtained by the models that used all haploblocks. Figures 4 and 5 compare the results obtained using QTL-haploblocks (blue lines) and randomly selected haploblocks (shaded areas and black lines). Figure 4 shows the results using Bayesian BLUP and Fig. 5 shows the results using the Bayesian mixture model. The random subset of haploblocks was repeated 10 times for each number of SNPs used to select haploblocks (1000 to 50 000), and predictions were performed for each subset. For both Bayesian BLUP and Bayesian mixture models, the mean reliability of the randomly selected haploblocks was lower than those achieved by the QTL-haploblocks, and most of the shaded area is below the blue lines. This confirmed that QTL-haploblocks are better explanatory variables for genomic prediction than haploblocks selected by a random subset of SNPs. It was expected that an advantage of QTL-haploblocks over randomly selected haploblocks would be observed, based on the use of selected individual SNPs for genomic prediction. As shown in Figs. 4 and 5, when a group of individual SNPs were selected based on their estimated effects, the genomic prediction obtained using this group was superior than would be observed if using a randomly selected group of individual SNPs. Table 5 presents P-values of the two-tailed Hotelling-Williams t-test to verify if the reliabilities obtained using QTL-haploblocks were different from those obtained using all haploblocks. The comparisons were made within each trait and each statistical model. Strictly, it is assumed that if the P-value is greater than 0.05, there is no evidence that the prediction reliabilities are statistically different. However, the closer to 1 the P-value is, the stronger the evidence that the prediction reliabilities are not equal. Table 5 shows that the P-values of the test became higher as the number of main SNPs used to select QTLhaploblocks increased. In general, P-values were high (0.632 to 0.999) when using 20 000 to 50 000 main SNPs, and low (0.063 to 0.268) when using 1000 to 3000 main SNPs. One interesting point about the selection of QTLhaploblocks was the relationship between the number of QTL-haploblocks selected by the main SNPs (Table 2) and the number of haploblock variables (Table 3). Figure 6 shows the average number of "alleles" per haploblock, for selection of QTL-haploblocks using 1000 to 10 000 main SNPs. For selection of QTL-haploblocks to predict fertility, the average number of "alleles" per haploblock was greater when using the Bayesian mixture model, than when using Bayesian BLUP. This difference was more accentuated at small numbers of main SNPs (1000 to 2000) to select QTL-haploblocks, and remained almost  unchanged up to 6000 main SNPs. For 7000 main SNPs and more, the average number of "alleles" per haploblock became similar, for both the Bayesian BLUP and the Bayesian mixture models. For mastitis, the number of "alleles" per haploblock is also greater when using the Bayesian mixture model up to 6000 main SNPs, except when using 2000 main SNPs, for which the average was slightly greater when using Bayesian BLUP. This difference is more accentuated when using 3000 to 5000 of main SNPs to select QTL-haploblocks. For 7000 main SNPs and more, the difference started to decrease until the ratios converged to the same value. For protein yield, the scenario observed was different, except when using 1000 main SNPs to select QTL-haploblocks, the number of "alleles" per haploblock was greater when using the Bayesian BLUP model, this was most pronounced up to 4000 main SNPs, then converged to the same value. Table 6 indicates the computing time required to perform the prediction with each model. It was observed that computing time increased according to the increase of main SNPs used to select QTL-haploblocks, hence the increase in predictive variables. It was clear that there was a reduction in computing time, when QTL-haploblocks  Prediction reliabilities obtained using Bayesian BLUP model with QTL-haploblocks and random selected haploblocks as covariates for protein yield. The values on the x-axis are the number of main SNPs used to select QTL-haploblocks to perform genomic prediction, and the y-axis indicates the reliability of predictions. The grey shaded area shows the range (minimum to maximum prediction reliabilities) and the black lines indicate the mean reliabilities obtained of the models using the randomly selected haploblocks were used to perform genomic prediction, and this is one relevant benefit of this method.

Discussion
Previous studies have already determined that haploblocks are able to better predict breeding values of economically important traits in dairy cattle, than individual SNPs [3,4,7]. Similar to results obtained in [7], the QTL-haploblocks used as predictors in this work may achieve predictions that are more accurate than those using all 492 057 individual SNPs, and as good as those achieved using all 76 062 haploblocks built from the genomic data.
When compared to the predictions with individual SNPs, the QTL-haploblocks were able to achieve higher prediction reliabilities for the three traits, when selected   Table 4 and clearly observed in Figs. 1, 2 and 3. This range of 20 000 to 50 000 main SNPs to select QTL-haploblocks was also found to result in prediction reliabilities either equal or very close to those obtained using all haploblocks for protein yield and fertility. In the prediction of mastitis, the reliabilities observed in this range were equal or even greater than those obtained using all haploblocks (up to 0.8 %). This particular result for mastitis was satisfying, taking into account that it is a trait with low heritability, and it is difficult to improve its prediction accuracy using haploblocks [4]. Furthermore, genetic progress is linearly related to accuracy of genetic evaluation. Considering a large dairy cattle population, a small improvement in reliability of predictions is considered important for cattle breeding.
The P-values of the Hotteling-Williams tests in Table 5 were used to compare the results of using only QTLhaploblocks in prediction reliability, with that when using all haploblocks. There was a strong interest in verifying how the increase of variables in the prediction model affected the evidence (P-value) favouring the hypothesis of equal reliabilities. For all the traits, the P-values indicated that the prediction reliabilities using QTLhaploblocks selected by at least 4000 main SNPs (i.e. the SNPs with the highest effects) were statistically not different to those obtained using all haploblocks. We could observe, moreover, that while we used 1000 to 10 000 main SNPs to select QTL-haploblocks, the Pvalues comparing the prediction reliabilities to the full haploblocks model increased regularly. This means that for up to 20 000 main SNPs, the more SNPs we use to select QTL-haploblocks, the stronger becomes our evidence that the prediction reliabilities are not different. In the range of 20 000 to 50 000 main SNPs used to select QTL-haploblocks, the P-values were high (greater than 0.6), suggesting that those models predict as well as or equally well as the model using all haploblocks. Fig. 6 Ratio between total number of haploblock variables used in the prediction models and total number of haploblocks containing the main SNP effects. The values on the x-axis are the number of main SNPs used to select QTL-haploblocks to perform genomic prediction, and the y-axis indicates the ratio (haploblock variables)/haploblocks One important feature of prediction models using QTLhaploblocks is the reduction in computing time. Compared to the individual SNP approach, the models with QTL-haploblocks take approximately only 20 to 25 % of the computing time, and 30 to 41 % when compared to using all haploblocks. Thus, although the increase in prediction reliability is small compared to the individual SNP approach (however still very important in cattle breeding), the increase in computational efficiency was considerable. Furthermore, in all our models, we used a MCMC with a fixed length of 50 000 iterations, and the first 20 000 were discarded as burn-in. Because the models that use QTL-haploblocks have significantly less explanatory variables, and because these variables are also less correlated to each other as are individual SNPs, it is expected that the number of MCMC iterations can be reduced. Consequently, a further reduction in computing time can be achieved.
For low or moderate density marker data, the computational gain provided by the use of QTL-haploblocks, from preparation of data and time to run prediction models, may not be relevant. However, the predictions using QTLhaploblocks will be of great importance when it comes to denser marker data, such as whole-genome sequence data. Hence, further studies on genomic prediction using haploblocks and QTL-haploblocks based on LD is a natural next step to evaluate the benefits from these predictors.

Conclusions
The results from this study suggest that when 20 000 to 50 000 main SNPs were used to select QTL-haploblocks, the use of QTL-haploblocks as predictors is generally sufficient to obtain reliabilities equal or higher than those obtained using all individual SNPs (up to 0.9 % increase for proteinyield, equivalent prediction for fertility and up to 1.3 % increase for mastitis, compared to the individual SNP approach), and similar to those obtained using all haploblocks.
In addition, the method presented here had a positive effect on computing time for prediction models using HD marker data. Compared to the computing time required for the models using all haploblocks, the model using the QTL-haploblocks containing 20 000 to 50 000 main SNPs took on average 40 % of the total time needed and obtained statistically similar results. The computing time for the models using QTL-haploblocks can be further reduced by using less MCMC cycles, since there are less explanatory variables. With denser marker data and whole-genome sequence data, the reduction in computing time would be an important issue in practical genomic prediction.