A panel of 186 durum wheat genotypes was used for this study (Supplementary Table 1). This panel has already been employed in a previous study published by Puglisi et al. (2025) to implement genomic prediction models for morpho-phenological traits [plant height (PH, cm); heading date (HD) and flowering time (FT), expressed as growing degree days (gdd)]. The same study also describes the genotyping procedures used for this panel, providing a comprehensive genomic characterization of the material.

The genotypes were grown across a few simulated Mediterranean conditions, including three sowing dates (early, optimal, and delayed) over three consecutive growing seasons (from 2021–22 to 2023-24) at the CREA Research Centre for Cereal and Industrial Crops, Foggia, Italy (41° 51′ N, 15° 52′ E, 60 m elevation).

In addition to morpho-phenological traits (i.e., HD, FT, and PH) previously reported (Puglisi et al., 2025), at the ripening stage (Zadoks stage 90) (Zadoks et al., 1974), the entire panel was manually harvested to evaluate yield-related traits for each sowing-by-season combination. The grain number per spike (GN), grain weight per spike (GW), number of spikelets per spike (NS), spike length (SL), and spike weight (SW) were measured in 25 randomly selected spikes from each plot across all sowing dates and growing seasons.

Adjusted means (Best linear unbiased estimates, BLUEs) were computed from raw measurements using the following linear mixed model:

where ysjr is the response variable, µ is the general mean, Gens is the fixed effect of the sth genotype and follow a normal distribution with mean 0 and variance equals to σg2, that is Gens∼NIID(0,σg2), the Repj is effect for the jth replication, Blockr(Repj) is the random interaction effect of rth block within the jth replication, finally ϵsjr is the residual terms normally distributed with mean 0 and variance equals to σϵ2, that is ϵsjr ∼NIID(0,σϵ2). This model (Equation 1) was implemented using the lme4 package (Bates et al., 2015) in R 4.0.3 statistical environment (R Core Team, 2020).

Pairwise correlations among traits were calculated for each environment and visualized using two complementary approaches in R 4.0.3 statistical environment (R Core Team, 2020): i) relationships among traits across all environments were explored using the GGally package (Schloerke et al., 2025) via the ggpairs function; ii) for each trait, correlations between environments were displayed as heatmaps using a custom ggplot2-based function (Wickham, 2016). The function reordered rows and columns based on hierarchical clustering.

Site regression (SREG) analysis was used to obtain more information on genotype-by-environment (GxE) interactions. SREG analysis was performed according to Samonte et al. (2005) using imputed values in GEA-R version 4.1 (Samonte et al., 2005; Pacheco et al., 2018), and scores of genotypes were plotted on a biplot.

Genomic Best Linear Unbiased Prediction (GBLUP) is one of the foundational models in modern genomic selection. It extends the classical Best Linear Unbiased Prediction (BLUP) framework by replacing the pedigree-based additive relationship matrix (A) with a genomic relationship matrix (G) derived from molecular marker data. This allows breeders to predict the genetic merit of unphenotyped individuals using genotypic information from dense genome-wide markers. GBLUP assumes that all marker effects contribute equally and independently to the total genetic variance, and it provides Best Linear Unbiased Predictions (BLUPs) of genomic breeding values (GEBVs).

Cross-validation (CV) is a conventional method for assessing predictive performance in genomic prediction studies. In genomic prediction, CV1 and CV2 refer to two common CV strategies used to evaluate genomic model performance across distinct prediction scenarios. These schemes simulate realistic breeding conditions and make predictions for either new genotypes or new environments. Here, we conducted 10 cycles of 5-fold random cross-validation to assess model predictability under SE, MT, ME, and MTME. The strategies MT, ME, and MTME were performed using both CV1 and CV2.

The CV for SE was performed within each year-sowing date combination separately. This consisted of randomly partitioning the population into five sets, with four used to train the model and the remaining set used to test it. This process is usually performed until all sets have been used for testing at least once. Then, this process was repeated for 10 cycles.

Additionally, both CV1 and CV2 are well explained by Lado et al. (2018) and Jarquín et al. (2014) for multi-trait and multi-environment, respectively. However, we provided brief descriptions of both cross-validation methods in the next paragraphs.

CV1 – Prediction of Unseen Genotypes (within known environments and traits).

The goal of CV1 is to assess how well the model predicts new lines (genotypes) that have never been phenotyped but are genotyped, assuming that the environments and/or traits are already represented in the training data. In CV1, the dataset is divided by genotypes: some are used for training and others (never seen during training) for testing. All environments are represented in both subsets. This setup mimics a “new candidate line prediction” scenario, where breeders aim to predict the performance of untested genotypes before field evaluation. CV1 primarily tests the model’s ability to generalize across genotypes, relying on molecular marker information and genetic relationships. Prediction accuracy in CV1 reflects how effectively genomic data alone can capture genetic value when environmental effects are already known.

CV2 – Prediction of tested genotypes (within known environments and traits).

The goal of CV2 is to evaluate how well a model predicts the performance of known genotypes in new environments (evaluated for other genotypes) for the primary (target) trait. In this setting, target genotypes are observed during training in some environments and/or for both primary and secondary traits, while their phenotypic records for the primary trait are missing in the target environment. Information from other genotypes evaluated in the target environment, as well as from secondary traits measured on the target genotypes, is available and used for model training. In CV2, the dataset is divided by environments and/or traits rather than by genotypes. The same set of genotypes is included in both training and testing, but phenotypic records from some environments are omitted during training. In this cross-validation strategy, the aim is to predict the performance of known genotypes in environments and/or traits for which they have no observed phenotypes.

Adjusted means (BLUEs) of the phenotypic data collected over three consecutive growing seasons (2021–2024) and multiple sowing-by-season combinations (environments) reveal substantial variability across traits and conditions (see diagonal plots of Figure 1; Supplementary Table 2).

The phenotypic distribution of HD and PH were strongly influenced by sowing date and growing season, respectively, as previously reported (Puglisi et al., 2025). Similarly, a wide range of variation was observed across environments for yield-related traits, indicating the influence of both genetic and environmental factors.

GN showed the widest phenotypic range, with values ranging from 9.89 to 83.65 grains, and the coefficient of variation ranging between 0.154 and 0.251. The broad-sense heritability (h²) across environments varied from 0.352 to 0.553, suggesting moderate genetic control.

GW ranged from 0.343 to 3.518 g, with coefficient of variation ranging from 0.132 to 0.282. The heritability values were moderate (h² = 0.302–0.533), reflecting a stable genetic influence under different environmental conditions.

NS ranged from 11.35 to 28.33, showing relatively low coefficient of variations (0.070–0.113), suggesting more stable expression across environments. Heritability values were consistently high, up to 0.642, pointing to a stronger genetic determinism.

SL values were between 4.765 and 14.555 cm, with coefficient of variation ranging from 0.127 to 0.162 and high heritability (h² = 0.674–0.736) across all evaluated environments, indicating that this trait is under strong genetic control.

SW showed a wide range (0.669 to 4.880 g) and moderate coefficient of variation (0.128–0.236). Heritability values ranged from 0.355 to 0.540, again reflecting moderate genetic contribution.

In general, traits such as SL and NS exhibited higher heritability, suggesting they are less affected by environmental variation and may respond better to selection. Conversely, GN and GW, while showing significant phenotypic variability, presented slightly lower heritability estimates, likely due to stronger genotype-by-environment interactions.

The prediction of morpho-phenological (PH and HD) achieved exceptionally high accuracies (~0.9), particularly when using MTME models (Supplementary Figures 2A, 3A; Supplementary Tables 8, 9).

Regarding the HD prediction, the MTME_CV2 structure consistently delivered the highest predictive abilities, ranging from 0.90 to 0.97 with an average of 0.95 across all sowing-by-season combinations (Supplementary Figure 2B). The MTME_CV2 models once again dominated the top rankings, occupying the majority of the Top-20 positions (Supplementary Figure 2C; Supplementary Table 8). The best-performing case was MTME_CV2 for environment 2_2022–2023 (r = 0.9669), followed closely by MTME_CV2 applied to 2_2023–2024 (r = 0.9654) and 3_2022–2023 (r = 0.9594). These high accuracies, all above 0.95, confirm the remarkable predictive performance and reproducibility of MTME structures for HD.

These strong results indicate that HD, as a highly heritable phenological trait, can be predicted with very high accuracy using models that integrate both trait and environmental covariation. Multi-environment (ME_CV2_G2) models followed closely, achieving r = 0.92–0.94 across several environments (2_2023–2024, 2_2022–2023, 1_2022–2023, and 3_2022–2023). These results highlight the advantage of including allelic information (G2) to improve cross-environment prediction. In contrast, MT and SE models achieved lower but still strong accuracies (r = 0.81–0.82), indicating reliable prediction even under simplified structures, maintained solid performance levels, consistent with the high heritability of HD. Across all Top-20 models, predictive ability ranged from 0.97 to 0.81, with most cases exceeding 0.90, confirming HD as one of the most accurately predicted traits in the study, reflecting the exceptional stability and robustness of MTME_CV2 and ME_CV2_G2 models. The overall robustness and precision of MTME_CV2 highlights its effectiveness in predicting phenological traits such as heading date in complex breeding datasets.

Conversely, the prediction of the morphological trait (PH) performed remarkably well for all models, with the MTME model standing out by achieving the highest predictive accuracies (Supplementary Figure 3B). The MTME_CV2 framework was the dominant performer, reaching r = 0.9723 in environment 2_2022–2023, followed closely by MTME_CV2 models for 2_2021–2022 (r = 0.9698) and 1_2022–2023 (r = 0.9688). These results, all above 0.96, indicate exceptional predictive reliability and confirm that PH is a highly heritable and stable trait across environments (Supplementary Figure 3C; Supplementary Table 9). ME_CV2_G2 models followed closely behind, typically achieving r = 0.94–0.95 across multiple environments. Their strong performance demonstrates that the integration of allelic combinations (G2) adds genetic resolution and robustness to multi-environment predictions. Standard deviations remained low (0.01–0.03), confirming model stability across CV folds.

These results highlight the consistency of both MTME_CV2 and ME_CV2_G2 frameworks in capturing G×E patterns across sowing-by-season combinations. MTME and ME models clearly outperformed MT and SE approaches, even for instance MT_CV1_G2 and SE_G2 models yielded slightly lower but still strong accuracies (r = 0.89–0.91), showing that accounting for both genetic and environmental structure is critical for accurate prediction of PH. Simpler models such as MT_CV1_G2 and SE_G2 remained competitive, confirming that PH prediction is robust even under reduced model complexity.

Overall, these findings validate MTME_CV2 as the most effective and consistent approach for predicting both morpho-phenological traits across diverse sowing-by-season combinations and environmental conditions.

Across all traits evaluated, MTME models consistently yielded the highest prediction accuracies, with performance gains of 5–20% over MT or ME models, and up to 40% relative to the SE baseline. The MTME framework leverages both genetic and environmental covariances, effectively borrowing information from correlated traits and environments (Montesinos-López et al., 2016, 2018b, 2018c, 2019a; Crossa et al., 2025a, 2025b). Such integrative modeling captures shared genetic determinants influencing multiple traits while accounting for plastic responses to environmental cues, thus providing a biologically realistic representation of durum wheat adaptation.

The improvement obtained under MTME-CV2 relative to CV1 reflects the model’s robustness in predicting genotype performance in unobserved environments for target lines, a critical feature for forward breeding and climate-resilient selection (Lopez-Cruz et al., 2015; Cuevas et al., 2018). Similar patterns have been reported in spring wheat (Sukumaran et al., 2017), maize (Jarquin et al., 2020), and barley (Bhatta et al., 2020), confirming that the joint exploitation of inter-trait and inter-environment covariance is a universally effective strategy across cereals. In durum wheat, where environments differ substantially between years and sites (e.g., Canada, Italy, Turkey, Mexico), these cross-links become even more relevant. Additionally, MT and MTME showed higher predictive performances, especially in CV2 designing. This is mostly due to the validation design used, which allows the inclusion of the performances of the target lines in some environments and of secondary and target traits in the training population, and predicts their target traits in other environments (testing), where the secondary traits are also available. This has widely been observed in previous studies (Lado et al., 2018; Vitale et al., 2024).

Finally, it is worth noting that while SE and ME analyses can be efficiently run on most commercial laptops, MT and, in particular, MTME analyses may be computationally demanding for such machines. Using a high-performance computing system, we observed a substantial increase in computational time for MT and MTME models compared with SE and ME, which were considerably faster.

The magnitude of prediction accuracy varied among traits, reflecting their genetic architecture, heritability, and correlation with other traits.

High predictive ability (r > 0.81) with MTME-CV2 was observed for all five yield-related traits investigated. This suggests that the MTME model effectively disentangled additive and non-additive sources of variation across seasonal fluctuations in temperature and water availability. The simultaneous modeling likely benefited from their shared physiological pathways in sink-source allocation, as previously shown in multi-trait genomic studies in bread and durum wheat (Lado et al., 2018; Montesinos-López et al., 2019b). Their predictable behavior within the MTME framework is consistent with the genetic stability of these traits and their moderate-to-high heritability (Taranto et al., 2023).

Importantly, heading date (HD) and plant height (PH) reached the highest accuracies overall (r > 0.95), confirming their simpler genetic control and robust expression across environments. These traits are highly correlated with flowering time and photoperiod sensitivity genes (e.g., Ppd-A1, Vrn-A1), whose allelic effects were captured in the G2 matrix incorporating known functional variants (Puglisi et al., 2025).

Taken together, these results confirm that combining multiple traits within an MTME framework not only increases predictive accuracy but also improves the interpretability of genotype responses across time and space. In durum wheat, where yield is the emergent outcome of multiple interacting developmental traits, such multi-trait modeling mirrors the true biological complexity of performance.

An innovative aspect of this study is the inclusion of both the conventional genomic relationship matrix (G) and an allelic matrix (G2) capturing functional gene combinations. This dual approach allowed quantification of how allelic information at key loci complements genome-wide additive effects. As previously reported for the SE and ME GP models (Puglisi et al., 2025), the G2-based models improved the prediction accuracy for morpho-phenological traits by 2–4%, also when using MT and MEME approaches, suggesting that known allelic variants contribute unique information not captured by SNP-based relationships. This is consistent with recent reports integrating gene-level annotations and known regulatory regions into genomic prediction frameworks (Crossa et al., 2025a).

Moreover, although the improvement was modest compared to morpho-phenological traits, the use of G2 increased the prediction accuracy for yield-related traits. This outcome supports its stabilizing effect under multi-environment conditions and highlights the importance of jointly modeling allelic effects, environmental factors, and trait correlations. The limited improvement may reflect the need to further investigate genes with stronger effects on these traits (Nadolska-Orczyk et al., 2017).

In durum wheat, functional markers linked to Ppd, Rht, and Vrn genes have well-documented effects on developmental timing and architecture. Incorporating these alleles into GP models enhances prediction accuracy and biological interpretability, particularly under contrasting photothermal regimes typical of Mediterranean environments. Beyond developmental traits, the approach provides a blueprint for future integration of multi-omic data (e.g., transcriptomic, metabolomic, and spectral information) to capture both static and dynamic determinants of performance (Wang et al., 2019; Wu et al., 2022, 2024; Crossa et al., 2025a).

The clear superiority of MTME-CV2 models underlines the need to shift from traditional one-trait, one-environment evaluation pipelines to integrated, multi-layered prediction systems. In practical terms, breeders can use MTME predictions to (i) identify lines with consistent performance across environments, (ii) predict untested genotypes in new sowing cycles, and (iii) optimize selection indices incorporating correlated traits (e.g., GN, GW, SL). When combined with economic weights or multi-trait selection indices the MTME framework can directly inform decision-making to maximize multi-trait genetic gain per cycle.

Moreover, the high accuracy of different sowing date predictions supports the implementation of sowing-specific genomic selection, where training populations are updated annually with the latest phenotypic and environmental data. Such dynamic updating has been shown to accelerate genetic gain by maintaining model relevance under changing climate and management conditions (Crossa et al., 2017). For durum wheat breeding programs in the Mediterranean Basin, where season-to-season variability is pronounced, this adaptive strategy could significantly enhance selection efficiency and stability.

While the present study focuses on genomic and phenotypic data, future expansions should integrate multi-modal inputs (e.g, environmental covariables, near-infrared spectra, time-series phenotyping data) to further improve model realism and predictive power. Multi-modal deep learning (MM-DL) frameworks (Montesinos-López et al., 2024; Crossa et al., 2025a, 2025b) can capture nonlinear relationships and trait-environment interactions that conventional kernels might overlook. For durum wheat, integrating environmental covariables aligned with phenological stages (e.g., sowing–heading, heading–flowering, flowering–harvest) has proven particularly valuable in previous CIMMYT studies (Basnet et al., 2019; Crossa et al., 2022; Cuevas et al., 2025).

Ultimately, the convergence of MTME statistical models and multi-modal data integration represents the frontier of prediction-based breeding. By combining the interpretability of kernel-based frameworks with the representational power of deep learning, breeders can more precisely identify genotypes with stable, high performance across variable environments—an essential requirement for ensuring food security in the face of climate change.