In this study, 335 spring wheat genotypes (including 325 European and 10 exotic genotypes from CIMMYT, China and Africa) were used to assess FHB resistance under field and controlled conditions. Initially, a set of 300 genotypes (breeding lines and varieties) from different origins, such as Norway, Estonia, Latvia, and Lithuania were genotyped and phenotyped for other agronomic traits, except for FHB resistance (Aleliūnas et al., 2024; Lin et al., 2024). Furthermore, for this study, we included 30 advanced Lithuanian breeding lines, the Slovenian cultivar PS Perlička, and the Chinese resistant landrace Wangshuibai that showed better performance in FHB resistance during one or two field inoculation cycles to the given set. During genetic quality control, three genotypes were excluded. Finally, the set comprised 332 genotypes, of which 65 originated from Nordic countries, 181 from Baltic countries, 76 from Central and Western Europe, and 10 were exotic (foreign origin). There were well-known check cultivars with high resistance and susceptibility to FHB in the group of foreign origin group. For instance, the well-known Chinese highly resistant cultivar Sumai 3 and landrace Wangshuibai were used as resistance checks, while the South African cultivar Gamenya was used as a check for susceptibility. Other exotic genotypes, such as CIMMYT FHB resistant lines SHA3/CBRD and MILAN/SHA7 and the Chinese FHB resistant line N894037, also have well-known resistance and were included as references of resistance. The field experiments were carried out at the Institute of Agriculture, Lithuanian Research Centre for Agriculture and Forestry (LAMMC) in 2022 and 2023. Genotypes were grown in Alpha experimental design in a 1.5 m² plot with two replicates. Spray inoculation method was used in the field and greenhouse. The spikes were covered with plastic bags for 48 h to maintain moisture and were later assessed for type I resistance to FHB. For the spawn grain inoculation method, the genotypes were grown in pots and kept in a nursery where sprinklers were installed to create a humid environment for the plants. To inoculate plants, spawn grains were applied at a rate of 30 g/m² at the stem elongation time, early booting stage, and at the end of the booting stage (Tekle et al., 2018). The same technique was followed for the point inoculation method in the greenhouse, where plants were grown in pots at 20-23°C temperature. When the plants reached the anthesis stage, the inoculum was injected into the central florets through glumes to provide mock disease pressure for each genotype to evaluate them for type II resistance. To conserve moisture, the spikes were covered with plastic bags for two days. The details of the trial design, preparation and concentration of inoculums, inoculation techniques and meteorological conditions were described in detail previously (Syed et al., 2024).

For phenotypic data, disease severity, FHB index and FHB incidence were recorded.

Genotyping of the 335 wheat genotypes performed using 25K SNP chip array at TraitGenetics GmbH (Germany) as described previously (Aleliūnas et al., 2024; Lin et al., 2024) under the framework of the NOBAL wheat project. Raw marker data were obtained from a previous study (Aleliūnas et al., 2024; Lin et al., 2024). Additionally, 35 genotypes which were selected as resistant in one of two cycles of pre-breeding for FHB resistance were included in this study, and genotypic data curation and quality control were performed separately. Monomorphic markers were removed from the dataset along with completely failed markers. Afterwards, the filtering was done to remove the markers with more than 20% missing data and minor allele frequency (MAF) greater than 0.05. After the filtering, a subset of 18,417 high-quality genetic markers were selected for the analysis. Overall, 332 genotypes passed quality control and were included in GWAS.

GWAS was performed using a Bayesian-information and Linkage-disequilibrium Iteratively Nested Keyway (BLINK) model with two the first principal components as covariates (Huang et al., 2019) in R package “GAPIT3” (Wang and Zhang, 2021). To find out the optimal number of principal components (PCs), the Elbow method was applied ( Supplementary Figure S1 ). The scree plot was built using “ggplot2” and “ggrepel” R packages for visualization. The genomic Inflation Factor (lambda) was closer to one for all studied traits when two PCs were included as covariates, confirming that two PCs were the optimal number.

The significance Level (α) was adjusted using the Bonferroni correction, dividing the original significance level of 0.05 by the total number of tested SNPs (0.05/18417 = 0.0000027149) and was converted into -log₁₀ (p) of 5.565 to improve data visualization.

The pairwise LD of SNP markers were calculated separately for each chromosome applying the full-matrix option in “TASSEL 5” software package (Bradbury et al., 2007). The results from TASSEL were used to estimate the LD decay over physical distance (Mbp) in R. The nonlinear model of Marroni et al. (2011) was applied to summarize the relationships between LD decay and physical distance (Marroni et al., 2011). Subsequently, a half-decay distance based on the maximum LD R² value and decay distance using a fixed LD threshold of 0.2 were found. QTL determination: MTAs were grouped as one part of the same QTL region when the LD (R² ) coefficient was greater than 0.5.

GWAS was conducted using the Best Linear Unbiased Estimates (BLUEs) from individual trials and BLUEs across different years, as well as among different types of trials, to assess the stability and genetic pattern of resistance of found significant marker–trait associations (MTA).

To describe the effects of SNPs, we used three characteristics/metrics: allelic effects (β), PVE, and the coefficient of determination R² (adjusted).

Allelic effects (β) were computed using the BLINK model, and PVE was calculated using the following formula:

where:

MAF: Minor Allele Frequency of the SNP.

β: Estimated effect size of SNP.

σ²: Total phenotypic variance of the trait.

The coefficient of determination R² (adjusted) was estimated by fitting the Ordinary Least Squares (OLS) regression model using the lm() function in R. Furthermore, R² (adjusted) values were extracted and converted to percentages by multiplying by 100. The adjusted coefficient of determination is denoted as R² throughout the text.

Population analysis and identification of genotype clusters were performed using Principal Component Analysis (PCA) and hierarchical clustering based on 18,417 polymorphic SNPs. PCA was performed using the GAPIT package and BLINK method with default settings in R. The centering of genotypes was done without scaling using the prcomp function. Missing data was imputed with the middle values (Huang et al., 2019). The PCA plot was drawn using ggplot2 package. Hierarchical clustering was based on the first two principal components of PCA. To perform hierarchical clustering, a distance matrix was computed using the Euclidean distance metric. Hierarchical clustering was performed using the hclust function with Ward’s method, and the cutree function was used to assign genotypes to clusters. Linkage distances were normalized to a percentage scale. The factoextra package was used to draw a dendrogram plot (Kassambara and Mundt, 2020).

Best linear unbiased estimates (BLUEs) were calculated for individual trials and across different combinations of trials using the META-R software (Alvarado et al., 2020). In the linear mixed model, genotypes were fitted as fixed effects, while environments, replicates, and genotype-by-environment interactions were treated as random effects. Broad-sense heritability (H²) was estimated using the following formula:

The details of calculating broad-sense heritability (H²) mentioned in the article (Syed et al., 2024).

Calculating the R² (phenotype explained variation), the average means of the markers were used to replace NA, since replacing them with a major allele could slightly distort the effect of alleles.

To compare the genotypic responses between different clusters, violin plots with boxplots, jitter, and means were created using the ggplot2 package in R. Pairwise Wilcoxon tests were performed to compare the genotypes responses to FHB between clusters. The p-values were adjusted using the Bonferroni method (Hollander et al., 2013).

Multiple regression analysis was performed using the lm() function from the base R package to determine the effects of MTAs/QTL on wheat phenotypic resistance (BLUEs across eight trials). Furthermore, the proportional contributions of major and minor alleles to the explained phenotypic variation were calculated using the LMG (Lindeman–Merenda–Gold) method via the function calc.relimp() from the package relaimpo in R (Grömping, 2007).

The visual displaying of major QTL on wheat chromosomes was built using an online tool “MG₂C” for drawing genetic maps (Chao et al., 2021).

To identify putative candidate genes associated with FHB resistance, a list of genes in the genomic region within the identified QTL was extracted from the published spring wheat “Chinese Spring” RefSeq Annotation v1.1 (International Wheat Genome Sequencing Consortium (IWGSC), 2018) using a custom R script. The GO terms for the identified genes were obtained using the DAVID online tool (Sherman et al., 2022). GO term definitions were downloaded from https://geneontology.org/docs/download-ontology/ (The Gene Ontology Consortium et al., 2023). Combined protein families, domains, and the functional annotations (InterPro descriptions) were collected from EnsemblPlants (https://plants.ensembl.org/) connected to “Plants Genes 60” and “Triticum aestivum Refseqv2 genes” databases.

Principal Component Analysis (PCA) was applied to investigate the population structure and detect the genetic patterns that explained the variation based on allele variation. The first principal component explained the largest amount of variation (9.2%). The first two principal components explained 15.1% of the variation, and the first three principal components explained 20.2%. The Elbow method demonstrated that two PCs capture the major genetic variation in our population and can be applied as covariates to adjust for population structural alleles ( Supplementary Figure S1 ). Apart from that, the genomic inflation factor (λ) was closer to one when two principal components (PCs) were included as covariates, confirming that using two PCs as covariates resulted in sufficient correction for the population structure. Adding a third PC might not improve the findings, but could lead to overcorrection, potentially removing true associations.

The set of studied genotypes consisted of 10 exotic genotypes and 322 genotypes of European origin. To determine the number of clusters within the population, all genotypes were distributed according to their positions on the first two principal components, and the Elbow method was applied. According to the Within-Cluster Sum of Squares in the Elbow method, there were at least three subgroups within our population ( Supplementary Figure S2 ). However, a cut-off point was selected at a 10% level of dissimilarity, assigning the genotypes into 13 clusters/subgroups for a detailed comparison of genetic relatedness between exotic resistant genotypes and adapted resistant genotypes ( Supplementary Figure S3 ).

Analysis of the genetic structure of the population based on 18,417 SNP markers demonstrated that Cluster 12 was the most distinct. Some distinctness was found for Clusters 13, 8, and 11 and, to a lesser extent, for Cluster 3. The other 8 clusters were relatively closely centered, indicating their close genetic relatedness ( Figure 1 ). Among the 13 clusters, 2, 7, 8, 12, and 13 did not comprise any resistant or moderately resistant genotypes, while resistant fell closer to the middle of the population ( Figures 1 , 2 ). Interestingly, exotic genotypes were intermixed with European genotypes and fell into two Clusters, 4 and 5, which were located in the middle of the population ( Figure 1 ; Supplementary Table S1 ).

The most resistant genotypes were found within 4 Clusters: 3, 4, 5, 6, and 10. There was no clear relationship between the origin of genotypes, their resistance and genetic distinctness (based on 18,417 polymorphic SNPs). Resistant genotypes and genotypes of exotic origin did not form separate clusters ( Figure 2 ). The clusters that comprised the most resistant genotypes also contained susceptible genotypes. Lithuanian genotypes prevailed among the most resistant, and few of them were of exotic and Latvian origin ( Figure 2 ). Lithuanian breeding line DS-1577-8-DH does not have exotic relatives in the pedigree, had comparable levels of resistance to Sumai 3 (Syed et al., 2024) and belonged to the large genetic Cluster 10, which comprised both highly resistant and susceptible genotypes. Exotic resistant genotypes such as Wangshuibai, Gen020 (SHA3/CBRD), Gen066 (N894037), and the Lithuanian resistant breeding line DS-1401-6-DH and moderately resistant DS-570-2-DH and DS-1401-6-DH belonged to the same Cluster 4 ( Supplementary Table S1 ). Clusters 3, 6 and 10 did not contain any exotic germplasms but comprised several resistant or highly resistant local genotypes. For example, Cluster 6 contained several resistant genotypes of Latvian origin. Clusters 6 and 10, comprised several resistant genotypes of Lithuanian origin. Genotypes which belonged to Cluster 1 had the highest susceptibility on average. According to the BLUEs values across the eight trials, the resistance of the genotypes in Cluster 1 differed significantly from that in Clusters 4, 6, 10, 11, and 12. Cluster 12 comprised mostly genotypes of Latvian origin with similar level of moderate resistance ( Figure 2 ).

The LD pattern of the 332 spring wheat genotypes is depicted in a scatter plot of pairwise LD (r² ) over physical distance in Mbp ( Supplementary Figure S4 ). The average half-decay distance was 0.57 Mbp (when r² decreased to half of the maximum value) with the longest distance in the B genome (0.76 Mbp) followed by the D and A genomes (0.57 and 0.39 Mbp). The fastest LD decline was observed in chromosomes 2A and 6A (0.24 Mbp) and the longest LD was observed in chromosome 3B (1.03 Mbp) ( Supplementary Table S2 ). The half-decay distances indicate that approximately half a million base pairs may belong to one LD block. Thus, on average, one SNP marker may be associated with about half a million base pairs. The level of genetic recombination was moderate in our population and consistent with other studies conducted on diverse wheat populations (Ghimire et al., 2022).

Significant MTAs were identified across all individual trials and among different combinations of trials. In total, 58 MTAs with known positions and seven with unknown positions were identified ( Supplementary Table S3 ). The lowest number of MTAs was found in the greenhouse experiment, where point inoculation was used, and the largest number (12 MTAs) was found in the combinations of different trials ( Tables 2 , 3 ).

Among 65 SNP markers, 12 explained more than 10% of the phenotypic variation according to the adjusted coefficient of determination (R² ). Among them, 7 MTAs had R² >10% and were stable (were detected in more than 1 trial) ( Supplementary Table S3 ).

MTAs were found across 15 chromosomes (1B, 2A, 2 B, 2D, 3A, 3B, 4A, 4D, 5A, 5B, 5D, 6B, 6D, 7A, and 7D) in the three wheat genomes. Among them, 25 MTAs were found in the A genome, 27 in the B genome, and six in the D genome. Cumulatively, across all the trials, not within a single trial, they explained 338%, 321%, and 37% of the variation (R² ), respectively. The A and B genomes harbored the largest number of MTAs and explained a major part of the variation. Individual SNP analysis was performed for all significantly associated SNPs across the trials, where the associations were significant at the adjusted significance level of 5.565 (-log₁₀). Comparing the BLUE values of groups with different alleles, the Wilcoxon test indicated that there were significant differences for 56 SNPs and no significant differences for nine SNPs. The results are presented in Supplementary Table S3 . The strongest and most stable associations were found between the SNP marker Tdurum_contig91519_224 and FHB resistance. The allele effect of this marker ranged from 4.41% to 6.94%, PVE from 0.97% to 7.18%, and R² from 8.86% to 15.18%. The resistance to FHB was associated with minor allele “T” which had frequency of 17% in our panel, while most genotypes had “G” allele associated with susceptibility (83%). The difference in disease severity between the two groups ranged from 12.2% to 18.7% ( Supplementary Table S3 ). The significant effects of this SNP marker were identified in nine different combinations of trials, demonstrating stable effects. However, a significant effect of this marker was identified in only one individual trial, namely field spray trial in 2023. A visual comparison of the groups with different alleles across eight individual trials and nine trial combinations is shown in Figure 3 . Comparing the average values between groups with “G” and “T” alleles the Wilcoxon test demonstrated a significant effect of alleles in all trials and combinations, even in the trials where significant associations were not detected by GWAS (under a threshold of -log₁₀(p) = 5.565). Significance at the level of p ≤ 0.05 between two groups was found in the following trials: spawn severity (2023), greenhouse point inoculation (2022), and incidence (2023); in other trials, the differences were significant at p ≤ 0.001 ( Figure 3 ).

Up to 12 alleles significantly associated with resistance were found in the combinations when BLUE values were calculated across five, six, and seven trials and from one to eight alleles in the individual trials ( Tables 2 , 3 ). However, many of these alleles provided minor effects ( Supplementary Table S3 ). To select the most meaningful ones and analyze their cumulative effect an R² threshold of 5% was applied. Furthermore, to investigate the pyramiding effect of resistant (R) alleles, the wheat genotypes were divided according to the number of R alleles. The alleles associated with FHB resistance were labeled as “R alleles” and the number near the letter “R” indicates the number of associated alleles ( Figure 4 ). Comparison of the average BLUE values of the groups according to the number of R alleles demonstrated a relationship between the number of R alleles and the level of resistance ( Figure 4 ). The Kruskal-Wallis test (p < 0.05), followed by Dunn’s test for multiple post hoc pairwise comparisons with Bonferroni correction, were conducted among the genotype groups to determine if the differences were significant between groups. Across a two-year field spray trial, four MTAs (four R alleles) with R² effects above 5% were identified, and the level of resistance gradually increased depending on the number of R alleles. The average disease severity mean of the null group was 75.0%, the group with one allele had disease severity of 66.5%, while two, three and four alleles reduced disease severity to 58.6%, 50.1%, and 40.0%, respectively. The group that contained four R alleles was significantly more resistant than the groups with 0, one and two R alleles. The groups with two and three R alleles significantly differed from the groups with 0 and one alleles, and the group with one R allele significantly differed from the group with no R alleles ( Figure 4A ).

Across five trials of overall resistance, four MTAs with R² effects greater than 5% were detected. According to the average means of each group, an increased number of MTAs per genotype resulted in an improved resistance to FHB. The average severity level of the null group was 55.8%, while the severity levels for groups with increasing numbers of alleles were as follows: one allele (49.6%), two alleles (44.6%), three alleles (39.5%), four alleles (30%), five alleles (20.8%). The average means of groups 5R, 4R, and 3R significantly differed from the groups with 0, one, and two R alleles, and the groups 0R, 1R, and 2R significantly differed from each other ( Figure 4B ).

Across all trials, five significantly MTAs with R² effects greater than 5% were identified. Comparing the average means of each group, a clear relationship between the number of accumulated R alleles and disease resistance levels was observed. The average severity of the null group was 60.0%, while the severity levels for groups with increasing numbers of R alleles were as follows: one allele (53.2%), two alleles (46.5%), three alleles (41.0%), four alleles (31.8%), and five alleles (9.8%). The groups with two, three, four and five alleles differed significantly only from the groups 0 and 1R, while the average means of the groups 0 and 1R significantly differed from each other and the groups with two, three, four, and five R alleles ( Figure 4C ). The group that accumulated 5 alleles likely did not differ significantly from the other groups, because it was represented by only one genotype.

In total, 58 significantly associated MTAs were found across different individual trials and combinations of trials. Furthermore, the relationships between the number of these MTAs and BLUE values across eight trials were determined. Additionally, the genotypes were grouped according to the number of MTAs and their geographic origin ( Figure 5 ). Multiple regression analysis was performed to determine the effect of 58 MTAs and partial contributions of 10 major and 48 minor MTAs. Collectively, 58 MTAs explained 55% of the variation in resistance. A larger part of this variation was explained by 10 major MTAs (31%), while 48 minor MTAs were responsible for the remaining 24% ( Supplementary Figure S5 ). Most of these MTAs were unstable, identified in single trials (or across a specific combination of trials), and demonstrated only minor effects ( Supplementary Table S3 ). Nevertheless, a relatively strong correlation (r=0.69, p < 0.001) and linear regression coefficient of determination (adjusted R² = 0.55) were found between the number of resistant alleles and BLUEs across eight trials ( Figure 5 , Supplementary Figure S5 ). The most resistant cultivars derived from Baltic countries, and a few genotypes had an exotic origin ( Figure 5 ). Additionally, the associations were determined only for major MTAs (using a threshold of R² >10%) to determine the overall effect of major MTAs ( Figure 6 ). Ten major-effect MTAs explained the largest part of resistance (r=0.61, p < 0.001; R² = 0.31) ( Figure 6 , Supplementary Figure S5 ). However, considering only major QTL, some genotypes fell into the tails of violins as most resistant, but belonged to the second group possessing only 4–6 MTAs ( Figure 6 ). In contrast, when all MTAs were considered, the most resistant genotypes (in the tails of violins) belonged to the Group 3 with the highest number of MTAs ( Figure 5 ). Overall, 48 minor MTAs explained 24% of the variation, indicating that small-effect MTAs also contributed significantly to resistance.

Across all trials and their combinations, we identified a total of 65 MTAs, 58 of which had known physical positions within the genome ( Supplementary Table S3 ). To define QTL regions associated with FHB resistance, half-decay chromosome-specific distances were applied ( Supplementary Table S2 ). If the LD values were higher than 0.5 between nearby MTAs, they were assigned to one QTL. Overall, 55 QTL were defined ( Supplementary Table S4 ), 10 of which had R² coefficients higher than 10% ( Table 3 ). Among them, 51 were related to overall resistance, and only one QFHB-2D was associated with Type II resistance (R² value was 6.2%). Three QTL were associated with Type I resistance FHB-2B.3, with an R² value of 13%, and QFHB-7A with an R² value of 6%, while the third QFHB-2BS was related to both Type I resistance and overall resistance ( Supplementary Table S4 ).

To assess the correspondence between the identified QTL and previously published ones, we compared their physical positions. For comparison, we used a set of published QTL from 2000 to 2020 compiled by Zheng et al. (2021). The authors collected 625 FHB associated QTL from 113 publications and determined their physical positions by aligning flanking and peak markers to the Chinese Spring reference genome, RefSeq v1.1 (International Wheat Genome Sequencing Consortium (IWGSC), 2018) (Zheng et al., 2021).

Comparing the physical positions of QTL identified in this study with previously published, we found that 46 QTL were collocated with previously reported QTL ( Supplementary Table S5 ), whereas 8 QTL did not match any published QTL, 4 of them (AX-89483131, TA001900-1836, AX-158584923 and IAAV6297) had R² values above 5% and the difference between the two groups of alleles was significant at p ≤ 0.001. AX-95100505 and AX-158585120 had R² values of 2.87% and 2.65%, respectively. AX-111556997 and TA016804–1075 had R² values close to 0% and there was no significant difference between the two groups of alleles. Meanwhile, all 10 QTL with phenotypic variation (R² ) higher than 10% collocated with previously reported QTL ( Table 4 ).

Additionally, we reviewed recent GWAS that reported associations with wheat resistance to FHB. Overall, 1760 MTAs were found in 13 articles ( Supplementary Table S5 ). To evaluate the potential collocation of our QTL with previously published MTAs, we manually added ±20 Mb to the positions of the markers. It was found that 42 QTL were potentially related to the MTAs from previous GWAS studies. When both published QTL and MTAs were compared with the QTL from this study, only 2 QTL (QFHB-4DS and QFHB-7DS.1) were not collocated in physical positions and are presumably novel ( Supplementary Table S7 ). QFHB-4DS was found only in one individual trial (spawn 2022) with an R² value of 7.68%, and the difference between the two groups of different alleles was 11.67% at the 0.001 significance level, while QFHB-7DS.1 was found in one individual trial (FHBindex_2023) with an R² value close to 0%, and no significant difference between the two groups was found.

Using annotated wheat genome data for high-confidence genes from the International Wheat Genome Sequencing Consortium (IWGSC) RefSeq v1.1 (International Wheat Genome Sequencing Consortium (IWGSC), 2018), 705 gene names and their physical positions were extracted from the physical regions of 55 identified QTL. Furthermore, 12590 matching GO terms were found for 543 genes ( Supplementary Table S8 ). Alongside, significantly up- or downregulated genes from previous transcriptomic studies were collected, in which gene expression was analyzed by comparing resistant and susceptible genotypes/groups in response to inoculation with Fusarium graminearum, to assess their correspondence with identified genes in this study (Buerstmayr et al., 2021; Pan et al., 2018; Seifi et al., 2023). Subsequently, 52 genes could be referred as candidate genes, since they matched both our identified genes and those reported in transcriptomic studies as significantly associated with FHB. These 52 genes were found within 25 QTL identified in our study ( Supplementary Table S9 ). Most genes were downregulated, 36 exhibited significantly decreased expression, while only 16 were upregulated, which is consistent with findings from previous studies (Buerstmayr et al., 2021; Seifi et al., 2023). Notably, several genes were found within the same QTL. For example, eight candidate genes were found within QFHB-1BS.2, seven within QFHB-1BS.1, six within QFHB-1BS.3, four within QFHB-3BL.4, three within QFHB-4DS and QFHB-2BL.1, and two within QFHB-2AL.5, QFHB-2DL, and QFHB-3AL.5 ( Supplementary Table S9 ).

Some candidate genes were inside or within close proximity to SNP markers, such as TraesCS1B02G324300, TraesCS1B02G046300, TraesCS3A02G527200, and TraesCS1B02G075200, which were merely 122, 1129, 1452, 1787 bp away from the SNP markers ( Supplementary Table S9 ). The identified candidate genes encoded proteins that play key roles in pathogen detection and defense signaling. For example, TraesCS3B02G598200, TraesCS3B02G598400, TraesCS1B02G047000 and TraesCS1B02G324300 contribute to cell wall reinforcement. TraesCS5A02G367100, TraesCS2B02G536500, TraesCS1B02G046000, TraesCS1B02G046300, TraesCS1B02G114100, TraesCS6B02G052400, TraesCS3A02G530200 and TraesCS2B02G296300 participate in pathogen recognition that can be referred to as Type I resistance. TraesCS1B02G114300, TraesCS1B02G113600, TraesCS1B02G113800, and TraesCS5B02G426400 were associated with regulation of ROS and programmed cell death that hints to Type II resistance. Other genes were associated with Type III resistance, such as DON detoxification (TraesCS3B02G598200 TraesCS3B02G598400, TraesCS1B02G047000, TraesCS1B02G114300). Moreover, some genes (TraesCS5A02G367100, TraesCS1B02G114100, TraesCS3B02G596200) play putative roles in susceptibility pathways, such as calcium signaling or receptor-like kinases (RLKs). Functional description of all the identified genes are presented in Supplementary Table S9 .

Fusarium head blight (FHB) is a challenging threat to wheat production worldwide. Resistance to FHB is regulated by a complex network of small-effect genes (Buerstmayr et al., 2021). In previous studies, plant breeders and researchers relied on genetic tools that enabled them to work with a limited number of genes or QTL. Traditional genetic methods were limited to screening large genomic regions, which is crucial for investigation of complex traits such as FHB resistance. However, after the development of high-density genotyping techniques such as high-density SNP arrays and genotyping-by-sequencing (GBS) combined with Linkage Disequilibrium analysis, research on FHB resistance has gained momentum, opening new possibilities to develop improved disease resistance in wheat (Mir et al., 2023; Myles et al., 2009; Saini et al., 2021). Using LD analysis, which can be performed directly in a diverse breeding population, a large part of the wheat genome can be quickly analyzed for associations to FHB resistance without the need for the whole-genome resequencing. Furthermore, physical regions that include the responsible genes can be defined by performing LD analysis. Although LD does not provide direct information about genetic linkage, it provides statistical evidence of linkage between alleles based on their actual distribution. That is not completely random, the frequency can be affected not only by genetic linkage but also by genotype pedigree, adaptability to specific environments, and breeding selection (Flint-Garcia et al., 2003; Sahoo et al., 2022; Scherer and Christensen, 2016). The use of LD strength between markers in diverse populations provides quick LD-based detection of QTL regions; however, the physical regions are not fine-mapped, do not have precise boundaries, and are based on selected statistical probability levels.

Traditionally, QTL identification is performed using mapped populations such as recombinant inbred lines (RILs), near-isogenic lines (NILs), or doubled haploids (DHs) (Zheng et al., 2021). Alternatively, GWAS, which is based on high-density genotyping of diverse populations, is usually used to identify MTAs without explicitly defining the specific QTL regions. In GWAS, when genetically diverse populations are used, the structural alleles can increase the risk of detecting false associations (Mir et al., 2023). When mapping populations in which all progenies are genetically related are used for analysis, the detection of associations between markers and traits is more reliable. However, the development of mapping populations requires significant time and effort. Because of the lower frequency of recombination events in such populations, the linkage disequilibrium decay is slower, and defined QTL regions are relatively large. Conversely, the application of diverse populations makes it possible to use the existing breeding populations. The higher frequency of recombination events in diverse populations results in higher LD decay. In our study, LD-based detection of QTL region was used. The main purpose of defining approximate QTL regions based on chromosome-specific half-decay distance was to set the searching boundaries for candidate genes and to check the collocation between relatively small QTL identified in diverse populations and previously detected QTL from mapping populations.

Overall, 58 MTAs with known positions and seven with unknown positions were identified in the studied wheat panel. A comparable number of identified MTAs was found in previous GWAS studies of FHB resistance in wheat, when diverse populations were studied (Ghimire et al., 2022). In contrast, mapped biparental populations do not possess such diversity of favorable MTAs/QTL, and the number of identified QTL is usually several times smaller in biparental populations than in diverse populations (Zheng et al., 2021). Applying half-decay chromosome-specific distances, 55 QTL were determined in the spring wheat population. Three pairs of MTAs were merged because they were strongly associated, LD values between them varied from 0.94 to 1.0. By comparing the physical collocation of the QTL identified in this study with previously detected QTL from mapped populations using chromosome-specific half-decay distances, we surprisingly found that 46 out of 55 were collocated with known QTL. This might indicate that the significantly associated markers were relatively closely positioned to the causal genes, and the LD approach of defining boundaries of QTL can be applied for approximate determination. Almost all identified QTL were associated with overall resistance. Only QFHB-2D (AX-94872625) was associated with Type II resistance after point inoculation under controlled conditions, with an R² value of 6.2%. Three QTL associated with Type I resistance: QFHB-2B.3 (AX-94456169) with R² value of 13%, and QFHB-7A (AX-158591608) with R² value of 6%, which are presumably related to the initial regulation of wheat defense. The third QFHB-2BS (BobWhite_c8113_532) was connected not only with initial resistance but was associated with later resistance and found in ten different combinations of trials ( Supplementary Table S3 ). The remaining 51 QTL were related to overall resistance. Crossing wheat cultivars that possess different types of resistance might be an effective strategy to combine the different components of resistance. For example, the breeding line Gen323 (DS-1401-6-DH) possesses three QTL of Type I resistance and nine QTL of overall resistance (across five trials of overall resistance) but does not contain QTL of Type II resistance. Crossing this breeding line with other breeding lines that already possess Type II resistance allele, along with a complex of overall alleles (e.g. breeding line Gen322), might result in the development of transgressive segregates with improved FHB resistance.

Some minor QTL from this study collocated with QTL with relatively high effects in other studies. For example, QFHB-2DL had R² value of 6.18 in our study, while in other studies, its effect ranged from 12.6 to 18 ( Supplementary Table S5 ). It is noteworthy that the most stable and effective QTL, QFHB-2AL.1, collocated with only one published QTL by Zhang et al. (2014). However, the authors defined a large region from 341.9 Mbp to 534.9 Mbp indicating 193 Mbp length/window at 2A chromosome and used mapping population, while in our study we defined regions from 451.6 Mb to 452 Mbp using only 0.48 Mbp window. Moreover, no genes related to FHB resistance were found within QFHB-2AL.1 and QFHB-3AL.1, according to the GO terms ( Supplementary Table S8 ) and InterPro descriptions ( Supplementary Table S10 ). Presumably, resistance genes were located outside these defined QTL. When the region of QTL is defined based on LD, the causal genes are more likely to be within this region; however, they can still be out of the defined region. Therefore, additional analysis of QTL regions is required.

In total, 705 genes were extracted from the physical regions of 55 identified QTL. Matching significantly up- and downregulated genes from the transcriptomic studies (Buerstmayr et al., 2021; Pan et al., 2018; Seifi et al., 2023) with genes located within determined QTL, 52 candidate genes were found. Among them, 36 (69.2%) were downregulated and 16 (30.8%) were upregulated in response to inoculation with FHB ( Supplementary Table S9 ), which is consistent with the findings of previous studies (Buerstmayr et al., 2021; Seifi et al., 2023). We did not find any functional validation of these 52 candidate genes, such as gene silencing or overexpression, in published literature to prove the true relationship between these genes and FHB resistance. However, some identified candidate genes encode wheat proteins, the roles of which have been well studied in wheat defense against FHB, supporting their potential contribution to FHB resistance (Ma et al., 2022; Sirangelo, 2024). For instance, TraesCS3B02G598200 encodes a glycosyltransferase that plays a vital role in DON detoxification by modifying it into DON-3-glucoside, cell wall reinforcement, and hormone regulation (Ghimire et al., 2022; Gottwald et al., 2012; He et al., 2018; Lemmens et al., 2005). TraesCS1B02G114300 encodes Glutathione S-transferase (GST) that is also a well reported protein for FHB resistance, associated with defense activation, detoxification, and oxidative stress (ROS). Additionally, it is an antioxidant that contributes to minimizing PCD (Gullner et al., 2018) by forming DON-glutathione conjugates, and aids in DON detoxification (Edwards et al., 2000; Gardiner et al., 2010; Seifi et al., 2023). TraesCS5A02G367100 and TraesCS2B02G536500 are encoded for receptor-like protein kinases (RLKs). These proteins are paramount for plant growth, strengthening of the cell wall, oxidative stress responses, and activation of defense signaling. Protein kinases play a crucial role in signaling when a pathogen is identified (Romeis, 2001; Tang et al., 2017). Furthermore, TraesCS7A02G051800 is identified as a 2-oxoglutarate/Fe(II)-dependent oxygenase (2OG oxygenase) superfamily enzyme involved in secondary metabolite synthesis (Jia et al., 2017). Inspection of the functions of candidate genes might provide additional insights into their resistance mechanisms and facilitate the proper pairing of QTL to attain durable FHB resistance. In our analysis, it was found that some QTL contributed to basal defense, while others participated in the detoxification of FHB mycotoxins ( Supplementary Table S9 ). Seifi et al. (2023) has previously postulated that if PCD is inhibited but DON production is not limited, the visual symptoms may be reduced, nevertheless, a high level of mycotoxins would be present. Conversely, if DON production is limited via the activation of ROS scavengers, and plant cell death process is not limited, visual symptoms would be visible; however, the levels of DON would be low in the grain. Therefore, breeders should aim to develop such resistance when both DON production and PCD are limited (Seifi et al., 2023).

The results of the current study confirm the commonly accepted statement that wheat resistance to FHB is controlled by multiple number of small-effect QTL, and the determination of which is hampered by environmental effects (Mesterhazy, 2024; Miedaner et al., 2024). In this study, the GWAS was conducted using phenotypic data obtained from different inoculation methods under field and controlled conditions. Moreover, GWAS was carried out not only within individual trials but also across different combinations of trials that extended the number of identified significantly associated QTL. Altogether, 55 QTL were identified, most of which had small effects and were found only in one individual trial or in one combination of trials ( Supplementary Table S3 ). The association between the number of QTL and FHB resistance was relatively strong, considering that FHB resistance is highly complex, some QTL related only to one type of resistance (Type I, Type II or overall), and were determined only under specific inoculation methods or environments. A relatively strong Pearson correlation (r=0.69, p < 0.001) and linear regression coefficients of determination (adjusted R² = 0.55) were found between the number of resistant alleles and the BLUE values of resistance across eight trials ( Figure 5 ). Consequently, the effect of 10 major QTL was 31% (R² = 0.31), while the effect of 45 minor QTL was 24% (R² = 0.24) on the phenotypic variation ( Supplementary Figure S5 ). The results demonstrate that even not stable “hardly detected” minor QTL are essential for wheat resistance to FHB and should not be neglected. The range of accumulated favorable QTL varied from 10 to 40 per genotype, while the number of major QTL varied from 0 to 10. The results demonstrate that these major QTL explained a large part of the variation and could potentially facilitate wheat breeding for FHB resistance. Additionally, this might indicate that modern cultivars possess plenty of favorable small-effect QTL, and new superior genotypes can be developed due to the transgressive segregation using selection under high disease pressure. The development of improved FHB resistance through transgressive segregation was observed by Mesterházy et al. (2018) and Zhang et al. (2018). Furthermore, the findings indicate that the inheritance of resistance to FHB had an accumulative nature. The more favorable QTL stacked, the greater the resistance to FHB. The positive additive effect of pyramiding small-effect QTL for improving wheat resistance has been reported in many studies (Haile et al., 2023; Kirana et al., 2023). For example, in the study by Ghimire et al. (2022) investigating the resistance of 278 elite winter wheat breeding lines, 11 QTL with R² greater than 10% were detected using GWAS. The severity of FHB was reduced by 21.6% when breeding lines possessed three alleles associated with FHB, a combination of four alleles resulted in a reduction of 32.6%, and a combination of five alleles caused a reduction of 42.6% compared to the null group (Ghimire et al., 2022). However, it should be noted that in this study some genotypes possessed close to the maximum number of favorable alleles, whereas they had an average or even below-average level of resistance ( Figure 5 ). These results might suggest that there were other factors influencing resistance, such as phenotyping inaccuracies caused by environmental effects, the presence of additional unidentified FHB-related QTL, epistatic gene interactions, or epigenetic regulations that contributed to the expression of FHB resistance (Brar et al., 2019; ElDoliefy et al., 2024; Kumar et al., 2020; Mierziak and Wojtasik, 2024).

The largest limitations in wheat breeding for FHB resistance are the lack of adapted sources of resistance and the genetic linkage between resistance and undesirable agronomic traits. The widely utilized Chinese sources of resistance, such as Sumai 3 and Wangshuibai, possess inferior agronomic traits and do not perform well in the environments of the Baltic and Nordic countries. Therefore, the search of resistance for ‘native’ resistance within locally adapted germplasm is an alternative breeding strategy for FHB resistance (Brown-Guedira et al., 2008; Mesterházy et al., 2018; Thambugala et al., 2020). Phenotyping a panel of 332 spring wheat varieties and breeding materials demonstrated that FHB resistance can be found in adapted genotypes to these regions (Syed et al., 2024). The medium resistance was found in breeding materials from general breeding and the most resistant breeding lines were detected in a subset of resistant breeding lines (breeding lines from Gen303 to Gen336) (Syed et al., 2024) which was previously selected after one or two cycles of screening under artificial inoculation. Some of these lines had resistance comparable to that of Sumai 3 (Syed et al., 2024). In our study, GWAS revealed ten major favorable alleles in wheat genotypes adapted to the environments of Baltic and Nordic countries. Although most of the identified QTL were found within previously published genomic regions, it was demonstrated that resistance to FHB can be found in the European gene pool, and moderate to high levels of resistance are controlled by the additive effect of major and small-effect QTL. A GWAS conducted on a large collection of adapted genotypes to Baltic and Nordic countries enabled identification of QTL, which presumably do not possess a negative linkage with crucial agronomic traits because they were selected within adapted breeding material. The accumulation of QTL identified in this study and the utilization of breeding lines carrying resistance alleles will facilitate further genetic improvement of wheat resistance to FHB in the environments of the Baltic and Nordic countries. These findings demonstrate that FHB resistance in wheat cultivars for these regions can be further improved by pyramiding these R alleles.

In conclusion, 55 QTL for wheat resistance to FHB were identified in the European gene pool. Among 55 QTL associated with FHB resistance, 10 had major effects (R² >10%), two QTL (QFHB-4DS and QFHB-7DS.1) have not been previously published and are presumably novel. The majority of QTL (52 QTL) were associated with overall resistance, one QTL (QFHB-2DL) was linked with Type II resistance, and two QTL (QFHB-7AS and QFHB-2BL.2) were associated with Type I resistance. According to the principal component analysis performed on the basis of 18,417 SNPs, resistant genotypes did not fall into separate clusters, but were placed in Clusters 3 and 4 in the middle of the population. Cluster analysis demonstrated that exotic resistant genotypes did not form a separate cluster, resistance was fairy well distributed within the population, and it can be found within the European wheat gene pool. Pyramiding of three and more major QTL resulted in improved resistance. A total of 52 candidate genes were identified by analyzing genes near significantly associated SNPs in combination with published transcriptome data. All identified QTL were found in elite, adapted breeding materials, and can be utilized in wheat breeding for improving FHB resistance in the Baltic and Nordic regions.