Our study used 250 rice accessions, including indica and japonica ecotype accessions from Asian countries (Supplementary Table S1). Accessions from China are mainly landrace indica varieties mostly cultivated on Cd-polluted soils in Guangdong province, China (Long et al., 2014).

The procedures followed for growing the 250 accessions and determining concentrations of trace elements in rice grain were as previously described by Liu S. et al. (2021). Briefly, seeds from the 250 accessions were first cultivated in pots filled with soil collected from the experimental station of the Agricultural Genomics Institute at Shenzhen, China. Next, germinated seeds were selected and cultivated in seedling trays for 4 weeks. Healthy seedlings were then transferred into pots containing soil amended with an initial concentration of Cd of 0.5 mg kg⁻¹. Finally, all the seedlings were planted in an augmented randomized complete block design with two replicates of 25 accessions from 20 July 2019 to 2 October 2019. Our study was limited to a single environment. Multi-environment data would be essential for understanding the environmental correlations and their stability and genotype effects by environment interactions (GxE). To determine grain concentrations of Zn, Fe, Cu, Mn, and Cd, grain samples were first peeled and dried at 65°C for 3 days. The dried samples were then crushed, wet-digested in concentrated nitric acid (HNO₃) at 120°C for 30 min, and further digested with perchloric acid (HClO₄) at 180°C until the samples became transparent. The samples were then diluted with ultrapure water. Finally, the grain concentration of each trace element was determined using the Inductively Coupled Plasma Mass Spectrometry (ICP-MS) machine (Houk et al., 1980).

The 250 accessions were genotyped following the re-sequencing and variants-calling procedures of the rice 3K project as reported by Wang et al. (2018). The following steps were implemented for all the genotypes to merge the variants-calling: First, raw reads were aligned to the R498 reference genome (Du et al., 2017) using the program bwa-mem alignment software (Li, 2013). Next, the PCR duplicates were identified with Picard software, version 2.9.0 (http://broadinstitute.github.io/picard/), and discarded. Following that, the GATK HaplotypeCaller engine (McKenna et al., 2010), with the option “-ERCGVCF,” was used to call genotypes at each site. The resulting genomic variants called format (gVCF) for each genotype were combined using the GATK Genotype GVCFs engine. Next, the GATK hard filter pipeline was used to individually call SNP and INDEL variants from the population variant file. All the variants within 5 bp of an INDEL were discarded. A variant was confirmed if at least one genotype supported it with a QUAL parameter greater than 30. After that, VCF tools indicated (Danecek et al., 2011) sites for which genotypes were not called in at least 20% of the used genotypes. The above procedures yielded 30,089,814 bi-allelic SNPs for the 250 genotypes. SNP quality control steps were implemented using PLINK software (Chang et al., 2015) with standards that remove SNPs with 1) minor allele frequency lower than 0.5, 2) call rate less than 0.9, and 3) pairwise LD (r ²) greater than 0.01. Finally, 36,171 SNPs were available for the 250 accessions.

The mixed linear model was used to estimate genomic heritability as follows: y = 1 n μ + g + e Where y is the vector of concentration of trace element under consideration, 1 n is the vector of ones, n is the number of genotyped cultivars, μ is the intercept, g is the vector of the genetic effects of accessions, e is the residual vector. g and e are assumed as random effects, respectively, following g ∼ N ( 0 , G σ g 2 ) and e ∼ N ( 0 , I σ e 2 ) where G is genomic relationship matrix estimated following (Yang et al., 2010) and I is identity matrix. The genomic heritability was estimated as h g 2 = σ g 2 ( σ g 2 + σ e 2 ) Where σ g 2 is the additive genetic variance component and σ e 2   is the residual variance component. GCTA software (Yang et al., 2011) was used to compute the genomic relationship matrix and genomic heritability.

The genomic correlation between traits was estimated using the formula: cor = co v i , j σ i σ j where c o v i , j is the genetic covariance between i^th trait and j^th trait, σ i and σ j are the square root of the genetic variance of i^th trait and j^th trait. The genetic covariances and variances were estimated using the R package MTM (De los Campos and Grüneberg, 2016).

A hierarchical cluster analysis based on the Euclidean’s distance matrix computed with the SNP genomic profiles was performed to inspect the genetic diversity among the 250 genotypes. In addition, a heat plot based on the cluster analysis was drawn to visualize the genetic dissimilarities.

The genomic prediction models used were the ST and MT models. The ST model only captured additive genetic effects, while the MT models accommodated both additive and local epistatic effects. ST model was the commonly used genomic best linear unbiased prediction (GBLUP) model and same as the mixed model estimating genomic heritability: y = 1 n μ + g + e , where y , 1 n , μ and ε were exactly as afore denoted, g is the vector of additive genetic effect in genotype-based model or additive plus local epistatic effects in the haplotype-based model. In the genotype-based model, we assumed g ∼ N ( 0 , G σ a 2 ) , where G is an n × n-dimensional additive genomic relationship matrix, and σ a 2 is the additive genetic variance component. In the haplotype-based model, we assumed g ∼ N ( 0 , H σ h 2 ) , where H is an n × n-dimensional haplotypic relationship matrix derived from the haplotypic profile matrix with values 0, 1, and 2 indicating the number of copies for a specific haplotypic allele (Jiang et al., 2018; He et al., 2019). To obtain the haplotypic profile matrix, the genotypic data of SNPs were phased using software SHAPEIT (Delaneau et al., 2012) with default augment settings. The phased genotypic data was recoded to haplotypic profiles using a fixed-length haplotype of 2, 3, 4, or 5 SNPs. G and H were established using software GCTA (Yang et al., 2010; Yang et al., 2011) based on the genotypic and haplotypic data, respectively.

For MT models, we used two approaches considering no correlation between traits; the Bayesian multi-output regressor stacking (BMORS) proposed by Montesinos-López et al. (2019); the MT-GBLUP; and two methods accommodating correlation between traits; the factorial analytic (FA) model and the unstructured (UN) model. For genotype-based approaches, the MT-GBLUP model was formulated as y = μ + u + ε , where y = ( y 1 ′ , y 2 ′ , … , y m ′ ) ′ , μ = ( μ 1 ′ , μ 2 ′ , … , μ m ′ ) ′ , u = ( g 1 ′ , g 2 ′ , … , g m ′ ) ′ , ε = ( e 1 ′ , e 2 ′ , … , e m ′ ) ′ , m is the number of traits included in the model. We assumed u ∼ N ( ( 0 ′ , 0 ′ , … , 0 ′ ) ′ , Ψ u ⊗ G ) , ε ∼ N ( ( 0 ′ , 0 ′ , … , 0 ′ ) ′ , Ψ ε ⊗ I ) where Ψ u = ( σ g 1 2 ⋯ 0 ⋯ 0 ⋮ ⋱ ⋮ ⋱ ⋮ 0 ⋯ σ g j 2 ⋯ 0 ⋮ ⋱ ⋮ ⋱ ⋮ 0 ⋯ 0 ⋯ σ g m 2 ) , Ψ ε = ( σ ε 1 2 ⋯ 0 ⋯ 0 ⋮ ⋱ ⋮ ⋱ ⋮ 0 ⋯ σ ε j 2 ⋯ 0 ⋮ ⋱ ⋮ ⋱ ⋮ 0 ⋯ 0 ⋯ σ ε m 2 ) , σ g j 2 and σ ε j 2 are respectively the genetic and residual variance of jth trait, and ⊗ denotes Kronecker product of matrices.

BMORS was a two-stage process. The first stage is the same as the MT-GBLUP, but instead of directly using the GBLUP predicted values as the final output, BMORS implemented a second stage that integrated the GBLUP predicted values of each trait in the first stage and fitted a ridge regression model. In this way, the prediction of a single trait could be corrected by the predictions of other traits in the first stage using the second stage model (Spyromitros-Xioufis et al., 2012; Spyromitros-Xioufis et al., 2016; Montesinos-López et al., 2019). The FA model was also based on the formula of the MT-GBLUP model but assuming a covariance structure between traits, that is Ψ u = ( σ g 1 2 ⋯ c o v g 1 g j ⋯ c o v g 1 g m ⋮ ⋱ ⋮ ⋱ ⋮ c o v g j g 1 ⋯ σ g j 2 ⋯ c o v g j g m ⋮ ⋱ ⋮ ⋱ ⋮ c o v g m g 1 ⋯ c o v g m g j ⋯ σ g m 2 ) = ( Λ Λ ′ + Π ) = F A ( k ) where k is the number of latent factors, Λ is a j × k dimensional matrix containing trait loadings, Π is a j × j diagonal matrix (Burgueño et al., 2012). Theoretically, the FA model requires at least three traits to be simultaneously included in the model. The UN model (Burgueño et al., 2012; Cuevas et al., 2017) tried to estimate all variances and covariances in Ψ u , i.e., σ g j 2 and c o v g i g j , i , j ∈ { 1 , … , m } , which may cause convergence problems when a large number of traits are considered. The haplotype-based approach was only implemented in the MT-UN model by replacing the relationship matrix G by H .

The ST-GBLUP model was implemented in R (R Core Team, 2016) using the BGLR package (De los Campos and Pérez-Rodríguez, 2015). The MT-GBLUP, FA, and UN approaches were realized using the R package MTM (De los Campos and Grüneberg 2016). BMORS was fitted using the R package BMTME (Montesinos-López et al., 2019). The number of iterations of all models was set to 20,000, and the first 12,000 were discarded as burn-in.

Four different CV schemes, referring to those reported by Lado et al. (2018) and Arojju et al. (2020) were used in our study (Table 1). CV1 was applied to both ST-GBLUP and MT models, referring to a scenario where the target trait was predicted without the support of auxiliary traits (ST-CV1) or with auxiliary traits only available in the training set (MT-CV1). CV2 and CV3 were only assessed for MT models. Under CV2 scheme, genotypes in both training and test sets had phenotypic data for all the auxiliary traits. Under CV3, phenotypes of the auxiliary traits were only available in the test set.

To assess the genomic prediction accuracy across the above CV schemes, the entire population of 250 genotypes was randomly divided into five equal-sized folds. Four folds collectively constituted the training set, and the remaining fold was the test set. Stochastic partitioning of training and test sets was repeated 20 times, yielding one hundred times (5 folds × 20 replicates) calibrations and predictions. The genomic prediction accuracy of the target trait was estimated using the Pearson correlation coefficient between the genomic predicted genetic values and the observed phenotypic values of 250 accessions when incorporated in the five test sets of each repeat of CV. The Student t-test was used to test the statistical difference in genomic prediction accuracies among the prediction models.

To investigate whether the predictive ability of MT genomic prediction on rice grain trace elements concentration could be boosted by optimizing the SNP marker panel, we applied a (GWAS) to identify the trait-associated SNPs (TA-SNPs) and establish the selective marker panel. The CV scenarios in which the MT-UN haplotype-based models disregarding the length of haplotypes (two to five SNPs) constantly showed statistically significantly (p < 0.05, t-test) higher prediction accuracies than their genotype-based counterparts and the ST-GBLUP model were used to validate the efficacy of using the TA-SNPs to train the genomic prediction models. In more detail, GWAS using the total SNP marker panel was performed in the training set of the designated CV scenarios. First, the TA-SNPs with p values less than 0.01 were recorded. Then each chromosome was divided into bins spanning 300 kb (the bin size is decided by the LD decay, with the physical distance between pairs of SNPs based on the total population). Finally, the most trait-associated SNPs with the lowest p-value in each bin was picked together with the TA-SNPs (p < 0.01) and recorded to constitute the selective SNP marker panel of each repeat of CV. The MT-UN genotype-based model was implemented using the genotypic data of the selective SNPs. Contrastingly, the adjacent selective SNPs located within 300 kb (highly possible as the position of the selective SNP from each bin is unfixed) were combined to compile the haplotypes using the phased genotypic data as the LD decay implied a non-negligible LD among them. The remaining SNPs without close neighbors within 300 kb were maintained, and their genotypic profiles were used. Therefore, the MT-UN haplotype-based model took advantage of both haplotypes and genotypes. The GWAS was implemented using GCTA software (Yang et al., 2011; Yang et al., 2014). The additive genomic relationship matrix was exclusively used in the GWAS model to account for the relatedness between accessions.

The LD decay distance between all SNP markers for the 250 accessions was ∼250–300 kb when the cut-off value (r ²) was set at 0.1, assuming non-negligible SNP pairwise correlation (r = 0.3) (Figure 1A). The kinship between accessions was determined based on pairwise Euclidean distances. Pairwise distances among accessions ranging from 0 to 0.2 accounted for less than 5% of all the pairwise distances. Pairwise distances from 0.6 to 0.8 were the most frequent and accounted for 12%–23% of all the pairwise distances (Figure 1B). Further, no genetically structured sub-populations were observed among the 250 varieties used in this study. However, several families were detected (Supplementary Figure S1).

The distributions of phenotypes (adjusted phenotypic means) based on the 250 accessions varied among the five traits studied. The distribution of Zn was almost symmetrical. The skewness was high and negative for Cu, and moderate and negative for Fe, Cd, and Mn (Supplementary Figure S2). The genomic heritability for all studied traits ranged from low (Zn: 0.14 and Cu: 0.21) to medium (Mn: 0.35) and high (Fe: 0.5 and Cd: 0.62) (Table 2). The genetic correlation estimated with the MTM model was highest between Fe and Cd (0.95) and Cu and Zn (0.95) and lowest between Mn and Cd (0.39) and Mn and Fe (0.44). Zn had the highest genetic correlations with all the other studied traits, ranging from 0.67 to 0.95 (Table 2).

The average prediction accuracy with the traditional ST-GBLUP model under the CV1 scheme was the highest for Cd (0.52), followed by Fe (0.39), Mn (0.36), Zn (0.23), and Cu (0.21) (Figure 2). Also, under the CV1 scheme, prediction accuracies of MT models were not statistically significantly (p < 0.05) superior to those of ST-GBLUP irrespective of the models and traits studied (Supplementary Tables S2–S6).

As compared, when phenotypes of the auxiliary traits were made available in both training and test sets (MT-CV2) or merely in the test set (MT-CV3), the MT models, namely FA or UN, significantly (p < 0.05) outperformed the ST-GBLUP model (Figures 3, 4). For most of the studied traits, the highest performance of MT models was observed under the MT-CV2 scheme ( Figures 3, 4 ). However, MT-GBLUP and BMORS MT models were not significantly (p < 0.05) superior to the ST-GBLUP model for all the CV schemes studied (Supplementary Tables S7–S11).

We further compared scenarios where the prediction of the target traits was assisted with a single auxiliary trait or a combination of multiple auxiliary traits in MT models. Supporting the prediction of Zn with one of its correlated traits (Cu, Fe, Mn, or Cd) was sufficient to significantly (p < 0.05) increase the prediction accuracy MT-UN model relative to ST-GBLUP model in MT-CV2 (Figure 3A). Cu was the best single auxiliary trait for predicting Zn. Incorporating observations from Cu in MT-UN model (under MT-CV2) significantly (p < 0.05) increased the prediction accuracy of Zn by 82.6% (0.23–0.42) relative to the ST-GBLUP model. However, the highest increase in prediction accuracy (126% or 0.23–0.52) of MT models was observed when observations from Mn, Fe, and Cu were combined as supporting traits for Zn under the MT-CV2 scheme (Figure 3A). Under MT-CV3, the MT-UN model outperformed the ST-GBLUP model only after multiple auxiliary traits were used to support the prediction of Zn (Figure 3A).

Similarly, compared to ST-GBLUP, the prediction accuracy of Cu by the MT-UN model significantly (p < 0.05) increased by 95.2% (0.21–0.41) and 38% (0.21–0.29) in MT-CV2 and MT-CV3, respectively, when Zn was used as a single supporting trait (Figure 3B). Yet, after including other traits in MT-CV2 (Mn and Zn) and MT-CV3 (Zn and Cd), the prediction accuracy improved by 109.5% (0.21–0.44) and 57.1% (0.21–0.33) relative to ST-GBLUP, respectively.

Similarly, when Fe was the target trait, MT-UN model accounting information from Cd significantly (p < 0.05) outperformed ST-GBLUP by 23% (0.39–0.48) and 12.8% (0.39–0.44) under MT-CV2 and MT-CV3 schemes, respectively (Figure 3C). Nevertheless, considering phenotypes from more auxiliary traits in MT-CV2 (Zn and Cd) and MT-CV3 (Zn and Cd) provided 33.3% (0.39–0.52) and 23% (0.39–0.48) gains in the prediction accuracy of Fe with the MT-UN model, respectively (Figure 3C).

Furthermore, the prediction accuracy of Cd (the most heritable trait) with the MT-UN model was significantly (p < 0.05) improved by 3.8% (0.52–0.54) when Fe, its strongly correlated trait, was used as a single auxiliary trait under the MT-CV2 scheme (Figure 4A). As observed for the other traits, 7.6% (0.52–0.56) and 11.5% (0.52–0.58) gains in prediction accuracy were attained after using combined information from multiple auxiliary traits (Mn, Fe, and Cu) in MT-UN and MT-FA models under MT-CV2, respectively (Figure 4A). Similarly, under MT-CV3, MT models did not significantly outperform ST-GBLUP models in scenarios where a single auxiliary trait was used. However, when information from Fe and Zn or Fe and Cu was considered, an improvement of 5.7% (0.52–0.55) in the prediction accuracy of the MT-FA model over the ST-GBLUP model was observed (Figure 4A).

Finally, when Mn was the target trait, MT-UN with a single auxiliary trait failed to improve its prediction accuracy under both MT-CV2 and MT-CV3 schemes (Figure 4B). However, considering information from additional traits (Cu, Zn, and Cd), using the MT-FA model significantly (p < 0.05) improved the prediction accuracy of Mn up to 11.1% (0.36–0.40) over ST-GBLUP under MT-CV2 (Figure 4B). On the other hand, under the MT-CV3 scheme, MT-UN or MT-FA did not significantly outperform ST-BLUP even after using multiple auxiliary traits.

We further investigated the benefits of accommodating local epistatic effects on the prediction accuracy of MT models by using haplotypes instead of genotypes in the UN model. Comparing to the genotype-based UN model, the observed largest and significant (p < 0.05) increment of prediction accuracies using haplotype-based models was 3.8% for Zn (0.52–0.54), 4.6% for Cu (0.43–0.45), and 3.5% (0.56–0.58) for Cd under MT-CV2. For Zn, the above improvement in prediction accuracy was achieved with a haplotype length of 3 SNPs, and when Mn, Fe, Cu, and Cd were collectively used as auxiliary traits (Figure 5A). For Cu, the observed gains were realized with a haplotype length of 4 SNPs and when auxiliary traits Mn and Zn were used together (Figure 5B). For Cd, the gains were from the MT-UN model with a haplotype length of 2 SNPs, and when Mn, Fe, Cu, and Zn were combined as the auxiliary traits (Figure 5C). Under MT-CV3, the haplotype-based UN model was significantly (p < 0.05) superior to the genotype-based UN model by 12.5% (0.32–0.36) for Zn (Figure 5A) and 6% (0.33–0.35) for Cu (Figure 5B).

Compared to the ST-GBLUP model, the haplotype-based UN models were significantly (p < 0.05) superior with an increment of prediction accuracy of 134.7% (0.23–0.54), 114.2% (0.21–0.45), 23% (0.39–0.48), and 11.5% (0.52–0.58) for Zn (Figure 5A), Cu (Figure 5B), Fe (Supplementary Figure S3), and Cd (Figure 5C), respectively.

With this study, we sought to investigate whether the prediction accuracy of Zn, Cu, and Cd with the haplotype-based model can be improved by using the selective marker panel made by TA-SNPs derived from GWAS (Supplementary Table S12). We purposely selected Zn and Cu because for both traits the haplotype-based MT model performed superiorly for several scenarios irrespective of the lengths of haplotypes (2–5 SNPs) (Figure 5). We also investigated Cd in addition to Zn and Cu since successive significant (p < 0.2) TA-SNPs were observed in GWAS based on the total population (Supplementary Figures 4A–C). The non-negligible LD (r ²≥0.1) observed between the TA-SNPs, especially in Cu and Cd, underpinned the necessity of modelling local epistatic effects among TA-SNPs (Figures 6A–C). The haplotype-based UN model accounting for TA-SNPs significantly (p < 0.05) outperformed their genotype-based counterparts; however, it was significantly (p < 0.05) inferior to the model using all genome-wide markers for all the traits and scenarios evaluated (Figures 7A–C).

In this study, MT models did not significantly outperform ST-GBLUP under the CV1 scheme for all the scenarios evaluated (Supplementary Tables S2–S6). Earlier studies also reported insignificant differences in the prediction accuracies of MT-CV1 and ST-CV1 (Calus and Veerkamp, 2011; dos Santos et al., 2016; Bhatta et al., 2020), implying that MT models are not always robust over ST models, especially when information on auxiliary traits is only available in the training set and the unobserved accessions are predicted only based on genotypic data. In contrast, when phenotypes of the auxiliary traits were present in the training and test set (MT-CV2) or merely in the test set (MT-CV3), the prediction accuracy of MT models (MT-UN and MT-FA) for the unobserved target traits (Zn, Cu, Fe, and Cd) was significantly improved relative to ST-GBLUP (Figures 3, 4). Previous studies attributed the predictive performance of MT models to both higher heritability of the auxiliary trait and strong genetic correlation between the target and auxiliary traits (Sun et al., 2017; Fernandes et al., 2018).

Using a single auxiliary trait in the MT-UN model significantly (p < 0.05) improved the prediction accuracy of target traits Zn, Cu, Fe, and Cd relative to the ST-GBLUP model (Figures 3, 4). When a strong genetic correlation exists between target and auxiliary trait, the prediction accuracy of MT models could still be improved under MT-CV2 or MT-CV3 regardless of trait heritability. For instance, supporting the prediction of Cu with Zn, its strongly correlated trait (cor _{Zn, Cu} = 0.95) but with lower heritability (h² _Zn = 0.14; h² _Cu = 0.21) significantly improved the prediction accuracy of Cu with the MT-UN model (Figure 3B). Also, supporting Cd with Fe, its strongly correlated trait (cor _{Cd, Fe} = 0.95) but with lower heritability (h² _Cd = 0.62; h² _Fe = 0.50), improved the prediction accuracy of Cd with the MT-UN model (Figure 4A). Arojju et al. (2020) also indicated that the genetic correlation was the main cause of the observed gain in prediction accuracy of MT models. The same study further showed that when a trait in strong genetic correlation with the target trait is used in the MT model, the predictive performance of the MT model was still superior to the ST model even after reducing the training population size by 50%.

Collectively accounting for phenotypes of multiple auxiliary traits further improved the predictive ability of the MT models compared to the MT models with a single auxiliary trait. For example, the highest increase in the prediction of Zn was 82.6% when a single auxiliary trait was used in MT-UN models. Yet, using multiple traits collectively in the same model improved the prediction of Zn by 126% compared to the ST-GBLUP model (Figure 3A). Also, MT models with one auxiliary trait showed no benefit over the ST-GBLUP model when predicting Mn, with relative medium heritability and no strong genetic correlation with any other studied trait. However, when auxiliary traits were collectively used in the MT-FA model, significant improvements in the prediction accuracy over the ST-GBLUP model were observed (Figure 4B). Multiple auxiliary traits would optimize MT models, though the assisting trait per se is neither strongly genetically correlated with the target trait nor highly heritable. Therefore, when no single auxiliary trait meets the criterion of heritability or genetic correlation, combining multiple auxiliary traits in the MT model could be an effective approach to enhance the predictive ability of MT models. These findings are concurrent with previous findings by Wang et al. (2017), indicating that the prediction accuracy of MT models was highest when eight different traits were used as auxiliary traits to predict grain yield in rice.

Previous studies demonstrated that accounting for local epistatic effects besides the additive effect in genomic prediction could improve the prediction accuracy of ST models (Akdemir and Jannink, 2015; Akdemir et al., 2017; Jiang et al., 2018; He et al., 2019). Here, we are the first to attempt to model the local epistatic effects in the context of MT genomic prediction. Accounting for the local epistatic effects in haplotype-based MT models significantly improved the prediction accuracy of Zn, Cu, and Cd relative to genotype-based MT models, only capturing additive effects (Figures 5A–C). Relative to ST-GLUP, the highest increase in prediction accuracy, 134.7% for Zn, was observed after incorporating the local epistatic effects into the MT-UN model (Figure 5A). These findings imply that the potential of MT models can be maximized by accounting for local epistatic effects besides additive effects in the model.

Using a selective marker panel based on approaches exploiting the trait biological and genetic background knowledge such as GWAS has been proven effective to improve the predictive ability of GS models (Owens et al., 2014; Wang et al., 2019). Our study did not show any improvements of prediction accuracy by using the TA-SNPs instead of all genome-wide SNPs for Cu, Zn, and Cd (Figures 7A–C). These findings could be attributed to the complex genetic architectures of the trace elements we studied (Supplementary Figures S4A–C). Our approach was slightly similar to previous methods using GWAS-derived TA-SNPs to construct the trait-associated matrix (Zhang et al., 2014; Ahmadi et al., 2021), except that we did not assign weights to haplotype- or genotype-based genomic relationship matrices. Though numerous studies reported improved gains from using the above strategy (Bhandari et al., 2019; Ahmadi et al., 2021), Veerkamp et al. (2016) showed that the proportion of total variance explained by the TA-SNPs combined in a GRM was considerably smaller than that explained by all variants in Holstein-Friesian cattle population. A potentially more promising way to use TA-SNPs would be to fit them as fixed effects in the GP model along with all other SNPs as random effects (Spindel et al., 2016). However, the latter approach is best suited for features with a few large-effect QTLs in a polygenic context (Poland and Rutkoski, 2016; Bian and Holland, 2017; Rice and Lipka, 2019). Therefore, the genetic architecture of the target traits must be studied before applying this strategy to a breeding program. The marked advantage of the haplotype-based UN model over their genotype-based counterparts using the TA-SNPs (Figures 7A–C) substantiates the existence of local epistasis in trace element traits (Sharma V. et al., 2021) and the merit of modelling local epistatic effects in MT-GP program.

Various factors affect the predictive ability of GP models used in GS (Crossa et al., 2017; Xu et al., 2021). In this study, the genomic heritabilities spanned a wide range from 0.14 to 0.62 (Table 2), which enabled the evaluation of the performance of MT models under contrasting levels of genomic heritability. The genomic heritability of Zn was the lowest (h² = 0.14, Table 2), which contradicts several previous studies that reported moderate to high heritability of Zn (Norton et al., 2010; Pinson et al., 2015; Naik et al., 2020). The poor heritability estimate of Zn in this study could be due to potential environmental effects. Unfortunately, our study does not include multi-environment trials and therefore does not provide insight into environmental factors and GxE interactions on genomic prediction of trace element traits.

Strong genetic correlations (Cor > 0.75) were observed between pairs of the studied trace element traits (Table 2). This was expected due to their overlapping genetic and physiological mechanisms (Sasaki et al., 2012; Cu et al., 2020). For example, transporter gene families like zinc-iron permease (ZIP), natural resistance-associated macrophage proteins (NRAMPs), and heavy metal transporting, ATPases (HMAs) have been associated with uptake and translocation of several trace elements in plants (Fernández-Paz et al., 2021; Vanderschueren et al., 2021; Zhang et al., 2022). As a result, borrowing information from correlated traits overall improved the prediction accuracy of the MT models.

The success of GS is also highly dependent on the LD between markers and unknown causal variants. The genetically distant training and test sets would have different LD decay patterns and consequently impede the prediction (Snelling et al., 2013; Desta and Ortiz, 2014; Thistlethwaite et al., 2020). Such a problem is typically prominent in germplasm accessions, limiting the power of GP (Crossa et al., 2017). This is also the case for the diverse rice population used in this study, as most accession pairs were distantly unrelated (Figure 1B). As a result, the predictive ability for ST-GBLUP was poor (Figure 2), particularly for traits with the lowest estimated heritabilities, namely Zn and Cu (Table 2). Adding related materials to the training population has been suggested to overcome the problem of low relatedness between training and test populations and improve the accuracy of genomic prediction (Arenas et al., 2021). Nevertheless, increasing relatedness will damage genetic gain in the long term because genetic variation will be limited or exhausted if related populations are overused (Jannink et al., 2010; Moeinizade et al., 2019). We show that MT models are powerful tools for predicting trace element traits in populations with diverse backgrounds. However, further studies with larger datasets are needed to elucidate the utility of different populations and marker optimization strategies in the context of MT genomic prediction.

Expediting genomic selection in gene banks to predict the genetic merit of the unobserved accessions would enable accurate identification of promising donor accessions without a comprehensive phenotypic test of all the accessions in the field (Pace et al., 2015; Yu et al., 2016; Crossa et al., 2017; Tanaka et al., 2021; He et al., 2022). In fact, as the traits of breeders’ interest are extensive, the genetic resources archived in gene banks would be evaluated for several traits. MT genomic prediction is an effective method to realize this comprehensive evaluation. Our study tested several MT models under three different CV schemes, implying different phenotyping layouts and costs. We found that high prediction accuracy of MT models can be achieved under prediction schemes MT-CV2 and MT-CV3. The MT-CV2 scheme requires more budget for phenotyping auxiliary traits in both training and test sets. Therefore, breeders may kindly MT-CV2 if phenotypes for the auxiliary traits can be inexpensively obtained. Otherwise, MT-CV3 is more cost-effective as phenotypes for the auxiliary traits are only required for the test set (i.e., 20% of the entire population). Since using multiple auxiliary traits collectively in the MT model can improve prediction even if the individual auxiliary traits do not fully meet the heritability and genetic correlation conditions, an ideal situation would be to phenotype less expensive and more manageable traits (e.g., root system architecture, 100-grain weight, data to heading, etc.) to support the prediction of expensive target traits with the MT models. Besides, accounting for local epistatic effects in MT models would help to improve the predictive ability. The different scenarios studied here and their respective potentials in terms of prediction accuracy and phenotyping cost are illustrated in Figures 8A,B.

To breed safe and nutritious crop varieties, further studies using the genomic selection index (Habyarimana et al., 2020), for example, are desired to provide a comprehensive understanding of the strategies to optimize essential nutrients and toxic metals such as Cd in food crops.