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-of-the-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 pedigree-based 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 region-specific variances, and 1 − π 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 Π = π = 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 multi-trait 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). Multi-trait 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 (Misztal and Legarra 2017). 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 (Su et al. 2014). 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.

Material and methods

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 (Goddard and Hayes 2009; 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 , where the sizes of and 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.

Table 1

Number of animals with genotype and phenotype in each generation

Generation	G	G&P	P	Total
Gen0	–	–	–	2200
Gen1	–	–	–	2200
Gen2	–	–	–	2200
Gen3	200 (M)	500 (F)	1500 (F)	2200
Gen4	200 (M)	500 (F)	1500 (F)	2200
Gen5	500	–	–	2200

G genotype, P phenotype, M male, F female

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 aswhere yL and yH are the vectors of phenotypes, 1 are vectors of ones, μL and μH are general means, ML and MH are the matrices of genotypes for k markers, L and H are the vectors of marker effects, and eL and eH 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, ML = MH, and these matrices will be denoted as M hereinafter, to simplify the demonstration. Residuals, , are typically assumed to follow a normal distribution, e | R0 ~ N(0, R0 ⊗ 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 involves the phenotypes of genotyped individuals only, M1, …, M are genotype matrices regarding genomic regions, and , …, (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) 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:where , , with being diagonal matrices consisting of SNP variances ( and ) or covariances (), and . 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 , , B, and R0 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, , and . 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 non-genotyped 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 is the vector of phenotypes for genotyped and non-genotyped individuals, , μ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 non-genotyped animals for the two traits, with , with and , . Zn and Zg are incidence matrices relating breeding values of non-genotyped and genotyped animals to their phenotypes, and Mg are matrices of imputed and observed genotypes for non-genotyped and genotyped animals, respectively, L and H are the vectors of allele substitution effects. The with and , where and are the vectors of imputation residuals. The e is a vector of random residual effects assumed to follow e | R0 ~ N(0, R0 ⊗ I), where , and I is an identity matrix. Vector of is assumed to follow | B ~ N(0, B) with , where B are diagonal matrices consisting of SNP variances (BL and BH) or covariances (BLH = BHL). Vector of is assumed to follow , where G0 is the additive genetic (co)variance matrix. The Ang, Agg and Ann 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, , is obtained with (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 = R0 ⊗ 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, , 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 and R0 were the same as in BayesN0. A multivariate normal prior, , was assumed for the vector of , and G0 was assigned an inverse Wishart prior, G0 | vG, VG ~ IW(vG, VG). Full conditional distributions of *, , B, , G0, R0 can be obtained after some algebra:where y*, and S are as defined before, , Nn is the number of non-genotyped individuals, and .

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 where , , and , which were derived from the mean of an inverse Wishart distributed random variable, and v = v = v = 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., , where is the first diagonal element in .

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 , and inverse Wishart priors for the (co)variance components were replaced with scaled inverted chi-square priors , 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 (Habier et al. 2010a). That is, e.g., . 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 and used in priors for multi-trait analysis, and genetic and residual variances 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 multi-trait BayesN0 were obtained fromThe 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.

Results

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 multi-trait 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.

Table 2

Accuracies for genotyped individuals using single- and multi-trait models in scenario G9

Trait1	Region size2	Single-trait3				Multi-trait
		BayesN0	ssSNPB1	ssBayesN0	ssSNPB2	BayesN0	ssSNPB1	ssBayesN0	ssSNPB2
L4	1 SNP	ab0.349e	ab0.452cd	ab0.470d	ab0.478c	b0.437d	b0.536b	b0.554b	a0.574a
	100 SNPs	a0.365d	a0.460c	a0.478c	a0.479c	a0.481c	a0.559b	a0.590a	a0.590a
	1 Chr	b0.335c	b0.434b	bc0.444b	bc0.445b	b0.402b	c0.493a	c0.497a	b0.499a
	WG	b0.335d	b0.433b	c0.433bc	c0.433b	c0.362cd	d0.461ab	d0.472a	c0.473a
H	1 SNP	b0.587g	b0.683e	b0.689de	b0.700bc	b0.593f	b0.688cd	b0.699b	a0.712a
	100 SNPs	a0.611f	a0.698d	a0.716b	a0.716b	a0.622e	a0.707c	a0.725a	a0.725a
	1 Chr	c0.552e	c0.650bd	c0.651cd	c0.651abcd	c0.558e	c0.655ac	c0.657ab	b0.657ab
	WG	d0.538d	c0.642c	c0.642bc	c0.643bc	d0.543d	d0.644abc	d0.646ab	c0.646a

1L and H: low (0.1) and high (0.4) heritability traits, respectively

2Chr chromosome, WG whole genome

3ssSNPB1 and ssSNPB2: Single-step SNPBLUP, for which the variance components were obtained from BayesN0 and ssBayesN0, respectively

4Different 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

Table 4

Accuracies for genotyped individuals using single- and multi-trait models in scenario N5

Trait1	Region size2	Single-trait3				Multi-trait
		BayesN0	ssSNPB1	ssBayesN0	ssSNPB2	BayesN0	ssSNPB1	ssBayesN0	ssSNPB2
L4	1 SNP	a0.314d	a0.432b	a0.432b	a0.433b	a0.362c	a0.470a	a0.469a	a0.470a
	100 SNPs	a0.313e	a0.428bc	a0.434b	ab0.429c	a0.367d	ab0.468a	a0.475a	a0.474a
	1 Chr	a0.309d	b0.419c	b0.420c	b0.419c	b0.341d	c0.447abc	b0.447b	b0.450a
	WG	a0.314c	a0.431ab	ab0.430b	a0.432a	b0.342c	bc0.447ab	ab0.459ab	ab0.460ab
H	1 SNP	b0.545d	a0.651c	b0.654bc	a0.658ab	b0.548d	a0.654bc	b0.657b	a0.662a
	100 SNPs	a0.554d	a0.654c	a0.663a	a0.662ab	a0.559d	a0.657bc	a0.666a	a0.665a
	1 Chr	c0.537b	b0.642a	c0.645a	b0.645a	c0.540b	b0.644a	c0.646a	b0.647a
	WG	c0.537b	b0.644a	c0.645a	b0.646a	c0.539b	b0.645a	c0.648a	b0.648a

1L and H: low (0.1) and high (0.4) heritability traits, respectively

2Chr chromosome, WG whole genome

3ssSNPB1 and ssSNPB2: Single-step SNPBLUP, for which the variance components were obtained from BayesN0 and ssBayesN0, respectively

4Different 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

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 multi-trait 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 non-genotyped animals, taking 100 adjacent SNPs as a genome region provided similar to or slightly higher accuracies 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 non-genotyped 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).

Table 5

Accuracies for non-genotyped individuals using single- and multi-trait models in scenario N5

Trait1	Region size2	Single-trait3			Multi-trait
		ssSNPB1	ssBayesN0	ssSNPB2	ssSNPB1	ssBayesN0	ssSNPB2
L4	1 SNP	a0.328c	a0.325c	a0.326c	a0.357a	a0.353b	a0.355ab
	100 SNPs	a0.327b	a0.327b	a0.326b	a0.357a	a0.355a	a0.355a
	1 Chr	a0.324c	a0.325c	a0.324c	ab0.349ab	a0.347b	a0.350a
	WG	a0.327b	a0.325b	a0.327ab	b0.343ab	a0.350b	a0.352a
H	1 SNP	a0.506d	b0.507cd	a0.510b	a0.508bc	ab0.509b	a0.512a
	100 SNPs	a0.507c	a0.511ab	a0.511ab	a0.509bc	a0.512ab	a0.512a
	1 Chr	ab0.503b	bc0.504ab	b0.504ab	ab0.505a	bc0.505ab	b0.506ab
	WG	b0.503b	c0.503b	b0.503b	b0.504ab	c0.506a	b0.506a

1L and H: low (0.1) and high (0.4) heritability traits, respectively

2Chr chromosome, WG whole genome

3ssSNPB1 and ssSNPB2: Single-step SNPBLUP, for which the variance components were obtained from BayesN0 and ssBayesN0, respectively

4Different 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

Discussion

Single- vs. multi-trait genomic prediction

Multi-trait analysis generally led to higher accuracies than their single-trait counterparts for trait L (h2 = 0.1), and similar to or higher accuracies than their single-trait counterparts for trait H (h2 = 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 single- and multi-trait analysis for traits with a genetic correlation of 0.7, and reported that accuracy for a low heritability trait (h2 = 0.1) was improved with multi-trait analysis, while accuracy for a high heritability (h2 = 0.8) trait remained unchanged. For a low heritability (h2 = 0.05) trait, which had incomplete data, Guo et al. (2014) showed that accuracy of genomic prediction was improved when a genetically correlated (rg = 0.5) trait with high heritability (h2 = 0.3) was available. 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 (h2 = 0.12), when the correlated trait was presence (or absence) of rust (h2 = 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 (t = L, H), where is the vector of QTL effects, VM is the (co)variance matrix of SNP genotypes, and VMQ 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 , for genome region s (s = 1, …, S). Because recombination rates vary over the genome, and SNPs are typically in imperfect LD with QTL, each 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 (Goddard and Hayes 2009), 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 (Hayes et al. 2009; 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

Implementation of our novel Bayesian multi-trait model (ssBayesN0) using the methodology of Fernando et al. (2014) yielded accuracies for genotyped individuals in the range of 0.47–0.59 and 0.45–0.48 for trait L, and 0.65–0.73 and 0.65–0.67 for trait H, in scenarios G9 and N5, respectively (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).

Table 3

Accuracies for non-genotyped individuals using single- and multi-trait models in scenario G9

Trait1	Region size2	Single-trait3			Multi-trait
		ssSNPB1	ssBayesN0	ssSNPB2	ssSNPB1	ssBayesN0	ssSNPB2
L4	1 SNP	ab0.351c	ab0.357c	ab0.361c	b0.402b	b0.406b	a0.418a
	100 SNPs	a0.355c	a0.363c	a0.363c	a0.412b	a0.426a	a0.427a
	1 Chr	b0.340b	b0.342b	bc0.342b	c0.378a	c0.377a	b0.378a
	WG	b0.337c	b0.337c	c0.336c	d0.361abc	d0.365b	c0.366a
H	1 SNP	a0.526cd	a0.528d	a0.531bc	a0.529bcd	a0.533b	a0.537a
	100 SNPs	a0.530c	a0.535b	a0.535b	a0.534b	a0.539a	a0.539a
	1 Chr	b0.513abc	b0.513c	b0.513bc	b0.516abc	b0.517ab	b0.517a
	WG	b0.511bc	b0.511c	b0.511bc	b0.513abc	b0.514ab	c0.514a

1L and H: low (0.1) and high (0.4) heritability traits, respectively

2Chr chromosome, WG whole genome

3ssSNPB1 and ssSNPB2: Single-step SNPBLUP, for which the variance components were obtained from BayesN0 and ssBayesN0, respectively

4Different 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

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 (Su et al. 2014).

For genotyped individuals, ssSNPB1 model yielded higher accuracies than BayesN0, at all region sizes in multi-trait 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.

Supplementary tables