The plant materials used in the current study comprised 300 Watkins wheat landraces and 16 non-Watkins landraces and elite cultivars obtained from the Germplasm Resource Unit (GRU) at the John Innes Centre in Norwich, United Kingdom (Supplementary Table S1). The 316 lines were grown in two consecutive seasons, namely 2020–2021 and 2021–2022, in four different agricultural research stations in Egypt with varying environmental conditions, including 1) Sakha Agricultural Research Station in the north Delta (30.0642 ° N, 30.5645 ° E), 2) Nubaria Agricultural Research Station in the West Nile Delta to represent newly reclaimed lands under surface irrigation (30.6973 ° N, 30.66713 ° E), 3) Gemmiza Agricultural Research Station in the Middle Delta to represent heavy and fertile soil (30.867 ° N, 31.028 ° E), and 4) Sids Agricultural Research Station in Middle Egypt, a region characterized by high yielding conditions in the Nile valley (29.076 ° N, 31.097 ° E) (Figure 1). Phenotypic data were collected as the number of days to heading, plant height, plant growth habit, and thousand kernel weight (1,000 KW). The phenotypes were obtained by growing each entry in 3-metre-long rows 30 cm apart. Fertilizers were applied by adding phosphorus fertilizer at 35 kgP ha⁻¹ one dose before sowing, and nitrogenous fertilizer at 180 kgN ha⁻¹ in three separate doses during sowing, 30 days post-sowing, and at the tillering stage. The sowing date was 15 November every season. The agronomic characteristics were recorded as days to heading (DH), the days to heading was recorded as heading when the spike emerged for 0.25 of its length in 50% of the plants (Zadoks DGS 53) and, plant height in cm (PH), was measured as PH, from the soil surface to spike top, and thousand kernel weight (1,000 KW; g) (Figure 2).

The screening traits were consecutive years to examine the phenotypic variation across four diverse ranges of growing conditions under natural field condition in Egypt in conductive 2 years. The mean was calculated as the average of ten plants in each row. The graphical analysis of the evaluated genotypes, including genotype-by-environment (GGE) biplots for days to heading (DH), plant height (PH), and thousand kernel weight (1,000 KW) across four environments over two consecutive seasons, was conducted using GenStat 19th Edition (VSN International Ltd., Hemel Hempstead, United Kingdom) following the methodology of Yan et al. (2000), Yan et al. (2007). The sowing date was 15 November every season. The agronomic characteristics were recorded as days to heading (DH), the days to heading was recorded as heading when the spike emerged for 0.25 of its length in 50% of the plants (Zadoks DGS 53) and, plant height in cm (PH), was measured as PH, from the soil surface to spike top, and thousand kernel weight (1,000 KW; g). A correlation matrix was calculated to determine the degree of correlation between the location and the phenotypes. The two stability parameters of superiority performance (Lin and Binns, 1988) and mean ranks (Nassar and Huehn, 1987) were used to quantify and rank wheat genotypes for good performance and stability where a genotype had the lowest values of the two parameters is considered the most stable one.

We downloaded the Watkins SNP genotype data from the Cereals Data Base (https://www.cerealsdb.uk.net/cerealgenomics/cgi-bin/display_varieties.pl?example=35K_breeders_array&submitter=Submit+Button). This data had previously been generated using the Affymetrix 35K Axiom® Wheat Genotyping Array and screened for genetic variations across the A, B, and D genomes (Allen et al., 2017). TASSEL software (version 5) (Bradbury et al., 2007) was employed to remove SNPs with a minor allele frequency (MAF) lower than 0.05 and a call rate below 90%.

The genetic structure of bread wheat was evaluated using samples from 33 diverse nations in Europe, Asia, and Africa. The samples included 56 landraces from India, 32 from China, 28 from Spain, and others (Supplementary Table S1). Discriminant analysis of principal components (DAPC) was used to identify the differences among populations and generate summarized features for cluster analysis. It was performed without prior information on individual populations using Adegenet v2.1.3 (Jombart, 2008). The find.clusters function was utilized to confirm the DAPC results and determine the optimal number of clusters ( k ) using three different plots. The first plot was based on the lowest value in the “elbow” of the Bayesian information criterion (BIC) curve, the second on a scatter plot of the discriminant functions (DAPC), and the third on a bar plot of each sample’s posterior probabilities of group assignment. To determine the genetic distance between countries, the dist function in the Stats v4.0.3 library (Pellegrini and Dusanter-Fourt, 1997) was utilized. The hclust (Hierarchical Clustering) function was used to determine cluster distance and sample linkage, and the As.phylo function in the ape v5.4.1 library (Paradis et al., 2004) reformat the dendrogram into tree format. The annotated dendrogram was plotted using ITOL software (Letunic and Bork, 2019).

The population structure analysis was employed as a covariate in the GWAS analysis. It was conducted to assess genetic admixture, which refers to the process or outcome of interbreeding between multiple isolated populations within a species (Yuan et al., 2017). The LEA R package version 3.2.0 was used to examine the genetic structure of the core collection. The Snmf (sparse nonnegative matrix factorization) function was utilized within the package, with the number of ancestral genetic groups ( K ) varied from 1 to 10 in ten repeat runs for each K value. The function provided least-squares estimates of ancestry proportions instead of maximum likelihood estimates (François and Durand, 2010). The cross-entropy function was utilized to determine the optimal number of ancestral populations that best explain the genotypic data. The results were presented in a bar plot showing the Q-matrix constructed by the Q() function.

In the current investigation, a GWAS was conducted using the Python-based tool vcf2gwas (Vogt et al., 2022). This tool implements the genome-wide efficient mixed model association (GEMMA) protocol and performs GWAS directly from a VCF file, while also allowing for multiple post-analysis operations (Zhou and Stephens, 2012). The GWAS analysis was executed using a linear mixed model (LMM) and the significant marker-trait associations were determined using the Wald test. The relevance of SNP-trait associations was determined based on the significance threshold of p -value ≤ 0.001 and False Discovery Rate (FDR) ≤ 0.1 displayed on the Manhattan plot. The results were visualized by generating Manhattan plots using the qqman and ggplot2 R packages. The TASSEL program (Bradbury et al., 2007) was used to calculate linkage disequilibrium (LD).

To study the phenotypic variation in the panel under Egyptian field conditions, the panel was grown for two seasons in four agricultural research stations in Sakha, Gemmiza, Sids, and Nubaria. We evaluated the panel for days to heading, days to maturity, plant height, and thousand kernel weight. We classified the correlations between the traits and locations, considering only values greater than 0.5 as significant due to the large sample size of the dataset. In the first season, a highly significant and positive correlation ( r > 0.6 ) was consistently observed for the number of days to heading across the four locations, indicating a similarity in climate during the heading stage across these regions. This positive and highly significant correlation persisted in the second season, particularly in Nubaria when compared to the other regions (Table 1). Regarding plant height, significant and positive correlations were found between Gemmiza and Sakha in the first season, and between Gemmiza and Nubaria in the second season (Table 2). However, for the 1,000 kernel weight trait, despite its statistical significance, the associations among this trait in the four locations were weak ( r < 0.5 ) (Table 3). This suggests that the trait is influenced by environmental factors that may differ between the four locations, such as temperature variation during the grain filling period and the severity of diseases. In summary, a significant correlation was generally observed between the traits as well as across the four field trial locations, except for thousand kernel weight (Table 3).

The two stability parameters of superiority performance [20] and mean ranks [21] were used to quantify and rank wheat genotypes for good performance and stability where a genotype had the lowest values of the two parameters is considered the most stable one (Lin and Binns, 1988; Nassar and Huehn, 1987). Genotype stability itself is not enough as a selection parameter for the aimed genotype unless its performance is good. Bearing in mind that the shortest and earliest genotypes plus that had the heaviest 1,000 kernel weight were desirable. In addition, the stability phenomenon was diagrammatically plotted using Genotype x Environment (G × E) biplot graph. In the case of heading date and plant height, as mentioned earlier the elite genotype must record the lowest values of these traits (below-grand mean that is located on the left side of the origin point in the GGE biplot graph) while the elite genotype regarding 1,000 kernel weight that surpassed the grand mean and laid out on the right side of the GGE biplot graph. Estimates of superiority performance and mean ranks of 30 selected wheat genotypes for heading date, plant height and 1,000 kernel weight are summarized in (Supplementary Tables S2–S4). Results shown in Supplementary Table S1 - Supplementary Figure S7 indicated 30 wheat genotypes that recorded the earliest heading date (less than the grand mean 132.5 days) as well as they recorded the lowest values of the two stability parameters being superiority performance and mean ranks. It is noted that 27 out of 30 wheat genotypes were stable using both stability parameters plus GxE biplot graph. Because the number of genotypes is more than 300, it is difficult to read the results from the crowded biplot graph.

Regarding plant height, there were 26 out of 30 wheat genotypes were stable using both numerical stability parameters and graphical method of GxE biplot graph (Table 1; Supplementary Figure S8). For 1,000 kernel weight, 25 genotypes were characterized by stability using the three methods of stability (Table 2; Supplementary Figure S9). Nine genotypes were revealed stability toward all tested agronomic characters (PH, DH and 1,000 KW) across location, namely, BecardKachu, CIMCOG_3, CIMCOG_32, CIMCOG_47, CIMCOG_49, CIMCOG_53, Reedling, Super 152 and Waxwing. While six wheat genotypes BAJ, CIMCOG_56, MISR1, Pfau, Weebill and Wyalkatchem were revealed stability toward plant height and days to heading in all different locations.

The Watkins wheat landrace collection is composed of 826 worldwide wheat genotypes collected during the 1930s before the onset of intensive breeding (Wingen et al., 2014). In the present study, we used a core set of 300 Watkins accessions previously selected to maximize the genetic diversity of the collection (Arora et al., 2023) and supplemented this with 16 additional entries (Supplementary Table S1). To elucidate the genetic diversity, structure, and relatedness of these 316 genotypes, we analyzed publicly available genotype data previously obtained with the Affymetrix 35K Axiom Wheat Breeder’s Genotyping Array. We filtered this data to retain 9402 SNP markers with a minor allele frequency (MAF) of less than 0.01 and a genotyping call rate of at least 90%. The number of SNP markers within 1,000 Mbp ranged from 36 (chromosome 4D) to 855 (chromosome 2B) (Figure 3).

The linkage disequilibrium (LD) analysis provided using SNP demonstrated how closely SNPs are have a high LD values (Figure 4). This pattern aligns with the expectation that nearby loci are more likely to be in LD. The green line denotes a significance threshold for LD, with values above this line being considered significant, while the blue line provides the average r 2 value for comparison, which was 35 Kbp in our study. The average R 2 values across different wheat chromosomes reveal significant variation. Chromosomes 1A–7A and 1B–7B show LD values ranging from very low (0.000) to high (up to 1.000). Notably, chromosome 2A has a higher median LD (0.0516) compared to chromosome 4A (0.0104), indicating generally higher LD on 2A. Conversely, chromosomes 1B–7B have lower median LD values (0.0174–0.0243) (Figure 4). Chromosomes 1D to 7D exhibit a broader range of LD values, with chromosome 1D having a maximum LD of 1.000. Chromosome 5D has a relatively high 3rd quartile (0.1718) among D chromosomes. The D genome generally shows higher median and maximum LD values compared to A and B genomes, suggesting differing recombination rates or selection pressures. Some chromosomes, like 7A and 4B, demonstrate notable variability, with 7A having a high 3rd quartile (0.4852) and 4B a much lower one (0.0329) (Figure 4).

The GWAS results showed the presence of 202 significant SNPs distributed across all chromosomes, with 3B having the highest frequency (30 SNPs) and 3D and 7D having the lowest (1 SNP) (Supplementary Table S5). Among the four growth indices, plant height was found to be most associated with 77 SNPs, followed by date of heading with 70 SNPs, thousand kernel weight with 44 SNPs, and date of maturity with 10 SNPs (Figures 9–11). The SNP variants presented themselves in several forms of nucleotide changes. The most common change observed was T (thymine) to C (cytosine), with a total of 44 occurrences and an even distribution across all traits. The second most frequent change was C (cytosine) to T (thymine), which was recorded 43 times. The association analysis of thousand kernel weight found significant SNPs on all chromosomes except for chromosomes 3A, 1B, 4B, 1D, 4D, and 5D. The highest number of SNPs affecting this trait was observed on 4A (7 SNPs), followed by 5A (5 SNPs). The trait was also associated with SNPs on chromosomes 2A, 2B, 3B, and 5B (4 SNPs each), chromosomes 6A, 6B, and 7B (3 SNPs each), chromosomes 2A and 6D (2 SNPs each), and chromosomes 1A, 2D, 3D, and 7D (1 SNP each) (Figure 9). All chromosomes except chromosomes 5B, 3D, 4D, 5D, and 7D had SNPs that associated with plant height. The highest number of SNPs associated with this phenotype were found on chromosome 3B, with 20 SNPs, followed by 12 SNPs on chromosome 2B, 8 SNPs on chromosomes 5A and 1B, 6 SNPs on chromosome 1D, 4 SNPs on chromosomes 6A and 7A, 3 SNPs on chromosomes 1A and 6B, 2 SNPs on chromosomes 3A and 7B, and 1 SNP on chromosomes 2A, 4A, 4B, 5B and 2D (Figure 10). All chromosomes except chromosomes 2D, 3D, 4D, 6D, and 7D were found to influence days to heading. Chromosomes 4A and 1D had the highest frequency of SNPs (11 and 10, respectively), while chromosomes 6B and 3D had the lowest (1 SNP each). The remaining SNPs associated with this trait were evenly distributed across the genome (Figure 3). These genotypic variations resulted in phenotypic changes that also influenced plant maturity. One SNP marker associated with this trait was found on chromosomes 3A, 6B, 1D, and 5D, while chromosomes 5A and 2B contained 3 SNPs (Figure 11). We observed only DH, have SNPs shared across different years and stations. Specifically, SNP AX-94406783 was present in the Gemmiza station in 2020 and again in Nubaria in 2021, while SNPs AX-95070278 and AX-94922585 were identified in both Gemmiza in 2020 and Sakha in 2020, highlighting their stability across diverse growing conditions. Additionally, SNP AX-94779755 was detected at Gemmiza in 2021 and Sakha in 2021, further emphasizing its potential relevance to the DH trait across different environments.

The evaluated Watkins landraces demonstrated remarkable stability for all tested agronomic traits, including PH, DH, and 1000KW, across diverse growing conditions. Nine specific accessions: BecardKachu, CIMCOG_3, CIMCOG_32, CIMCOG_47, CIMCOG_49, CIMCOG_53, Reedling, Super 152, and Waxwing—consistently performed well across multiple locations. Additionally, six wheat genotypes—BAJ, CIMCOG_56, MISR1, Pfau, Weebill, and Wyalkatchem—showed stability in plant height and days to heading across all trial sites. Several crosses were initiated between Egyptian wheat cultivars and stable accessions. CIMCOG_32 and CIMCOG_47 were selected for the 2023–2024 crossing block. CIMCOG_32 was used as a parent in multiple crosses, including with the Egyptian cultivar Sakha 95 and advanced exotic wheat lines. These efforts underscore the potential of CIMCOG_32 and Watkins landraces in breeding programs to enhance wheat stability and performance across different environments.

Population structure is a statistical technique used to clarify the genetic composition of individuals within a population as well as the ancestry ratio, demonstrating genetic variance among populations (Patterson et al., 2006). It is mainly used for identifying Subpopulations clusters within a larger population, understanding genetic relationships, and depends on PCA analysis to reduce the dimensionality of genetic data and visualizes genetic variation among samples (Patterson et al., 2006). We assessed genetic structure in a collection of 316 global bread wheat genotypes (Triticum aestivum) representing mostly landraces. The panel was genotyped with the 35K Axiom(B) Wheat Breeder’s Genotyping Array (Affymetrix product ID 550524). The array contains 35,143 SNPs selected to be informative across a diverse global collection of elite and landrace varieties of hexaploid and tetraploid wheats (Winfield et al., 2018). We filtered the markers based on minor allele frequency and missing data and found that our shortlist of 9,402 SNP markers had a genome distribution consistent with earlier research including a two to five times lower marker density in the D genome (Alipour et al., 2017; Allen et al., 2013). The number of anchored markers is largest in the B genome, followed by the A genome, and lowest in the D genome.

There is significant variation in average R 2 values across different wheat chromosomes. Chromosomes 1A–7A and 1B–7B show LD values ranging from very low to high, with 1B–7B generally having lower median LD values. Chromosomes 1D to 7D exhibit a broader range of LD values, with higher values observed in the D genome compared to the A and B genomes, indicating differing recombination rates or selection pressures. Notable variability is seen in chromosomes like 7A and 4B. This observed trend is consistent with previous wheat studies. For instance, Mourad et al. (2020) reported that the D genome had the highest significant LD, with an average R 2 of 0.887, followed by the A and B genomes. Roncallo et al. (2021), found that 13.4% of marker pairs exhibited significant LD, with a high LD ( r 2 > 0.7 ) in only 0.94% of comparisons, and an LD decay of 11.8 Mb. Additionally, Aleksandrov et al. (2021), noted that in an old germplasm collection, significant LD ranged from 6% to 27% across different chromosomes, while in a modern collection, it ranged from 10% to 45%. These findings align with our results, which indicate considerable variation in LD across different wheat chromosomes and higher LD values in the D genome compared to the A and B genomes.

A cluster is a group of objects that are closer and more similar to each other than those outside of the group (Grabowski et al., 2018). The DAPC analysis revealed three distinct clusters within the populations, with a total diversity score of 11.6 across both axes. This suggests that while the DAPC could identify relatively low diversity ratios between samples, there is substantial genetic diversity among them. This value helps us understand the extent of variation captured by the analysis. This ratio is relatively low compared to the 60.1% of total variance reported by Fiore et al. (2019) for their study involving twenty-seven durum wheat varieties and one bread wheat Sicilian landrace. This suggests that our observed diversity, as indicated by the sum of 11.6, is lower relative to their findings, potentially indicating less pronounced clustering or different levels of genetic variation in our sample. This value is relatively high compared to the second and third principal components reported by Alemu et al. (2020) in their study of Ethiopian durum wheat (Triticum turgidum ssp. durum). Their study highlighted different patterns of variance, suggesting that our results reflect a higher level of diversity or variability in the principal components analyzed. Afghanistan, Iran, India, China, Iraq, and Burma comprise the first cluster (the majority of the Asian cluster). The second cluster includes most European countries, including Italy, Spain, Portugal, Turkey, and Greece. The third cluster includes African countries such as Egypt, Algeria, and Ethiopia, which are not in the same cluster but are close. The Bayesian information criterion (BIC) is a well-known and commonly used method in statistical model selection. It was applied to approximate a transformation of the Bayesian posterior probability of a candidate model (Neath and Cavanaugh, 2012). To identify the optimum cluster (K)(knee), we displayed BIC in three forms, depending on the DAPC; the first plot was a boxplot, which showed that three and four were significant values for K. The second was a scatter plot that classified the samples into three categories. The third was a bar plot that performed with K = 2 and K = 3 and identified those three as the best group. The clustering findings are consistent with the population structure, which revealed that while K = 4 is the best, K = 2 and K = 3 are also significant. It should be highlighted that the current results of population structure analysis are consistent with our previous findings of the entire set (804 accessions) of hexaploid wheat (Winfield et al., 2018). This consistency reaffirms the reliability of our present findings, indicating a robust and persistent population structure across diverse datasets. Our findings on the clustering of the Watkins landrace populations could provide complementary insights to Pont et al. (2019) identification of selection footprints and evolutionary history in modern wheat. Both studies emphasize the complex genetic landscape of wheat and contribute to a broader understanding of its genetic evolution and diversity.

We conducted a comprehensive GWAS to explore the influence of single nucleotide polymorphism (SNP) markers located on genes controlling four important traits in wheat: days to heading, days of maturity, plant height, and 1,000 kernel weight. The identified SNPs have potential as marker-trait associations (MTAs) to guide breeding programs in the agricultural sector. The analysis of 203 SNP markers yielded 31 genomic regions associated with the studied traits. Notably, three SNPs on chromosome 5A were found to be linked to 1,000 kernel weight, with a physical distance of less than 10 Mb. Additionally, four clusters evenly distributed on chromosomes 3B and 1D were identified as influencing plant height. Thirteen SNPs in clusters on chromosomes 2A, 5B, 6B, and 1D were associated with days to heading. For the trait of maturity date, three SNPs were located on chromosomes 2A, 2B, and 6B. Interestingly, six SNPs were found to be linked to both 1,000 kernel weight and days to heading. These results provide valuable insights for breeders seeking to improve these important traits in wheat and lay the foundation for further functional studies on the identified markers.

To map the detected probes onto the wheat genome and identify candidate genes and encoding protein domains influencing the traits, the GrainGenes database (https://graingenes.org/GG3/) was utilized. Notable associations were observed between specific SNPs and genes controlling the traits. For example, SNP AX-94462177 with a p-value of 0.0002 and 0.4 FDR, associated with 1,000 kernel weight, was detected within the Kinesin-like protein domains (NPK1/TraesCS5A02G317000), known for its role in seed development in rice; these domains play crucial roles in almost all biological processes in plants (Li et al., 2011). SNPs AX-94643695 and AX-94814458 with p-values of 0.00004 and 0.00008, and FDRs of and 0.3, associated with wheat date of heading, were in the ATP-sulfurylase PUA-like gene (TraesCS2D02G031800), indirectly influencing plant development (Xiao et al., 2022; Khan et al., 2007). SNPs AX-94446101 and AX-95231601 with p-values of 0.0004 and 0.0005, and FDRs of 0.904116745, also associated with date of heading, mapped to the Serine/threonine-protein kinase D6PK-like gene (TraesCS5B02G252500) which serves as a lipid domain-dependent regulator of root epidermal planar polarity in Arabidopsis (Stanislas et al., 2015). A cluster of loci (AX-94510523, AX-94425015, and AX-94632604), with p-values of (0.0003, 0.00008, and 0.00057), and FDRs of (0.7,0.3, and 0.6) associated with date of heading was located in the wheat C2H2 ZINC finger transcription factor gene TraesCS4B02G003500 encoding a C2H2-type zinc finger protein, reported as the best candidate gene in the QTL Qhd.2AS controlling wheat growth and development (Li et al., 2023). SNP AX-94983266 with a p-value of 0.0006 and FDR of 0.5, associated with date of maturity, located to the HAUS1 (HAUS Augmin Like Complex Subunit 1/TraesCS3A02G380800) gene. In Arabidopsis, the AUGMIN complex impacts spindle and phragmoplast microtubule arrays during sexual reproduction (Oh et al., 2016). Lastly, SNP AX-94409249 (p-value of 0.0007, and FDR of 0.6), associated with wheat yield and days to maturity, mapped to the RING finger domain gene TraesCS3B02G139600 on chromosome 3B (Supplementary Table S6). In Arabidopsis and tobacco, RING zinc finger genes are involved in seed development and stress resistance (Xu and Li, 2003; Zeba et al., 2009; Han et al., 2022).