To a large extent we follow the notation of Runcie and Cheng (2019), assuming observations on traits Y₁, …, Y_p+1, where each Y_j is a column vector. The first one (Y₁ = Y_f) is the focal or target trait, for which genomic predictions are required; Y₂, …, Y_p+1 are the secondary traits. Ys=(Y2t,…,Yp+1t)t is the column vector containing all secondary traits; similarly, Y=(Y1t,…,Yp+1t)t is the column vector containing all traits. We have in total n = n_t + n_o genotypes, including n_o genotypes for which the target trait is observed (the training set), and n_t for which it is to be predicted (the t referring to test set). We will use subscripts t and o to indicate that we take the subset of values on the test, respectively training set, for example Y_o and Y_{f, o}.

The secondary phenotypes are either observed only on the training set (the CV1-scenario, using the terminology of Runcie and Cheng, 2019), or also for the test genotypes (CV2). Since our focus here is on variable selection and dimension reduction (rather than different cross-validation schemes), we will refer to these simply with scenarios 1 and 2, respectively. The n × n genetic relatedness matrix K is partitioned as:

where the n_t × n_o matrix K_to defines the relatedness between new (test) and observed (training) genotypes. We will also write K_t· = [K_tt K_to] and K_o· = [K_ot K_oo]. Similarly, we can decompose the genetic and residual covariance matrices Σ^u and Σ^e as

where the scalars Σffu and Σffe are respectively the genetic and residual variance of the focal trait, and the matrices Σssu and Σsse contain the genetic and residual (co)variances of the secondary traits. The row-vectors Σfsu and Σfse contain the genetic and residual covariance between the focal and the secondary traits.

The joint distribution of Y = (Y₁, …, Y_p+1) is assumed to be

where

The genetic covariances (Σfsu) quantify the degree of overlap among genetic signals, based on which multivariate methods can potentially improve genomic prediction. The residual covariances (Σfse) are important when traits are measured on the same individuals; if measured on different individuals (typically, in a different experiment), Σ^e can assumed to be diagonal. Σ^u and Σ^e are usually unknown, and need to be estimated from the data. For p larger than 5 − 10, this usually requires approximations. Below we describe several dimension reduction approaches, which reduce the dimensionality of the secondary phenotypes to 1, and exact REML-estimates of Σ^u and Σ^e can be obtained with standard software.

The main objective is the prediction of the genetic effect U₁ = U_f, i.e., the breeding values for the focal trait, in particular for the test set (U_{f, t}). In our simulations we assess prediction accuracy in terms of the Pearson correlation (r) between the simulated and predicted genetic effects, on the test set. For real data, we consider the correlation between the predicted genetic effects and the trait values observed on the test sets. Although it is well-known that this is a biased estimator of the true accuracy (i.e., the correlation with the unknown genetic effect), the bias is likely to be constant among methods, as long as the target and secondary traits are observed on different plants (Runcie and Cheng, 2019).

The univariate GBLUP for U_{f, t} is defined by

where U^f,o(uni) is the GBLUP for the training set, and REML-estimates of β_f and the variance components Σffu and Σffe are obtained from a univariate mixed model for Y_f. This is the best (univariate) linear unbiased predictor, at least given the true values of the variance components.

The multivariate GBLUP in scenario 1 is

where U^f,o(m1) is the GBLUP for the training set, and REML-estimates of β and the variance components (matrices) Σ^u and Σ^e are obtained from the multivariate mixed model for Y_f and Y_s. As pointed out by Runcie and Cheng (2019), U^f,t(m1) and U^f,t(uni) have the same form, but the “input” Û_{f, o} differs.

The multivariate GBLUP in scenario 2 is

where 0 denotes a n_t × n_o matrix of zeros. This differs from the CV2 prediction in Runcie and Cheng (2019), who described a two-step approach.

Expressions (4) and (5) are valid regardless whether there is just a single secondary phenotype, or multiple ones. However, when the dimension of the secondary phenotype (p) is larger than 5 − 10, estimation of the required genetic covariances quickly becomes challenging and often infeasible (Zhou and Stephens, 2014; Zwiernik et al., 2017). Moreover, even if estimates of genetic covariance are available, the resulting predictions may be prone to overfitting. Reducing the dimension of the secondary phenotype appears to be a relevant strategy to deal with these issues.

Here we propose the dimension reduction S = ĥ(Y_s), where ĥ(Y_s) is a prediction of Y_f based on Y_s, obtained either with LASSO or random forests. Genomic prediction in scenarios 1 and 2 is then performed using (4) and (5), with S = ĥ(Y_s) as secondary trait. We will refer to the resulting genomic predictions using LS-BLUP and RF-BLUP, depending on whether the dimension reduction was achieved by respectively LASSO or random forests. In a GWAS context, such dimension reductions have been used by van Heerwaarden et al. (2015) and Melandri (2019). The intuition behind this dimension reduction is that some of the secondary traits may have a causal effect on Y_f (Figure 1, left). Genomic prediction with LS-BLUP and RF-BLUP may then work well if Ŷ_f captures most of the relevant genetic correlations. In our simulations described below, we also consider the situation where genetic correlations are not the result of a causal effect of Y_s on Y_f (for example, as in Figure 1, right panel). Because of the relatively small size of the populations considered here, the dimension reduction is computed on the same training set that is used for genomic prediction. This is of course not essential for this approach, and various sample splitting techniques may be of interest for larger populations; see the discussion section below.

When using RF-BLUP in the simulations described below, we used the R-package randomForest, with the default settings. Often however, a more accurate dimension reduction can be achieved by tuning various hyperparameters (like the number of trees), which we explore for the real data.

In addition to the notation Y_s for the column vector containing all secondary traits, we will now also use Y_s(j) for the column-vector containing the jth secondary trait, the dimension being either n_o × 1 (scenario 1) or n × 1 (scenario 2). We will use Ys(i) for the row-vector containing all secondary traits for genotype i. Recall that the individual secondary traits are still labeled Y₂, …, Y_p+1, Y₁ being the target trait.

A well-known alternative dimension reduction approach is to use a selection index S=∑j=1pγjYs(j), which is a linear combination of secondary traits, with coefficients such that the resulting index best predicts the genetic effect of the target trait (Falconer and Mackay, 1996). Assuming independent genetic effects (i.e., ignoring population structure), the p × 1 vector γ of coefficients is obtained by minimizing, for each individual i, the expectation of (Uf[i]-Ys(i)γ)2. The minimizing γ then equals the inverse variance-covariance of Y_s times the vector of genetic covariances between Y_s and Y_f, i.e., γSI=Σs-1Σsfu. To estimate γ^SI one could plug in estimates Σ^s and Σ^sfu, where Σ^s=Σ^ssu⊗Koo+Σ^sse⊗Ino is the estimated variance-covariance matrix of the secondary traits on the training population, and Σ^sfu contains estimates of genetic covariances with the target trait. However, when the dimension (p) is large, Σssu and Σsse are difficult to estimate, and the selection index is likely to overfit, as some elements in Σsfu may be large by chance, and receive too much weight.

To address these issues, Lopez-Cruz et al. (2020) proposed penalized selection indices, minimizing instead E(Uf[i]-Ys(i)γ)2+λJ(γ), where λ > 0 is the penalty and J(γ) is either ∑j=1pγj2 (ridge penalty) or ∑j=1p|γj| (LASSO penalty). λ = 0 gives the classical (unpenalized) SI. In case of a ridge penalty, the penalized SI is given by

We will follow the implementation by Lopez-Cruz et al. (2020) in their R-package SFSI, where Σsfu is estimated with MANOVA on the individual plant or plot-level data, and Σssu is estimated using the sample covariance matrix of the secondary traits. We emphasize that no multi-trait mixed-model of the form (1)–(2) is fitted. Moreover, the regularization only controls how Σ^s affects Σ^sfu; the estimates Σ^s and Σ^sfu themselves are not regularized.

Following again (Lopez-Cruz et al., 2020), we use internal cross-validation within the training set to choose an appropriate value of λ, maximizing h(S)ρ_G(S, Y_f). After selecting a value for λ, genomic prediction in scenarios 1 and 2 is performed using (4) and (5), with a single secondary trait, i.e., the selection index ∑j=1pγj(λ)Ys[j]. We will use SI-BLUP to refer to the genomic prediction obtained this way.

Another alternative to selection indices is to model the secondary traits using random effects (see e.g., Riedelsheimer et al., 2012; Van De Wiel et al., 2016; Xu et al., 2016; Schrag et al., 2018; Xiang et al., 2019; Azodi et al., 2020). In addition to the genetic relatedness matrix K, these models use an additional relatedness matrix M derived from the secondary phenotypes, and assume that

where Uf(gen)~N(0,σK2K) and Vf(sec)~N(0,σM2M). We will call this the Multi-BLUP model (not to be confused with Speed and Balding, 2014, where the same type of model is used, but where genomic regions are represented by different relatedness matrices). The variance components σK2, σM2, and σE2 can be estimated with REML or with cross-validation. For simplicity we consider only one type of secondary phenotypes. Similar to the equivalence between GBLUP and SNP-BLUP, the effects Vf(sec) can be written as Y_sb_s, for a vector b_s of independent random effects with N(0,p-1σM2) distribution. Hence, similar to the LS-BLUP and RF-BLUP, the Multi-BLUP approach implicitly assumes a causal effect of Y_s on Y_f (Figure 1, left), which is assumed to be linear, with random coefficients. The usual “genomic” prediction based on model (7) is

i.e., the sum of the BLUPs for the genetic and secondary trait effects. We put genomic between quotes because (8) is partly a phenotypic prediction: instead of the genetic component of the secondary traits, it directly relies on these traits themselves, which are assumed to be available on the test set. As a consequence, the use of (8) is limited to scenario 2.

To overcome these limitations we propose the GM-BLUP:

where b^s is the vector of predicted random coefficients obtained from the Multi-BLUP model, and U^s(gen) is the matrix of GBLUPs for the secondary traits (either univariate or multivariate). These GBLUPs can of course also be computed in scenario 1. Apart from being the “genomic analogue” of (8), (9) can also be motivated by a causal model of the form

as considered by Töpner et al. (2017) and Grotzinger et al. (2019). In contrast to the Multi-BLUP, GM-BLUP only depends on the genetic components of the secondary traits.

Finally, following many other authors (e.g., Riedelsheimer et al., 2012; Xu et al., 2016) we will also compute a prediction based on the secondary traits alone, using the model

and define the MBLUP

Again, this is to some degree a phenotypic prediction, and since the direct effects of the SNPs are ignored, the estimated effects b^s will differ from those obtained from model (7).

We first compare the different methods on simulated data, with p = 300 secondary traits. We used existing genotypic data, from the Arabidopsis RegMap, containing 1, 307 accessions genotyped with 214, 051 SNPs (Horton et al., 2012). For each data-set we randomly selected 500 accessions, from which we randomly sampled a test set of 100 accessions. We randomly selected 1, 500 SNPs with a minor allele frequency of at least 0.3. For each data-set we first simulated direct genetic effects (g_i) and residuals (r_i) for each accession i, and the final trait values were obtained using a structural equation model, describing functional relations between traits. More specifically, for each individual i, the (p + 1) × 1 vector of trait values is defined by y_i = y_iΛ + g_i + r_i, Λ being the (p + 1) × (p + 1) matrix of structural coefficients. The (k, l)th entry of Λ contains the effect of trait k on trait l, and the vectors g_i and r_i have zero mean Gaussian distributions with covariance matrices Σ^g and Σ^r, respectively. The joint distribution of all n(p + 1) trait values is then as in (1), with Σ^u = Γ^tΣ^gΓ and Σ^e = Γ^tΣ^rΓ, where Γ = (I − Λ)⁻¹ (Gianola and Sorensen, 2004; Töpner et al., 2017; Kruijer et al., 2020).

The target trait is defined as Y_f = Y₁ = λ(Y₂ + Y₃ + Y₄) + G₁ + R₁, and we do not assume any functional relations among the secondary traits. Hence, if λ ≠ 0, there is a causal effect from Y₂, Y₃, and Y₄ on Y₁, but the algorithms under consideration do not know which of the 300 secondary traits are the actual causal ones. We consider λ values on the grid {−1, −0.5, 0, 0.5, 1}. Σ^g has diagonal elements (0.2, 0.7, …, 0.7), i.e., the variances of the direct genetic effects are 0.2 for Y_f and 0.7 for each of the secondary traits. The off-diagonal elements corresponding to Y₁ vs. (Y₂, Y₃, Y₄) are ρG0.2·0.7, where we choose ρ_G ∈ {−0.5, 0, 0.5}. Similarly, Σ^r has diagonal elements 0.8 for Y_f and 0.3 for the secondary traits, and the off-diagonal elements between Y₁ and (Y₂, Y₃, Y₄) are ρE0.8·0.3, with ρ_E ∈ {−0.5, 0, 0.5}. The other off-diagonal elements in Σ^g and Σ^r are zero.

For the special case λ = 0 we have Γ = I, Σ^u = Σ^g and Σ^e = Σ^r, and Y_f will have a heritability of 0.2. The secondary traits will have heritability 0.7, and there is no causal effect of (Y₂, Y₃, Y₄) on Y₁. Genomic prediction for Y₁ can however still benefit from the genetic correlation between these traits (which is present when ρ_G ≠ 0). When λ ≠ 0, the causal effect of (Y₂ + Y₃ + Y₄) on Y₁ will introduce additional genetic and residual covariance in Σ^u and Σ^e.

For each of the 125 combinations of λ, ρ_G and ρ_E we simulate 50 data-sets; for each of them we predicted the simulated genetic effects for the test set, with the different methods.

To test the methods on real data, we consider four data-sets with various target and secondary phenotypes. To assess accuracy, each data set was randomly split into training (70%) and a test genotypes (30%). This was repeated 160 times, and we report accuracy averaged over the 160 test sets. Because of the required computing time, only 50 test sets were analyzed for RF-BLUP with hyper-parameter-optimization (for the Arabidopsis data-sets), and 30 test-sets for the maize data (for all methods). With one exception (mentioned below), the target and secondary phenotypes were measured on different plants; therefore, all bivariate mixed models were fitted with diagonal residual covariance (i.e., diagonal Σ^e in Equations 4 and 5).

The first two data sets were measured on the A. thaliana HapMap population, where 36 metabolites from Fusari et al. (2017) were used as secondary phenotypes and the kinship matrix was estimated based on one million imputed SNPs (Arouisse et al., 2020). Dataset 1 contains three target traits related to biotic and abiotic stress, from Thoen et al. (2017). In dataset 2, the target is the rosette fresh weight, measured in of the experiments of Fusari et al. (2017). This is the only dataset for which the residual covariance is non-diagonal.

In the third data set, we predicted the grain yield, plant height (PH) and flowering time (FT) of 388 inbred maize lines (Z. mays), using 5, 760 transcripts (Azodi et al., 2020) as secondary traits. In this case, we selected for each data-set a subset of transcripts using the LASSO on the training set, following Azodi et al. (2020). In other words, the transcripts selected by LS-BLUP were also used for the other methods.

The data that support the findings of this study are available at:

https://doi.org/10.1105/tpc.19.00332 (Maize data)

https://doi.org/10.1105/tpc.17.00232 (A. thaliana Metabolite data)

https://doi.org/10.1111/nph.14220 (A. thaliana Phenotypes)

https://doi.org/10.1111/tpj.14659 (A. thaliana SNP data)

All data-sets (except the maize transcriptomics) are included in an Rdata file available at: https://figshare.com/s/5d01062711ce33bb327e.

The required computing time is mainly driven by the complexity of fitting either a bivariate mixed model with a single relatedness matrix, or univariate mixed models with either one or two relatedness matrices. For the datasets considered here, each bivariate mixed model took between 20 and 50 s to fit, the univariate mixed models taking at most a few seconds. For complexity as function of n and p we refer to Zhou and Stephens (2014).

R-code for all methods is available at https://figshare.com/s/5d01062711ce33bb327e, where we mostly relied on asreml-R (Butler et al., 2009). Several open source alternatives are however available; in particular sommer (Covarrubias-Pazaran, 2016) for bivariate mixed models, and gaston for univariate mixed models. Using gaston's lmm.diago.likelihood function, the (univariate) GBLUP for large numbers of traits can be computed in only a few seconds, which is useful for the GM-BLUP method. For the dimension reduction in LS- and RF-BLUP we used the R-packages glmnet (Friedman et al., 2010), caret (https://cran.r-project.org/package=caret), and randomForest (Liaw and Wiener, 2002). For the maize data, LASSO and random-forest regression were performed in python, using the scikit-learn packages.

Figures 2, 3 show the estimated accuracy as function of λ, i.e., the size of the causal effects of Y₂, Y₃, and Y₄ on the target trait Y_f (i.e., Y₁). We focus on three cases, with different values for the correlations between the direct genetic effects on Y₁, …, Y₄, as well as the corresponding residuals (see section 2): (A) ρ_G = 0.5 and ρ_E = −0.5, (B) ρ_G = ρ_E = 0, and (C) ρ_G = 0.5 and ρ_E = 0.5. In scenario 1 (Figure 2) as well as scenario 2 (Figure 3), accuracies are generally higher when λ moves away from zero. This is expected, as the total genetic variance and heritability increase due to the causal effect, especially when ρ_G and λ have the same sign. When they have opposite sign, the lowest accuracy can occur at an intermediate value of λ [e.g., at λ = −0.5 in case of (A)].

The multi-trait benchmark with perfect information on the genetic and residual covariance between the target trait Y_f and secondary traits Y₂, Y₃, and Y₄ always outperforms univariate GBLUP, except when ρ_G = λ = 0, in which case accuracies are equal. When ρ_G ≠ 0, the benchmark always benefits from the genetic correlations between the target trait and the secondary traits, even if the latter do not have a causal effect on Y_f.

The accuracy of univariate GBLUP varied between r = 0.44 and r = 0.70, while the benchmark had accuracy between 0.50 − 0.70 (scenario 1) and 0.50 − 0.92 (scenario 2). The difference between scenario 2 (secondary traits observed on the test set) and scenario 1 (secondary traits only observed on the training set) was bigger for large values of |λ|. This is because for large |λ|, the total genetic correlation (which is also a function of ρ_G) between Y_f and the causal secondary traits (Y₂, Y₃, and Y₄) is larger.

In absence of a causal effect Y_s → Y_f (λ = 0) and residual genetic and residual correlations having opposite sign (case A), our simulation setup appeared to be too challenging, and none of the methods performed better than univariate GBLUP. Something similar occurred in case C, for λ = −0.5. On the positive side, for large values of |λ|, both SI-BLUP and LS-BLUP have near-benchmark accuracy, where the latter did not rely on plot-level observations. In scenario 2, RF-BLUP appeared to be an interesting alternative, with somewhat lower accuracy on the extreme sides, but relatively good performance at unfavorable values of λ.

Prediction based on the secondary traits only (M-BLUP; only available in scenario 2) is generally one of the least successful. The multi-kernel methods (Multi-BLUP and GM-BLUP) are somewhere in between, GM-BLUP often having an accuracy similar to that of RF-BLUP. GM-BLUP appears to be slightly better than Multi-BLUP, but in most cases the difference is smaller than the standard errors of the accuracy estimates.

Tables 1, 2 contain the accuracies for datasets 1–4 described above, averaged over randomly sampled test sets (see section 2). Because the original individual plant (or plot) data were not available, we could not compute the SI-BLUP here.

In scenario 1 (Table 1), none of the multi-trait methods performed consistently better than univariate GBLUP. For the second trait in data-set 1 (Salt5), RF-BLUP had accuracy 0.09, vs. 0.03 for univariate GBLUP; the latter had highest accuracy for the first and third trait in dataset 1 (fungus, and drought and fungus stress combined).

The remainder of this section we focus on scenario 2 (Table 2), in which there were more substantial differences among methods. For all datasets, methods based on multiple relatedness matrices (Multi-BLUP and GM-BLUP) had accuracies similar to single-trait GBLUP. As in the simulations, GM-BLUP gave only a minor (if any) improvement over Multi-BLUP. The approaches based on dimension reduction of the secondary traits (LS-BLUP and RF-BLUP) appeared to give a substantial improvement over univariate GBLUP, e.g., from r = 0.03 to r = 0.23 (LS-BLUP) for the Salt5 trait in data-set 1, or from r = 0.55 to r = 0.65 (RF-BLUP) for Maize yield in data-set 3, with transcriptomics as secondary traits.

LS-BLUP had the highest accuracy in all Arabidopsis datasets, with a small but consistent improvement over RF-BLUP (0.02–0.03 higher), also when optimized with the caret/scikit-learn packages. This hyperparameter optimization appeared to be rather important for the Maize data; using the default settings from the randomForest package (as in the simulations), accuracy was considerably lower (for yield and the transcripts for example, r = 0.65 vs. r = 0.51).

For the maize data, RF/LS-BLUP improved accuracy for yield from around 0.64 − 0.65 to 0.71 − 72 when plant height and flowering time were included as secondary phenotypes, together with the transcriptome data. None of the other methods could exploit the additional data, and accuracies were similar to those obtained with the transcripts alone. Prediction based on the secondary traits alone (M-BLUP) had around zero accuracy in all Arabidopsis data-sets, but r = 0.49−0.54 for the maize data, similar to GBLUP and multi-BLUP.