To conduct our comparative study, we used a large dataset widely recognized in the field of crop genetics, published by Juliana et al. (2019). This dataset, provided by CIMMYT (International Maize and Wheat Improvement Center, Mexico), includes a total of 44,400 wheat lines for which 78,606 SNPs were identified through genotyping-by-sequencing (GBS) and 50 phenotypic traits were evaluated. After thoroughly examining the CIMMYT dataset, we created two subsets for our analysis (Figure 1A): the first subset was selected for the trait plant height (PH) including 7,887 individuals evaluated in first-year yield trials conducted in Ciudad Obregón, Mexico, during the 2014–2015 season, with two replicates, resembling Juliana et al.’s panel 4. The second subset targeted the trait thousand kernel weight (TKW) and comprised 3,478 individuals derived from multiple trials evaluated in Ciudad Obregón, Mexico: 766 lines from the elite yield trial (EYT) 2013-2014, 775 lines from EYT 2014-2015, 964 lines from EYT 2015-2016, and 980 lines from EYT 2016-2017, resembling Juliana et al.’s panel 1 (but seven individuals less due to removal of individuals with missing values in the quality traits). For more information on the yield trials be referred to Guzman et al. (2016) and Juliana et al. (2019). These datasets were subsequently filtered using TASSEL.v5 (Bradbury et al., 2007) as described in Juliana et al. (2019): first, a filtering of unaligned markers (UN) was done resulting in a final set of 78,148 markers (cf. Figure 1 of Juliana et al.), followed by a second filtering step for markers with less than 40% missing data, a MAF > 0.05, and less than 5% heterozygosity. Missing markers were imputed using the LinkImpute (Money et al., 2015) module within TASSEL.v5, with the following settings: Linkage Disequilibrium (LD) was set to 60, and k-Nearest Neighbors (kNN) to 20. Using the tool´s internal masking procedure to assess accuracy, the imputation error was estimated at 8% (PH) and 9% (TKW), respectively. After the imputation step, SNPs were again filtered for minor allele frequency (MAF > 0.05), resulting in the final datasets as follows: the PH dataset included 7,890 markers and 7,886 individuals, and the TKW dataset comprised 7,287 markers and 3,478 individuals. In our study, the resulting final number of markers differed from those reported by Juliana et al. (2019), who identified for panel 1 6,355 markers. The number of markers of panel 4 is not explicitly stated but falls in the range between 9,171 and 9,704 markers. In addition to this ambiguity regarding panel 4, especially the lack of information on thresholds, settings, and imputation accuracy made it impossible to replicate the filtering procedures. Furthermore, additional factors contributing to this discrepancy may be attributed to differences in TASSEL software versions. A note to the file formats: the initial dataset was in HapMap format, however, in order to use the data for Gemma and GCTA, it was converted using PLINK 2.0 (Chang et al., 2015) into bed/bim/fam format.

Code including the tools sommer and GAPIT, ML algorithms, TSLRF and plotting was run in R Studio server version 2025.5.0.496 with R version 4.4.3 (R Core Team, 2025) on a Ubuntu 22.04.5 LTS server equipped with an Intel(R) Xeon(R) CPU E5–2637 v3 @ 3.50GHz (16 cores, 32 threads). Manhattan & QQ-plots were plotted with ggplot2 version 3.5.2. UpSet plot and Venn diagram were visualized using the R-package ggVennDiagramm version 1.5.2. GCTA, GEMMA and TASSEL were used over the command line. The code (separated in seven R scripts: EN.R [commit tag: 3fcdd80], FeatureSelection [3fcdd80], Gapit.R [3fcdd80], RF.R [2385148], Sommer.R [3fcdd80], TSLRF.R [09f8d51], XBR.R [c20e9f2]) is publicly available on https://github.com/MolinLab/BenchmarkingGWAS. The overall computational workflow is depicted in Figure 1B and is given in more detail below:

Linkage disequilibrium (LD) calculated with PLINK 2.0 (Chang et al., 2015) was used to define flanking regions around significant MTAs. Candidate QTLs and genes were identified based on the Triticum aestivum reference genome (IWGSC; Ensembl Plants, accessed 18 May 2025; Dyer et al., 2025). Within each defined window, genes were annotated based on gene ontology, functional classification, and prior literature to prioritize and report the most relevant candidate genes.

To explore ML for GWAS, we employed the R package mlr3, which provides a robust base set of learners that can be extended with additional packages, e.g. EN, RF, and XGB. SVM was implemented with package e1071, but we encountered a complex prediction error and therefore did not investigate this approach further within this study. Additionally, we evaluated TSLRF, which integrates population stratification using the FASTmrEMMA method (Wen et al., 2018) and feature selection via LARS, followed by a final prediction with a RF model. We acknowledge that our study does not cover all available methods and tools, but the chosen approaches provide a strong foundation for assessing how ML can enhance GWAS in wheat breeding. The selected models provide interpretable feature-importance metrics, which are essential for identifying MTAs. In addition, they represent complementary algorithmic approaches, enabling a diverse yet comparable set of importance metrics that can be directly contrasted with traditional GWAS outputs.

Using McNemar’s test (Supplementary Figures 14, 15), we assessed whether the different models identified significantly different sets of MTAs. For TKW, the results revealed a clear methodological separation among the approaches. The MLM-based approaches, together with TSLRF and EN, showed no significant pairwise differences, indicating that these methods detected largely overlapping MTAs. In contrast, FarmCPU, BLINK, RF, and XGB produced MTA sets that differed significantly from the MLM cluster. Within this non-MLM group, BLINK and RF as well as RF and FarmCPU did not differ significantly, whereas most other comparisons were significant. Notably, XGB and EN did not differ significantly, suggesting partially overlapping signals despite their methodological differences. Importantly, XGB differed from most MLM approaches only at α = 0.05, whereas FarmCPU, BLINK, and RF showed much stronger divergence at α = 0.001.

For PH, a similar pattern emerged. The MLM methods that detected MTAs showed no significant pairwise differences, while all other models yielded MTA sets significantly different from the MLM cluster, with the exception of TASSEL, which did not differ from XGB. Again, XGB and EN showed no significant difference, although EN differed significantly from all MLM approaches for PH, in contrast to its alignment with MLM methods for TKW. FarmCPU and BLINK did not show similarity to any other model, with the only non-significant comparison observed between RF and FarmCPU. Similar to the TKW results, XGB differed from MLM approaches only at α = 0.05 for PH, except for the comparison with TASSEL, which was not significant. Influence of population structure on MTAs.

In terms of computational efficiency on the Intel(R) Xeon(R) CPU E5–2637 v3 @ 3.50GHz (16 cores, 32 threads), the tools demonstrated considerable variation. It is important to note that the server was used in parallel by multiple users during these experiments, which may have introduced variability in the reported computation times due to resource contention. Under these conditions, TASSEL exhibited the longest total computation times, ~1,025 minutes (min) for TKW and 26,220 min (~ 18 days) for the PH. GAPIT showed large variation across models: for TKW, runtimes ranged from 1.5 (FarmCPU) to 1,851 min (CMLM), and for PH, from 4 min (FarmCPU) to 19,551 min (CMLM). Notably, FarmCPU and BLINK consistently ran significantly faster than other GAPIT models. Sommer completed in 13 min (TKW) and 311 min (PH). Command-line tools such as GCTA required ~2 min (TKW) and 10 min (PH), while GEMMA took 3 min (TKW) and 36 min (PH). For the ML models implemented in this study, computational times included both hyperparameter tuning and model training, and were strongly influenced by the chosen hyperparameters. In general, the models ranked from fastest to slowest as follows: EN (TKW: ~ 2 min; PH: ~ 12 min) < RF < XGB (TKW: 100 min; PH: 344 min). EN was the fastest, which may be attributed to its efficient regularized regression framework. RF and XGB, on the other hand, incurred higher computational costs, partly due to hyperparameters that were computationally intensive. For instance, we used 1,000 rounds of boosting for XGB, and reducing this number would decrease computational time but would result in less performance.

Even though most of the traditional tools primarily rely on a MLM, subtle differences in implementation and internal computations resulted in variations in significance levels. In most cases, these differences were minor and led to a few MTAs falling just above or below the significance threshold depending on the tool used. However, the biggest distinction across tools was computational performance and model flexibility. TASSEL provided an easy-to-use platform with several in depth filtering options, making it a solid choice for pre-processing, SNP filtering, and imputation, and has been used already several times in wheat GWAS (e.g. Sukumaran et al., 2015; Liu et al., 2017; Genc et al., 2019).

However, in our study, TASSEL exhibited notably long computation times, including a runtime of up to 2.5 weeks for the PH dataset. While this may partly reflect less efficient computation and parallelization constraints, it is important to consider that the server was shared among multiple users, which could have influenced the observed runtimes. TASSEL does offer additional modelling options, such as GLM and Fast Association, which can significantly increase computation time but may not always adequately account for kinship and population structure. Sommer, implemented in R, provides users with the flexibility to define complex MLM structures and random effects tailored to experimental designs. While its runtime was moderate in our tests, its customizability and compatibility for HapMap-formatted genotype data make it a valuable tool. GCTA and GEMMA, both command-line tools optimized for PLINK input, consistently demonstrated faster performance compared to TASSEL and R-based tools. However, converting HapMap files into the required input format may present a challenge for some users. Among the tools evaluated, GAPIT offered the widest range of model options, allowing users to select the best model based on genomic control or benchmark results. Computational times varied significantly across GAPIT models, with FarmCPU and BLINK running considerably faster than others. While CMLM was introduced to reduce computational time (Zhang et al., 2010), the process of calculating the optimal compression level, such as the 158 grouping variations for PH, substantially increased runtime compared to other GAPIT models. This observation might be data size-dependent and could be more pronounced with datasets containing several thousand individuals.

GAPIT’s variety also highlights how different model assumptions influence results. For instance, the basic GLM tests one marker at a time with only population structure covariates without kinship matrix, to account for relationships among individuals. As a result, GLM often produces many spurious associations when individuals are related, or population structure is complex (Yu et al., 2006). This behavior is reflected in our study and other experiments and simulations (Kaler et al., 2020; Sandhu et al., 2024). In SUPER, kinship is built from a subset of informative markers, and markers in LD with each test SNP are excluded to avoid confounding (Liu et al., 2016). However, using only a few pseudo-QTNs often fails to capture the full relatedness structure, reintroducing confounding. Thus, SUPER can still suffer high inflation if the pseudo-QTNs do not fully account for cryptic kinship, especially in structured wheat populations. The high amount of spurious associations we observed is in line with simulation studies (Wen et al., 2018; Kaler et al., 2020; Sandhu et al., 2024). MLM approaches i.e. MLM, CMLM, MLMM performed well in controlling for population stratification and relatedness which is reflected in the genomic inflation values. However, this strict control comes at the cost of power to detect true associations (Kaler et al., 2020; Huang et al., 2019). By modelling the polygenic background, MLM-based tests effectively attribute variance to kinship and leave little for each SNP. As a result, only very strong MTAs will pass strict significance thresholds (Liu et al., 2016). Simulations under high heritability showed that MLM/CMLM detected five out of 20 true QTLs, while MLMM was able to find 11 out of 20 (Sandhu et al., 2024). This result is mirrored in the results for the trait TKW in this study. No association was detected for PH using GAPIT’s MLM approach, which could be explained by overfitting of kinship. Although when using different clustering in CMLM, the model detected one MTA, showing a slight increase in power. FarmCPU and BLINK showed moderately inflated λ, which might reflect spurious signals. Empirically, FarmCPU often finds more associations than MLM but allows for more false positives, while BLINK efficiently detects associations often with less false positives, outperforming FarmCPU (Huang et al., 2019). When assessing potential confounding due to population structure, we observed that certain MTAs uniquely detected by FarmCPU and BLINK indeed showed higher informativeness values, indicating a greater likelihood of being influenced by population structure. This may reflect differences in the modeling approaches of these methods and their interaction with the specific characteristics of the analyzed dataset. Nevertheless, not all MTAs detected by BLINK and FarmCPU appear to be related due to population structure. The higher number of significant markers observed from these methods likely arises due to application of multi-locus tests in FarmCPU and BLINK, enabling detection of loci that single-marker tests might miss (Merrick et al., 2022). However, the potential of false positives, as indicated by the model inflation underscores the value of incorporating ML approaches as a benchmark to further validate and refine these findings.

The genetic architecture of the traits TKW and PH in wheat is clearly more complex than can be captured by purely linear, single marker models. In order to account for the well‐known limitations of traditional GWAS, namely their inability to capture non-linear and epistatic interactions, we employed ML-mediated GWAS to complement conventional methods. Yoosefzadeh-Najafabadi et al. (2022) demonstrated that ML-based GWAS can enhance the detection of MTAs in soybean, providing a broader and more comprehensive genetic insight than conventional methods alone.

Unlike traditional models, ML approaches do not estimate direct allele substitution effects. Instead, we focused on their capacity to identify MTAs, QTLs, and candidate genes. To further assess their predictive performance, we evaluated overall model fit using RMSE, where XGB outperformed other models for both traits, PH and TKW. TSLRF was developed with the same reasoning as MLMs, that control for population structure and polygenic background, which is crucial to reduce false positives, i.e. importance rankings of SNPs (Price et al., 2006; Sun et al., 2019). Previous work by Stephan et al. (2015) showed that RF and Lasso correcting for population structure outperformed their unadjusted counterparts. However, this stringent control might come at the expense of power, which is reflected in high false negative rates of MLM (Liu et al., 2016). This trade-off was evident in our results: for TKW, TSLRF only detected markers on chromosomes 2A and 6D, consistent with MLM findings. Notably, eight of the ten most important markers overlapped with MTAs from other models, suggesting that a less stringent cutoff might reveal additional meaningful associations, similar to traditional approaches. For PH no marker was significantly associated which is consistent with some MLM approaches used in this study and further supports the notion that stringent control may decrease detection power in traits with more complex architectures.

In contrast, RF without a control for population structure and relatedness, identified a larger number of MTAs. This notable difference to TSLRF indicates that some of the MTAs detected may partly reflect confounding due to population structure. However, when investigating marker informativeness along the PCs this does not seem to be the dominant cause. Another factor could be that, unlike EN and XGB, RF lacks feature regularization resulting in correlated markers receiving similar importance scores. This is more likely, especially because most MTAs identified in RF cluster at Chromosome 6A for TKW and 5A and 7A in PH. Investigation of these clusters revealed that most of the markers are in strong LD, thereby inflating the number of associations. This discrepancy to EN and XGB suggests that regularization might be very important for ML-based GWAS approaches. Despite these limitations, the unadjusted RF model was able to detect MTAs with showed strong biological relevance, with a pipeline that better accounts for correlated features, RF might be a powerful algorithm for GWAS. Previous studies showed the high potential of RF in detecting MTAs: both, RF and gradient boosting machines identified regions with multiple candidate genes and some previously not reported in cattle (Alves et al., 2022; Li et al., 2018). In soy, RF detected more regions compared to traditional models (MLM and FarmCPU) for yield component traits, with no overlap between the approaches (Yoosefzadeh-Najafabadi et al., 2022). In contrast, our study revealed partial overlap between RF and traditional approaches, suggesting that while RF uncovers additional loci, it can still capture some of the core signals detected by traditional GWAS methods.

XGB, designed for high-dimensional and complex data, is well-suited for GWAS applications. Gill et al. (2022) showed that XGB not only outperformed deep learning architectures on most traits but also detected loci associated with flower and seed coat color that are supported by previous research. For complex traits with non-additive gene action gradient boosting performed better than traditional RF and deep learning architectures (Abdollahi-Arpanahi et al., 2020). In our study, compared to mixed‐model single‐marker tests, XGB identified more MTAs, but found fewer than multi‐locus methods like FarmCPU and BLINK. Due to efficient regularization XGB was able to better able to rank markers compared to RF. For TKW, 13 MTAs were detected, seven of which overlapped with those detected by RF. Notably, three MTAs were exclusively identified by these ML algorithms, two of which were linked to biologically plausible gene candidates, highlighting the potential of ML to uncover possible non‐linear or epistatic interactions overlooked in single‐marker mixed models.

The regularized regression model GLMnet is commonly used for feature ranking (Dagasso et al., 2021) or feature selection (Cho et al., 2010; Waldmann et al., 2013). John et al. (2022) compared several ML algorithms on prediction and feature selection performance and showed that EN performed as one of the best selection algorithms with good prediction performance outperforming other ML algorithms. However, EN reliance on regularization can lead to unstable features and may fail to select all informative features. To address this, Meinshausen and Bühlmann (2010) proposed stability selection, a procedure that aggregates selection results across subsamples to improve consistency and control false discoveries. Alexander and Lange (2011) investigated this process on the Wellcome Trust Case-Control Consortium human dataset, however, in this case Lasso was used instead of EN, finding that stability selection effectively controls the family-wise error rate but suffers from a loss of power. Nevertheless, in our study, we implemented stability selection alongside EN because the ridge component of EN preserves groups of correlated predictors, resulting in a less sparse model, while stability selection further improves the robustness and interpretability of our marker rankings. Although EN does not model non‐linear or interaction terms directly, it uncovered several distinct MTAs linked to genes associated with the traits. For PH, a notable example was S7A_22981667, which neither tree‐based methods nor MLMs detected. These unique SNPs likely represent additive effects that were too subtle to surpass genome-wide thresholds in single-marker MLMs yet survived the regularization and stability‐selection process.

Overall, ML-based approaches recovered core MTAs identified by traditional GWAS while differing in their sensitivity to additional loci, with overlap patterns strongly dependent on the specific algorithm and trait. Interpretations involving MLM-based comparisons should nevertheless be made cautiously, as the small number of detected MTAs resulted in sparse contingency tables and limited power, despite the use of a continuity-corrected McNemar’s test.

To assess the added value of using ML approaches for GWAS and candidate gene discovery, we examined the MTAs previously identified by ML approaches for their biological function relevant for the analyzed traits. After removing MTAs in strong LD, we found 21 markers with relevant genes, 11 for TKW and ten for PH. First, we will discuss those MTAs that were found in traditional methods and further validated using ML-based GWAS approaches.

The TKW-associated marker S2A_143724068 was identified by all models except RF and XGB and points to the candidate gene TraesCS2A02G183900. This gene was previously reported by Liu et al. (2022) as an auxin-responsive differentially expressed gene involved in phytohormone metabolism and signaling in pairwise near-isogenic lines. Auxin plays a critical role in lateral root initiation and development. Root spatial configuration in the soil is recognized as a promising strategy to enhance crop yield (Lynch, 2022). Following this reasoning it might be a promising candidate for TKW. Marker S3D_181564815 (FarmCPU, BLINK, EN) co-localizes with three candidate genes potentially influencing TKW. The first one, TraesCS3D02G191300 encodes an expansin, a cell-wall–loosening protein promoting cell expansion in growing tissues such as the developing endosperm. At the molecular level, expansin expression is a key determinant of grain cell enlargement and thus final kernel size. Although 241 expansins have been annotated in wheat, only a subset has been functionally characterized, with some linked to grain traits (Mira et al., 2023). TraesCS3D02G189900, a Class III peroxidase, modulates cell-wall stiffening and has been associated with spikelet traits inversely related to TKW. Overexpression of the related gene TaPRX-2A reduced grain number per spike and spike length, indicating potential pleiotropic effects (Zang et al., 2025). Additionally, TraesCS3D02G189500, a putative nitrate transporter (NPF family), is implicated in nitrogen use efficiency which is a key determinant of grain filling and starch accumulation (Kumar et al., 2022). The gene TraesCS5A02G359900, annotated as a serine/threonine protein phosphatase, was the most likely candidate for marker S5A_561661320 (Blink, RF and XGB). Hu et al. (2012) identified a serine/threonine protein phosphatase gene as a candidate after fine mapping of a grain length, grain width and TKW QTL in rice. A promising candidate for marker S6A_481736675 (FarmCPU, RF, XGB) is TraesCS6A02G259000, a MADS-box transcription factor serving as a candidate gene for grain size and weight in Triticum turgidum (Mangini et al., 2025). Marker S6A_484065257 (FarmCPU, RF, XGB) lies near TraesCS6A02G260600 and TraesCS6A02G260500 that share a seed−storage helical domain and are annotated for lipid transport (GO:0006869). These non-specific lipid transfer proteins have been linked to cuticle formation and seed coat integrity and knock-out plants in rice showed decreased TKW and grain size (Li et al., 2023b). Located near marker S6A_531394936 (BLINK, RF, XGB) is GS1 (TraesCS6A02G298100), a central enzyme in the glutamine synthetase (GS)/glutamate synthase (GOGAT) cycle which is involved in N assimilation and represents a crucial metabolic step in regulating nitrogen use efficiency and grain yield (Fortunato et al., 2023). Additional studies on the GS1 gene family in several crops have shown that they play different roles in determining seed size and number, grain filling and development (Tabuchi et al., 2005; Bernard et al., 2008; Fortunato et al., 2023). In the LD interval also lie TraesCS6A02G296500 (Manjunath et al., 2023) and TraesCS6A02G296400 (Devate et al., 2023), which have been reported as QTLs for TKW. Marker S6D_82630628 (GCTA, GEMMA, sommer, RF, XGB, and EN) falls directly on TraesCS6D02G116200 (TaPRX-2A), encoding a Class III Peroxidase gene previously shown to regulate grain number per spike in common wheat (Zang et al., 2025). Because grain number per spike and TKW are often negatively correlated, variation in TaPRX-2A may indirectly influence TKW. Similarly, TraesCS7A02G136600, a COBRA-like protein showing a range of functions and participate in various developmental processes in cereals. In rice, for example, it has been associated with QTLs for grain number per spike (Mizuno et al., 2021). While this gene also showed no direct correlation with TKW, the study highlighted the complex genetic architecture in spike traits. S7A_90336792 linked with TraesCS7A02G136600 was significantly associated with FarmCPU, RF and XGB.

Regarding the trait PH, the marker S5A_563358779 (FarmCPU, BLINK, XGB, EN), which corresponds to a cluster of Cytochrome P450 (CYP450) genes (TraesCS5A02G362400, TraesCS5A02G361800, TraesCS5A02G361700, TraesCS5A02G363000, TraesCS5A02G362900, TraesCS5A02G362600, TraesCS5A02G362700), has been detected in our study. Some members of the CYP450 family play a key role in brassinosteroid biosynthesis, which influences plant height through the regulation of cell elongation. The biosynthesis and catabolism of brassinosteroids rely on the enzymatic activity of several CYP450 subfamilies (Li and Wei, 2020). Manipulating the expression of BR-related genes can result in pleiotropic effects, some of which may be undesirable in crop breeding. For example, changes in BR pathways often lead to altered plant architecture, including reduced or excessive plant height (Yamamuro et al., 2000). Furthermore, the markers S1B_568961025 (FarmCPU, BLINK, RF) and S7A_403940100 (BLINK, FarmCPU, MLMM, sommer, TASSEL_MLM, RF) lie close to a Growth-Regulating Factor (GRF) and zinc finger GRF, respectively. GRFs are plant-specific transcription factors known to play critical roles in plant growth, development, and responses to environmental stress (Wu et al., 2014; Van Der Knaap et al., 2000; Bazhenov et al., 2021), demonstrating the essential role of GRFs as candidate genes for PH.

The high amount of overlapping markers shows ML-driven GWAS can effectively validate and prioritize potential MTAs that are detected in traditional methods like FarmCPU and BLINK, which show a higher power to detect associations at a cost of more false positives. However, in addition to MTAs that co‐localize with genes of known function, we also identified several associations without any obvious candidate annotation. Two of these unannotated loci stood out for their consistency across methods: S6D_241296319 associated with TKW was identified by all approaches and S7A_335064522, associated with PH detected by six methods. These findings reflect the fact that GWAS is only as informative as the available functional annotation, and without high‐confidence gene models or known pathways, the biological interpretation remains uncertain.

Notably, the marker S2B_565068485, which was detected by FarmCPU, BLINK, and XGB, showed no obvious candidate genes in its vicinity but instead exhibited a signal pattern more closely aligned with population structure, suggesting potential confounding. A similar pattern was observed for S4A_679187682, identified by BLINK and both tree-based learners. These findings indicate that ML-based GWAS can be susceptible to residual confounding when stratification is not accounted for. However, some markers, for instance, S6D_143959985 or S6D_241296319 identified by most approaches applied in this study including the conservative MLM showed similar patterns. A likely conclusion might be that at least some of these signals may reflect genuine selection-driven associations rather than structural or noise-related artifacts.

While the overlap between ML-based and traditional GWAS highlights the robustness of the detected MTAs, we now focus on the ability of ML algorithms to detect additional, previously overlooked associations to extend and enrich conventional discovery pipelines, building on prior studies that have already demonstrated this potential (Mieth et al., 2021; Alves et al., 2022; Gill et al., 2022; Ghose et al., 2024).

Tree-based ML algorithms uncovered three novel MTAs that were not detected by traditional GWAS. For TKW, a candidate behind S2A_507030190 (RF, XGB) is TraesCS2A02G294600, blue‐light inhibitor of cryptochromes 1 (BIC1), functioning as a transcriptional coactivator that promotes brassinosteroid signaling and plant growth (Yang et al. 2021b). Since brassinosteroids might directly influence grain size through cell‐expansion pathways, this MTA represents a physiologically plausible contributor to TKW variation. For PH, S5A_519892082, found by RF and XGB, and S4D_455321401, solely found by XGB, showed several promising genes in their proximity: the first marker is near TraesCS5A02G307300, a PHD-type domain-containing protein, playing an important role in plant growth, development and response to abiotic stress (Pang et al., 2022) and two CYP450 genes (TraesCS5A02G307700, TraesCS5A02G307600) linked to PH. The second marker is associated with TraesCS4D02G283000, a MYB transcription factor, playing important roles in plant development, including internode elongation and stem architecture. Although relatively few studies have directly linked MYB transcription factors to stem development, available evidence suggests their involvement in regulating stem height and strength. For instance, overexpression of TaMYB18s led to a pronounced leaf rolling phenotype, a phenotype also related to drought stress (Kadioglu and Terzi, 2007), and a significant reduction in plant height in wheat (Zhang et al., 2016). Similarly, GmGAMYB has been shown to promote stem elongation via the gibberellin signaling pathway in soybean (Yang et al., 2021a), while OsMPH1 acts as a positive regulator of stem height in rice (Zhang et al., 2017). These findings support the potential involvement of MYB-related genes in regulating plant height across different species. On the other hand, by applying EN, most MTAs for TKW overlapped with markers that were found with traditional methods, while also finding two unique markers. One of them, S1A_540136972, is located near the gene TraesCS1A02G357800, which is involved in the regulation of flower development (GO:0009909). This gene belongs to the AP2/ERF transcription factor family, previously shown to increase grain yield due to regulation of inflorescence architecture and grain number per spike in wheat and soybean (Wang et al., 2022; Li et al., 2023a). The increase of grain number per spike is often linked with a decrease in grain weight directly influencing TKW.

For PH, EN uniquely identified 14 MTAs; out of this, the following five were associated to the candidate genes S1B_670929735, S3A_646672662, S4A_4593612, S7A_22981667, and S7B_695381826. The first one is located near a cluster of S-type anion channel genes (SLAH2 family: TraesCS1B02G456200, TraesCS1B02G456300, TraesCS1B02G456400, and TraesCS1B02G456500), which are associated with nitrate transport (Kumar et al., 2022) and regulation of nitrate flux between root and shoot (Maierhofer et al., 2014), thereby potentially influencing overall plant growth and development, including height. The other marker of relevance is S3A_646672662, which is located near TraesCS3A02G399800, a homologue of GA20ox-D3, which is catalyzing a late step in gibberellin biosynthesis and is responsible for unidimensional cell growth (Qin et al., 2013). This gene underlies classic “Green Revolution” height phenotypes in rice and wheat (sd1/Rht) and controls stem elongation in legumes (Spielmeyer et al., 2002). The other marker, S4A_4593612, encodes ABC transporter B family member 19, which is well-characterized as a major auxin exporter in Arabidopsis. The phytohormone auxin plays a critical role in plant development, in cell growth, division, and expansion (Zhao et al., 2013). This could directly influence the plant height in wheat, making this a promising candidate gene. Furthermore, S7A_22981667 is located close to the gene GRF5-A1 belonging to the GRFs, which have been mentioned already in the chapter above. And finally, S7B_695381826 is located close to TraesCS7B02G426100 which is encoding a F-box domain protein, known to be relevant for various functions (Rameneni et al., 2018).

Collectively, these novel ML-derived MTAs enrich our understanding of the genetic architecture of TKW and PH by revealing hormone-related, nutrient-transport, and protein-turnover genes that traditional GWAS alone failed to detect. Especially, penalized models efficiently rank markers and show promise as ML-based GWAS approaches. Moreover, several of the MTAs mapped to genes involved in biotic and abiotic stress responses, pathways we did not explicitly include, but which can profoundly influence resource allocation, developmental timing, and ultimately affecting grain weight (TKW) and stem elongation (PH) under field conditions.