Multi-trait single-step genomic prediction accounting for heterogeneous (co)variances over the genome

Widely used genomic prediction models may not properly account for heterogeneous (co)variance structure across the genome. Models such as BayesA and BayesB assume locus-specific variance, which are highly influenced by the prior for (co)variance of single nucleotide polymorphism (SNP) effect, regardless of the size of data. Models such as BayesC or GBLUP assume a common (co)variance for a proportion (BayesC) or all (GBLUP) of the SNP effects. In this study, we propose a multi-trait Bayesian whole genome regression method (BayesN0), which is based on grouping a number of predefined SNPs to account for heterogeneous (co)variance structure across the genome. This model was also implemented in single-step Bayesian regression (ssBayesN0). For practical implementation, we considered multi-trait single-step SNPBLUP models, using (co)variance estimates from BayesN0 or ssBayesN0. Genotype data were simulated using haplotypes on first five chromosomes of 2200 Danish Holstein cattle, and phenotypes were simulated for two traits with heritabilities 0.1 or 0.4, assuming 200 quantitative trait loci (QTL). We compared prediction accuracy from different prediction models and different region sizes (one SNP, 100 SNPs, one chromosome or whole genome). In general, highest accuracies were obtained when 100 adjacent SNPs were grouped together. The ssBayesN0 improved accuracies over BayesN0, and using (co)variance estimates from ssBayesN0 generally yielded higher accuracies than using (co)variance estimates from BayesN0, for the 100 SNPs region size. Our results suggest that it could be a good strategy to estimate (co)variance components from ssBayesN0, and then to use those estimates in genomic prediction using multi-trait single-step SNPBLUP, in routine genomic evaluations.


Background
Genomic selection was pioneered by the study of Meuwissen et al. (2001), and is rapidly becoming the state-ofthe-art genetic selection methodology in many breeding programs around the world. The models proposed by Meuwissen et al. (2001) include a BLUP model, where the variances of single nucleotide polymorphism (SNP) effects are assumed to be the same for all SNPs (SNPBLUP), or specific to each SNP (BayesA and BayesB). Under a series of assumptions, the SNPBLUP model is equivalent to a mixed linear model, GBLUP (Habier et al. 2007), which uses a relationship matrix (G) computed from genetic markers (Nejati-Javaremi et al. 1997) to model covariances between individuals' genetic effects (Stranden and Garrick 2009). This equivalency resulted in a widespread adoption of genomic prediction in genetic evaluations, because only an extra step of computation of G and its inverse is required for the traditional mixed model equations (Henderson 1984) used in animal breeding (Karaman et al. 2016). Moreover, it also allows all extensions of BLUP methodology, such as multiple-trait, random regression, or repeated measures to be easily implemented in genomic evaluations (Tiezzi and Maltecca 2015). The GBLUP model has been widely used to predict breeding values in animal species, such as cattle (Luan et al. 2009;Su et al. 2012b), pig (Lukić et al. 2015), sheep (Daetwyler et al. 2010a) and fish (Ødegård et al. 2014;Tsai et al. 2016), and accuracies from GBLUP were reported to be higher than those from traditional pedigreebased BLUP.
Although widely used in genomic evaluations, these BLUP-based genomic prediction models have some drawbacks. First, they ignore the fact that a large proportion of the SNPs may not have any influence on the trait of interest. Second, different loci or genomic regions may have rather different variances. The two models of Meuwissen et al. (2001), BayesA and BayesB, were proposed to overcome such drawbacks. Assuming SNP-specific variances, BayesA fits each of SNPs, while BayesB fits approximately 1-π of the SNPs, where π is the percentage of SNPs which have no influence on the trait of interest. When π = 0, BayesB is equivalent to BayesA. As pointed out by Gianola et al. (2009), both models are problematic as full conditional posteriors of the SNP-specific variances have only one additional degree of freedom compared to their priors regardless of the amount of data available. A simpler model that similarly fits approximately 1-π of the SNPs, but with a common variance, BayesC, was also proposed (Meuwissen 2009;Kizilkaya et al. 2010). Zeng et al. (2016) introduced a Bayesian partitioned regression model for genomic prediction, which involves the selection of genome regions followed by the selection of SNPs within those selected regions. The model fits approximately 1 − Π of the regions assuming regionspecific variances, and 1 − π s of the SNPs within the region s assuming a common variance for the SNPs in the region. Referring to this "nested" variable selection structure of the model, it was termed as BayesN. The special case of the partitioned regression model of Zeng et al. (2016), i.e., BayesN with Π = π s = 0 (hereafter, BayesN0), is equivalent to BayesA or GBLUP when a fixed region size is set at one SNP or the whole genome, respectively. We hypothesize that, at any other region size, but these two extreme sizes of genome regions, higher prediction accuracies can be obtained using BayesN0. Although it ignores the fact that a proportion of the genome regions, and therefore a proportion of the SNPs, may not have any influence on the trait of interest, prediction accuracy may increase compared to BayesA by benefiting from the increase in the accuracy in estimation of SNP variances, and compared to BLUP-based models by allowing SNPs in different regions to have different variances. Partitioning of the covariate matrix of marker genotypes, M, or in other words, assigning priors to genome regions rather than individual SNPs, was shown to influence the accuracy of genomic predictions (Brøndum et al. 2012;Gebreyesus et al. 2017;Karaman et al. 2018).
Many important traits in animal breeding have genetic correlations in varying sizes with one or more traits, and therefore, measurements of such correlated traits carry information for the genetic values of others. Several multitrait models have been proposed for genomic prediction (Calus and Veerkamp 2011;Jia and Jannink 2012;Hayashi and Iwata 2013;Gebreyesus et al. 2017;Cheng et al. 2018b), and simulations have shown that genomic prediction accuracies from multi-trait models are superior to those from single-trait models (Calus and Veerkamp 2011;Jia and Jannink 2012;Guo et al. 2014;Karaman et al. 2018). Multitrait genetic evaluation rely on the genetic association between the traits through the genetic variance and covariance structure. Models used for genomic prediction, therefore, should properly account for the makeup of these genetic (co)variance components to obtain the highest accuracy of prediction. When only a few genome regions explain a considerable amount of the variances and/or covariance in a two-trait analysis, models that account for the heterogeneous correlation structure over the genome may have advantages over the methods that assumes a constant correlation over the genome (Gebreyesus et al. 2017;Karaman et al. 2018).
The GBLUP model was extended to utilize all phenotypic, pedigree and genotypic information simultaneously, including phenotypic information on non-genotyped individuals, and termed as single-step GBLUP (ssGBLUP) (Christensen and Lund 2010;Aguilar et al. 2010). In ssGBLUP, the pedigree-based relationship matrix A and the genomic relationship matrix G are combined into a single matrix H. As for GBLUP, only an extra step for computation of H and its inverse is required for the traditional mixed model equations used in animal breeding . However, ssGBLUP also suffers from the same drawbacks of GBLUP. Fernando et al. (2014) proposed a class of single-step models, which not only unifies all available information as ssGBLUP does, but also accommodates any Bayesian whole genome regression model. This yields models of, for instance, ssBayesA or ssBayesN0, referring to the Bayesian whole genome regression model used in the single-step analysis. However, such an approach requires that all unknowns of the model to be estimated using Markov-chain Monte Carlo techniques which may be computationally infeasible especially in routine genomic evaluations. In genomic predictions using weighted GBLUP, it was shown that the use of the same SNP variances over a few years does not reduce prediction accuracy . Indeed, in routine evaluations, variance components are not updated for each round of evaluation, because they are expected to be relatively consistent over time (Calus et al. 2014). An alternative to the fully Bayesian approach in Fernando et al. (2014) could be a strategy, where all necessary parameters are estimated using a Bayesian whole genome regression model first, and mixed model equations are then solved given the "known" values of the variance components, leading to a single-step SNPBLUP (ssSNPBLUP) model.
The aim of this study was three-fold: (i) to introduce a multi-trait whole genome regression model that allows heterogeneous (co)variances, (ii) to compare accuracies from single-and multi-trait genomic prediction, and (iii) to investigate the use of region-specific estimates of (co)variances in genomic predictions using ssSNPBLUP.

Data sets and simulations
The genotype data were simulated for five generations (Gen1−Gen5) based on real haplotypes of 2200 Holsteins (Gen0), as described in Karaman et al. (2018). At each generation, the number of males and females were kept constant at 200 and 2000, respectively, and the mating ratio was 1:10. Mating was completely at random, and selection was not considered. Each sire was mated twice with one of the ten dams to keep the population size at 2200 at each generation. Only the single nucleotide polymorphisms (SNPs) (11,154) located on first five chromosomes were considered.
Phenotypic values of the two traits were simulated to have heritabilities of 0.1 and 0.4, which represents low (L) and high (H) heritability traits, respectively. Total number of quantitative trait loci (QTL) was set at 200, which were randomly selected from the SNP set, ensuring that the average minor allele frequency (MAF) of QTL is 0.15 (Karaman et al. 2018). The criterion for the MAF of the QTL was based on the assumption that they in general have relatively low MAF Kemper and Goddard 2012). The QTL were randomly assigned into three groups according to their causal relationships with the traits. This was done by assuming a percentage of the total QTL (82%) had pleiotropic effects on two traits, while one half of the remaining QTL had effect on one trait, and one half on the other trait.
Two scenarios, G9 and N5, were considered in terms of the distribution of QTL effects and correlations for the effect of pleiotropic QTL. In the scenario G9, the effects of the pleiotropic QTL were achieved by simulating two correlated gamma variables (Dvorkin 2012) with marginal distributions of G(0.4, 1.66), and a correlation of 0.9. The 78% of those QTL were assigned to a correlation between effects on two traits of 0.9, and 22% of −0.9 randomly. The correlation group of −0.9 was achieved by switching the sign of QTL effect for one of the traits at random. The QTL effects, which were assumed to have a correlation of 0.9, were assigned a negative or positive sign at random for both traits. In the second scenario, scenario N5, effects of all pleiotropic QTL were simulated from a bivariate normal distribution with a correlation of 0.5. Although fluctuated across the replicates, all scenarios lead to genetic correlations of about 0.45 at Gen0. The QTL SNPs were excluded from the final data set of SNP for the analysis. Random residual effects were sampled from N 0; Iσ 2 where the sizes of σ 2 e L and σ 2 e H were determined according to heritabilities of 0.1 and 0.4, respectively.
Final data (see Table 1) were created by masking genotypes and/or phenotypes of the animals as follows. For generations 3 and 4, it was assumed that males had no phenotypes, but genotypes, while all females had phenotypes, and some fraction of them had also genotypes. Those genotyped females were selected completely at random. Generation 5 was used as validation population, where 500 randomly selected animals were assumed to be genotyped. Pedigree was traced back to Gen0. Animals had phenotypes on both traits, or none of them. In total, 20 replicates were generated.

Models and methods
A novel multi-trait Bayesian whole genome regression model (BayesN0), single-step SNPBLUP and single-step Bayesian regression models introduced by Fernando et al. (2014) were compared for multi-trait genomic prediction. Single-trait analysis were also performed, but neither the models nor their theory were given in this paper, as the models are special cases of their multi-trait counterparts. In this section, we followed the notation in Fernando et al. (2014) as closely as possible.

Basic multi-trait model
A multi-trait mixed model including only general means as fixed effects and marker effects as random effects can be written as where y L and y H are the vectors of phenotypes, 1 are vectors of ones, μ L and μ H are general means, M L and M H are the matrices of genotypes for k markers, α L and α H are the G genotype, P phenotype, M male, F female vectors of marker effects, and e L and e H are the vectors of random residual effects, for traits "L" and "H", respectively. In our simulations, animals had records for both traits or none of them. Therefore, M L = M H , and these matrices will be denoted as M hereinafter, to simplify the demonstration. Residuals, e′ ¼ e ′ L ; e ′ H Â Ã , are typically assumed to follow a normal distribution, e | R 0~N (0, R 0 ⊗ I), where , and I is an identity matrix.

Multi-trait Bayesian partitioned regression (BayesN0)
The columns of M and vector α given in Eq.
(1) can be divided into S subsets in a conformable manner: where y ¼ y L y H ! involves the phenotypes of genotyped individuals only, M 1 , …, M S are genotype matrices regarding genomic regions, and α t,1 , …, α t,S (t = L, H for low and high heritability traits, respectively) are vectors of SNP effects for corresponding genomic regions. We assume that all SNPs j (j = 1, …, k s ) in the same genomic region s (s = 1, …, S) have the same (co)variance for the two traits: Likelihood of the model is given as: with B i being diagonal matrices consisting of SNP variances (B L and B H ) or covariances (B LH ¼ B HL ), and R ¼ R 0 I. The vector of fixed effects, μ, were assigned a flat prior, and other parameters of the model were assigned a normal or an inverse Wishart (IW) prior for conjugacy: Full conditional distributions of μ, α sj , B s , and R 0 can be obtained after some algebra: where, "." stands for all other parameters and y * , y * is the vector of phenotypes corrected for all other effects, i . This multi-trait whole genome regression model was referred to as multi-trait BayesN0 throughout this paper, as it is an extension of a particular form of partitioned regression model (BayesN) introduced by Zeng et al. (2016), to multi-trait case. Note that when the size of region is fixed at one SNP or whole genome, model becomes equivalent to multi-trait BayesA or GBLUP, respectively.

Multi-trait single-step SNPBLUP
In the following expressions, n stands for the nongenotyped animals, and g stands for the genotyped animals. Note that in our simulations, animals had records for both traits or none of them. In a multi-trait single-step SNPBLUP (ssSNPBLUP) analysis, the phenotypes are modeled as (Fernando et al. 2014): where y ¼ y L y H ! is the vector of phenotypes for genotyped and non-genotyped individuals, μ Ã ¼ 6 4 3 7 7 5 , μ L and μ H are the overall means of the two traits, μ g,L and μ g,H are the differences between breeding values of genotyped and nongenotyped animals for the two traits, . Z n and Z g are incidence matrices relating breeding values of non-genotyped and genotyped animals to their phenotypes,M n and M g are matrices of imputed and observed genotypes for nongenotyped and genotyped animals, respectively, α L and α H are the vectors of allele substitution effects. The U ¼ ϵ L and ϵ H are the vectors of imputation residuals. The e is a vector of random residual effects assumed to follow e | R 0Ñ , and I is an The A ng , A gg and A nn are submatrices of the pedigree-based relationship matrix, A, corresponding to the relationships between non-genotyped and genotyped individuals, among the genotyped individuals, and among the non-genotyped individuals, respectively. The matrix of imputed genotypes,M n , is obtained with A ng A À1 gg M g (Fernando et al. 2014).
The mixed model equations corresponding to the model in Eq. (2) is as follows.
The A −nn is the part of the inverse of pedigree-based relationship matrix, A, corresponding to the non-genotyped individuals, and R = R 0 ⊗ I.

Multi-trait single-step BayesN0 (ssBayesN0)
The single-step SNPBLUP requires the estimation of (co) variance components, and then use of these in mixed model equations to estimate breeding values. In contrast, Bayesian approach can be used to obtain the vector of fixed and random effect estimates, b μ; b α; b ϵ ½ ′ , the genetic and residual variance components, and the SNP (co) variances simultaneously, as in the original paper of Fernando et al. (2014). In principle, any Bayesian whole genome regression model can be incorporated in this single-step model, and BayesN0 was used here (ssBayesN0). Likelihood of the ssBayesN0 model is given as: where matrices and parameters are as specified earlier. A flat prior was assumed for μ*. Priors for α, e, B s and R 0 were the same as in BayesN0. A multivariate normal prior, ϵ j G 0 ; A $ Nð0; G 0 AÞ, was assumed for the vector of ϵ, and G 0 was assigned an inverse Wishart prior, G 0 , R 0 can be obtained after some algebra: where y * , S B S and S R are as defined before,

Statistical analysis
Single-and multi-trait models of BayesN0 and single-step BayesN0 (ssBayesN0) were fitted with varying region sizes (one SNP, 100 SNPs, a whole chromosome and the whole genome). The parameters of the priors for SNP, residual and genetic (co) variance matrices in the multi-trait models were which were derived from the mean of an inverse Wishart distributed random variable, and v B = v R = v G = 5. It is worth noting that inverse Wishart distribution imply a scaled inverse chi-square distribution for each variance with specific parameters (Wang et al. 2018). That is, e.g., Single-trait BayesN0 and ssBayesN0 models were special cases of their multi-trait counterparts, for which the multivariate normal priors for SNP effects, model residuals and imputation residuals were replaced with univariate normal priors e:g:; α L;sj $ N 0; σ 2 α L;s , and inverse Wishart priors for the (co)variance components were replaced with scaled inverted chi-square priors e:g:; σ 2 α L;s $ χ À2 df; S 2 L À Á , for conjugacy. Parameters for these scaled inverted chi-square prior distributions for SNP, residual and genetic variances were df = 4 and a scale parameter, derived from the expected value of a scaled inverse chi-square distributed random variable e:g: Habier et al. 2010a). That is, e.g., σ 2 α L;s $ χ À2 4;σ 2 α L;s 2 . Hence, not only the mean, but also the distribution of priors for the variances were consistent between the single-and multi-trait analysis, with only difference being the value of variance components used. The matrices of e G 0 and e R 0 used in priors for multi-trait analysis, and geneticσ 2 g and residual variancesσ 2 e À Á used in priors for single-trait analysis, were the estimates obtained by fitting single or multi-trait Ridge-Regression models at SNP level, respectively, using the JWAS (Cheng et al. 2018a) package in Julia (Bezanson et al. 2017).
Markov-chain Monte Carlo (MCMC) algorithm with Gibbs sampling method was used to obtain samples of each parameter from its full conditional posterior distribution. Chain length for the analyses using BayesN0 and ssBayesN0 consisted of 50,000 or 70,000 cycles, of which the first 30,000 or 50,000 cycles were discarded as burn-in, respectively. Convergence was tested by comparing results for the two chain lengths (50,000 vs. 70,000) on a random subset of the replicates and region sizes (Zeng et al. 2018). Every tenth sample of the post burn-in cycles were stored for posterior analysis, yielding 2,000 posterior samples. Mean value of the posterior samples was used as the estimate of each parameter. The change in accuracy of prediction was negligible for 70,000 compared to 50,000 cycles of Markov chain, and therefore, the results from the chain length of 50,000 were presented.
For single-and multi-trait ssSNPBLUP models, the genetic and residual (co)variances and SNP (co)variances were obtained as the mean values of the posterior samples from BayesN0 or ssBayesN0. The genetic (co)variances required in mixed model equations forε were computed as the mean of the (co)variances of the breeding values at each MCMC cycle for BayesN0, or directly as the mean of genetic (co)variances for ssBayesN0. Hereafter, analysis using the variance components from BayesN0 and ssBayesN0 will be referred to as ssSNPB1 and ssSNPB2, respectively. The ssSNPB1 and ssSNPB2 models were solved with the Conjugate Gradients method with diagonal preconditioning using the IterativeSolvers package in Julia, and convergence tolerance was chosen to be 10 −12 . All analyses were performed using self-written scripts in Julia.
The predicted breeding values of animals using multitrait BayesN0 were obtained from The predicted breeding values of animals using single-step models, ssBayesN0, ssSNPB1 and ssSNPB2, were obtained from: Prediction accuracy was assessed as the correlation between true and predicted breeding values of validation individuals. The bias of prediction was assessed based on the slope of the regression of true breeding values on the estimated breeding values of validation individuals. Accuracy for single-and multi-trait models with different region sizes were compared for each trait, and each model separately. Prediction accuracy for all methods was compared for each trait and at each scenario of region size. All comparisons were performed separately for genotyped and non-genotyped individuals using a two-sided paired t-tests, for which accuracies were paired across each replicate for the same validation population. A Bonferroni correction was used to control the Type 1 error rate of 0.05, caused by multiple comparisons.

Bayesian whole genome regression (BayesN0)
Prediction accuracies from single-and multi-trait BayesN0 models are given in Tables 2 and 4 for genotyped individuals in validation population, at varying sizes of genome region. Grouping 100 adjacent SNPs generally provided the highest accuracies for both single-and multi-trait models, with some exceptions in scenario N5. Accuracies for different region sizes were generally ranked as 100 SNPs > 1 SNP > 1 Chr > WG in scenario G9. When a multi-trait model was used in scenario G9, prediction accuracy for the region size of 100 SNPs were about 4 and 12 percentage points higher for low heritability trait (L), and about 3 and 8 percentage points higher for high heritability trait (H), compared to those for region sizes of one SNP (BayesA) and whole genome (GBLUP), respectively. Using multitrait BayesN0 with a region size of 100 SNPs resulted in higher accuracies than corresponding single-trait BayesN0 for both traits, though not always significant. Bias for predicting breeding values of genotyped individuals is shown in Supplementary Tables S1 and S3. Regression coefficients were generally closer to 1 for trait H in both scenarios. They were higher than 1 particularly for single-trait analysis of trait L in scenario G9 and single-and multi-trait analysis of trait L in scenario N5.

Single-step genomic prediction
Prediction accuracies from single-and multi-trait analysis are given in Tables 2-5. Similar to BayesN0, accuracies for different region sizes were generally ranked as 100 SNPs > 1 SNP > 1 Chr > WG in scenario G9. For single-trait analysis of trait L, accuracies from the region size of 1 SNP and/or WG were similar to, or even slightly higher than, that of region size of 100 SNPs in scenario N5. Using ssSNPB1 improved accuracies for genotyped individuals compared to using BayesN0, for both single-and multitrait analysis. Accuracies from ssBayesN0 were generally similar to or somewhat higher than those from ssSNPB1, particularly in scenario G9. Using single-step SNPBLUP with (co)variances obtained from ssBayesN0, i.e., ssSNPB2, yielded similar accuracies to the corresponding ssBayesN0 model. Accuracies from ssSNPB2 were similar to, though sometimes slightly higher in scenario G9, those from ssSNPB1 for non-genotyped animals. For nongenotyped animals, taking 100 adjacent SNPs as a genome region provided similar to or slightly higher accuracies   Different alphabets mean significantly different values at a Type 1 error rate of 0.05 with Bonferroni correction. Subscripts and superscripts stand for comparisons within column and row, respectively, for each trait than taking one SNP as a genome region, but higher accuracies than taking whole genome as a genome region, in scenario G9. For scenario N5, on the other hand, all region sizes generally lead to similar accuracies for nongenotyped animals. Regression coefficients were generally closer to 1 for trait H, but higher than 1 for trait L in scenario N5 (Supplementary Tables S1-S4).

Single-vs. multi-trait genomic prediction
Multi-trait analysis generally led to higher accuracies than their single-trait counterparts for trait L (h 2 = 0.1), and similar to or higher accuracies than their single-trait counterparts for trait H (h 2 = 0.4) (Tables 2-5). This was expected because the gain of accuracy from multi-trait over single-trait genomic prediction is more profound for low heritability traits that are genetically correlated with a high heritability trait (Jia and Jannink 2012;Guo et al. 2014). Hayashi and Iwata (2013) compared accuracies from singleand multi-trait analysis for traits with a genetic correlation of 0.7, and reported that accuracy for a low heritability trait (h 2 = 0.1) was improved with multi-trait analysis, while accuracy for a high heritability (h 2 = 0.8) trait remained unchanged. For a low heritability (h 2 = 0.05) trait, which had incomplete data, Guo et al. (2014) showed that accuracy of genomic prediction was improved when a genetically correlated (r g = 0.5) trait with high heritability   Cheng et al. (2018b) reported that the mean of the posterior probability that a marker has a null effect was higher (0.97 vs. 0.74) in multi-trait analysis (BayesCΠ) compared to single-trait analysis (BayesCπ) for gall volume (h 2 = 0.12), when the correlated trait was presence (or absence) of rust (h 2 = 0.21), in Loblolly Pine (Pinus taeda L.) (Resende et al. 2012). Beside heritability, another factor influencing accuracy is the absolute difference between genetic and residual correlations (Schaeffer 1984;Thompson and Meyer 1986). In this study, the simulated residual correlation was null and the genetic correlation was moderate (0.45), though the estimates of those correlations varied around the simulated true values. Averaged over the replicates, genetic correlations were generally overestimated, whereas the residual correlations were nearly zero and varied only after second decimal, in both scenarios and for all region sizes. Genetic correlations were 0.47 and 0.45 from BayesN0 with the region size of 100 SNPs, and 0.56 and 0.51 from GBLUP (BayesN0 with whole genome as one region), for scenarios G9 and N5, respectively (results not given elsewhere). Those were 0.49 and 0.47 for ssBayesN0 with the region size of 100 SNPs, and 0.54 and 0.48 for ssGBLUP (ssBayesN0 with whole genome as one region), for scenarios G9 and N5, respectively (results not given elsewhere). These small deviations of genetic correlations from their true values are expected to have little influence in variance of prediction error (PEV), and multi-trait models can increase the precision of breeding value estimates by reducing PEV compared to single-trait models (Schaeffer 1984). The PEV was additionally computed for BayesN0 and ssBayesN0, from the variance of posterior samples for breeding values of genotyped individuals in validation population. Averaged over region sizes, the mean reduction in PEV from multi-trait BayesN0 were about 2.5% for trait L and 0.5% for trait H, and 5% for trait L and 0.5% for trait H, in scenarios G9 and N5, respectively (results not given elsewhere). The mean reduction in PEV from multi-trait ssBayesN0 were about 9% for trait L and 0.9% for trait H, and 6% for trait L and 0.8% for trait H, in scenarios G9 and N5, respectively (results not given elsewhere). Bias for single-trait analysis was relatively high for trait L particularly in scenario G9 (Supplementary Tables S1-S4), however, it was generally reduced by using multi-trait models.
In multi-trait genomic prediction, correlation structures between the traits is central to gaining advantage in prediction accuracy over single-trait predictions (Gebreyesus et al. 2017). Our results showed that the improvement from multi-trait analysis over single-trait analysis were dependent on whether the genetic makeup of the (co)variance structure of the studied traits (Tables 2-5) were accounted for, and this will be discussed in detail in the later sections.
Accounting for heterogeneous (co)variances across the genome using BayesN0 Multi-trait genomic prediction rely on the genetic association between the traits through the genetic variances and covariances, which may vary across the genome. A few genome regions may explain a substantial proportion of the covariance, whereas others account for nearly no covariance between the traits (Sørensen et al. 2012). Moreover, covariances between particular traits may be positive for some regions and negative for others, while the overall genetic correlations are low/high (Li et al. 2017;Gebreyesus et al. 2017). This study investigated the affect of assigning priors to genome regions, which were defined as fixed number of SNPs (one SNP, 100 SNPs, one chromosome or whole genome), on accuracy in multi-trait genomic prediction.
Genomic prediction rests on the LD between QTL and SNPs (Meuwissen et al. 2001). Although the simulation settings in this study resulted in correlations of QTL effects that fall into different categories, it may be of a general question where does the heterogenity of (co)variances over the genome come from, or what does it refer to. It can be shown that the best linear predictor of SNP effects is α where γ t is the vector of QTL effects, V M is the (co)variance matrix of SNP genotypes, and V MQ is the covariance matrix of SNP and QTL genotypes (de los Campos et al. 2015). Note that for a QTL that affect only L (or H), corresponding row of γ H (or γ L ) is zero. Under some assumptions, (co)variance of the SNP effects are proportional to V M s Q s V ′ M s Q s , for genome region s (s = 1, …, S). Because recombination rates vary over the genome, and SNPs are typically in imperfect LD with QTL, each V M s Q s may be different (Wang et al. 2013), resulting in genome having a different (co)variance pattern at the SNP level than that of at the QTL level (de los Campos et al. 2015).
Multi-trait BayesA (BayesN0 with region size of one SNP) was able to account for the heterogeneous correlation structure across the genome to some extent, compared to multi-trait GBLUP (BayesN0 with whole genome as one region), which assumes a constant correlation across the genome (Tables 2 and 4). Accuracies were further improved when a group of 100 SNPs were allowed to have a common (co)variance. It should be noted that the choice of region sizes was arbitrary, and therefore, the region size of 100 SNPs may not be optimal. Alternatively, regions can be achieved by grouping SNPs based on fixed length of genomic region or LD information. Because the extent of LD is highly variable in different populations (Wang et al. 2013), and varies with respect to SNP density , the decision of optimal region size is crucial to obtain highest accuracy of genomic prediction (Gebreyesus et al. 2017).
Simulation studies have shown that Bayesian whole genome regression models, which allow variances of SNP effects differing among loci or genome regions, perform better than GBLUP model (Meuwissen et al. 2001;Lund et al. 2009;Karaman et al. 2018). In real-data applications, the accuracy of genomic prediction using the Bayesian whole genome regression models led to similar to or higher accuracies than methods assuming a constant variance structure (e.g., GBLUP) across the genome Habier et al. 2010b;Su et al. 2012a). The benefit from Bayesian whole genome regression models was larger for traits with simple genetic architectures (Coster et al. 2010;Daetwyler et al. 2010b;Clark et al. 2011;Karaman et al. 2018). Examples for such traits can be milk protein composition traits, in which a substantial proportion of the variance is explained by a few QTL (Heck et al. 2009;Schopen et al. 2011). Gebreyesus et al. (2017 reported that BaysesAS model resulted in higher prediction reliabilities than GBLUP for milk protein composition traits, when 100 SNPs were assumed to have a common (co)variance, based on a data set from 50 K SNP panel in Danish Holstein cattle.
For prediction of traits with large effect QTL, the GBLUP model, in which a selection of SNPs, i.e., SNPs identified in earlier genome-wide association studies (GWAS) or identified via GWAS using the current data (de novo GWAS, Spindel et al. (2016)), are considered as fixed effects, can provide accuracies as high as those from Bayesian whole genome prediction methods (Spindel et al. 2016;Lopes et al. 2017). Although the approach is relatively straightforward, it either requires a priori information about the SNPs for the traits of interest, or running a GWAS prior to genomic prediction. Depending on the choice of statistical method, the definition of the QTL region and the significance threshold, different sets of SNPs can be achieved even with the same data, and QTL regions that explain a substantial proportion of the variance may also not always be identified for all traits (Goddard et al. 2016;Lopes et al. 2017).
By applying BayesN0, one has the possibility of putting emphasis on genome regions with large effect, without requiring any prior knowledge on the QTL region affecting the trait(s), or without running a de novo GWAS (Lopes et al. 2017). For practical application in breeding programs, we think this is an advantage over GBLUP, in which "some" SNPs are considered as fixed effects. Our results for scenario N5 imply that the advantage of grouping SNPs in BayesN0 over GBLUP is not limited only to traits with a few QTL with large effect and many with small effects (scenario G9). In scenario N5, BayesN0 with 100 SNPs region size led to a similar accuracy to that from GBLUP in single-trait analysis of trait L, but to a higher accuracy than GBLUP in multi-trait analysis of trait L. Since using a multi-trait model may be beneficial for traits by increasing the amount of information, it can be argued that the accuracies from single-trait analysis of trait L would also differ among the region sizes, for the intermediate sizes of data (Karaman et al. 2016). For asymptotically large sizes of data, on the other hand, there might be little or no benefit of using more sophisticated methods compared to GBLUP (Karaman et al. 2016;Cheng et al. 2018b).
In a simulation study for single-trait genomic prediction, Zeng et al. (2018) showed that BayesN was superior to BayesB when the QTL had relatively low MAF, for a panel consisting of 50 K SNPs. It is, however, unclear if this was due to selection of regions at each cycle of MCMC, or due to reliable estimation of SNP variances by assuming common variance to SNPs in each region, rather than assuming a variance specific to each SNP. In that study, fitting ten SNPs per region also provided higher accuracies of prediction than fitting two SNPs per region. Hess et al. (2017) further allowed SNPs within a region to have different variances, in a study using 50 K SNP panel of an admixed cattle population in New Zealand. There was no advantage of BayesN over BayesB, for milk fat yield, live-weight and somatic cell score. They also showed that fitting all SNPs in a region resulted in slightly higher accuracies than fitting only two SNPs per region.
Accounting for heterogeneous (co)variances across the genome using single-step Bayesian regression  (Tables 2 and 4). In a single-step analysis using Bayesian regression (Fernando et al. 2014), taking one SNP as a genome region is equivalent to single-step BayesA (ssBayesA) and taking whole genome as one region is equivalent to single-step GBLUP (ssGBLUP). Our results indicate that ssBayesA can lead to higher accuracies than ssGBLUP in a multi-trait analysis, by exploiting the heterogeneous (co)variance structure across the genome. However, similar to the regular BayesA (Meuwissen et al. 2001), the information in the data that is utilized by ssBayesA is limited, due to its strong dependency on the prior for (co)variance of SNP effects (Gianola et al. 2009). This dependency on the prior was overcome to some extent by assuming a common (co)variance for 100 adjacent SNPs using ssBayesN0, which generally led to higher accuracies than ssBayesA and ssGBLUP for both trait L and H (Tables  2 and 4). Similarly, prediction accuracy for non-genotyped individuals were increased about 6 and 0.5 percentage points for trait L, and about 2.5 and 0.6 percentage points for trait H, for scenarios G9 and N5, respectively, when region size was changed from whole genome to 100 SNPs in multi-trait analyses (Tables 3 and 5).
For genotyped individuals, using multi-trait ssBayesN0 led to higher accuracies than using multi-trait BayesN0 (Tables 2 and 4). As other Bayesian whole genome regression models, BayesN0 can only use the phenotypes of genotyped animals. The ssBayesN0, on the other hand, simultaneously uses the phenotypes of genotyped animals (1000) and non-genotyped animals (3000, the genotypes were imputed) ( Table 1) to estimate the SNP effects in the MCMC procedure, while accounting for the error in imputation for non-genotyped individuals with phenotypes. This enhances the data size used in estimation of SNP effects, which has a key role to obtain reliable prediction of breeding values (Daetwyler et al. 2008;Goddard 2009;Karaman et al. 2016;Cheng et al. 2018b).
Practical implementation of single-step models using previously estimated (co)variance components In this study, the estimates of the (co)variances were obtained from BayesN0 or ssBayesN0. The former led to ssSNPB1 model for which the (co)variance components were obtained using only the information of genotyped individuals, while the latter led to ssSNPB2 model for which the (co)variance components were obtained using the information of genotyped and non-genotyped individuals. In practice, the (co)variance components can be estimated less frequently compared to routine genomic evaluations without harming the prediction accuracies .
For genotyped individuals, ssSNPB1 model yielded higher accuracies than BayesN0, at all region sizes in multitrait analysis (Tables 2 and 4). This was due to more accurate estimation of SNP effects by the use of phenotypes of non-genotyped individuals. The ssBayesN0 and ssSNPB2, where the (co)variance components from ssBayesN0 were used, generally yielded similar accuracies. This was not surprising, because similar to BayesC0 and SNPBLUP being equivalent models, ssBayesN0 and ssSNPB2, are also equivalent. The ssSNPB2 generally led to higher accuracies than ssSNPB1 in scenario G9, and similar to or slightly higher accuracies than ssSNPB1 in scenario N5, in multi-trait analyses.
Accuracies for non-genotyped animals were generally similar among the models, i.e., ssSNPB1, ssBayesN0 and ssSNPB2, in multi-trait analysis. Analysis using the model of Fernando et al. (2014) starts with an explicit imputation of markers for non-genotyped individuals, using pedigree information and genotypes of genotyped relatives. Then, marker effects and imputation residuals (ϵ) accounting for the part of breeding values, which cannot be modeled by imputed markers, are estimated (Gao et al. 2018). Imputation residual is added to marker-based breeding value (sum of individual SNP effects) of non-genotyped individuals to obtain their total breeding values (Fernando et al. 2014). The breeding value of genotyped individuals, on the other hand, is composed only of sum of individual SNP effects. Hence, a change in the accuracy of SNP effect estimates has less impact on the accuracy of breeding value estimates for non-genotyped individuals than for genotyped individuals (Zhou et al. 2018).
One way to account for heterogeneous (co)variance structure in single-step genomic prediction could be to construct weighted G matrices (Zhang et al. 2010), and in turn their weighted H matrix counterparts (Fragomeni et al. 2017) for ssGBLUP. In an earlier study (Karaman et al. 2018), we have shown that weighted multi-trait GBLUP can reach accuracies similar to that of the Bayesian whole genome regression model which was used to derive weights. This was expected, because those "weighted" relationship matrices are indeed implicit to Bayesian whole genome regression methods (Fernando and Gianola 2018;Karaman et al. 2018). A drawback of the approach using weighted relationship matrices is that it requires the computation of a number of relationship matrices which increase with the number of traits in a multi-trait model, and not only the computing time but also the storage of such H matrices might be impractical for genomic prediction using weighted ssGBLUP in routine evaluations. Moreover, compared to ssGBLUP, the equations needed to be solved for ssSNPBLUP does not grow with the number of genotyped individuals, and the inverse of the combined relationship matrix, H, is not needed (Fernando et al. 2014).
We did not focus on the computational (dis)advantages of ssSNPBLUP, nor its convergency properties. Both ssGBLUP and ssSNPBLUP, though, are known to have some computational challenges (Taskinen et al. 2017). Averaged over the scenarios and region sizes, ssSNPB1 and ssSNPB2 models achieved relative convergence of 10 −12 in 200 and 203 iterations (average of two traits) in single-trait analysis, and in 435 and 478 iterations in multi-trait analysis, respectively. It should be noted that these numbers apply only to the current data, and could vary with equivalent formulations of the models (Taskinen et al. 2017) or with a preconditioner other than diagonal used in this study.

Conclusions
In this study, a multi-trait whole genome regression model, BayesN0, was proposed. The model has its equivalent counterparts when the region size is set at one SNP (BayesA) or the whole genome (GBLUP). Our results showed that assigning priors to genome regions defined as fixed number of SNPs, e.g., 100 SNPs, may improve accuracies over BayesA and GBLUP by accounting for heterogeneous (co)variance structure across the genome efficiently. The model was also implemented in single-step (ssBayesN0) Bayesian regression approach, which unifies pedigree, phenotypes and genotypes in a single analysis. Highest prediction accuracies were obtained when 100 adjacent SNPs were assumed to have a common (co)variance in ssBayesN0. For routine genomic evaluations, it could be a good strategy to estimate (co)variance components from ssBayesN0, and then to use those estimates in genomic prediction using multi-trait single-step SNPBLUP. Such a strategy has the potential to provide reliable estimates of breeding values for both genotyped and nongenotyped individuals.

Data availability
Genotype and pedigree data can be found at https://doi.org/ 10.5061/dryad.v4126t4, along with a file including necessary SNP information (chromosome ID and base-pair position). The data and the methodology described previously are sufficient to reproduce the results of this study.