Characterization of the Genetic Architecture for Fusarium Head Blight Resistance in Durum Wheat: The Complex Association of Resistance, Flowering Time, and Height Genes

Durum wheat is an economically important crop for Canadian farmers. Fusarium head blight (FHB) is one of the most destructive diseases that threatens durum production in Canada. FHB reduces yield and end-use quality and most commonly contaminates the grain with the fungal mycotoxin deoxynivalenol, also known as DON. Serious outbreaks of FHB can occur in durum wheat in Canada, and combining genetic resistance with fungicide application is a cost effective approach to control this disease. However, there is limited variation for genetic resistance to FHB in elite Canadian durum cultivars. To explore and identify useful genetic FHB resistance variation for the improvement of Canadian durum wheat, we assembled an association mapping (AM) panel of diverse durum germplasms and performed genome wide association analysis (GWAS). Thirty-one quantitative trait loci (QTL) across all 14 chromosomes were significantly associated with FHB resistance. On 3BS, a stable QTL with a larger effect for resistance was located close to the centromere of 3BS. Three haplotypes of Fhb1 QTL were identified, with an emmer wheat haplotype contributing to disease susceptibility. The large number of QTL identified here can provide a rich resource to improve FHB resistance in commercially grown durum wheat. Among the 31 QTL most were associated with plant height and/or flower time. QTL 1A.1, 1A.2, 3B.2, 5A.1, 6A.1, 7A.3 were associated with FHB resistance and not associated or only weakly associated with flowering time nor plant height. These QTL have features that would make them good targets for FHB resistance breeding.


INTRODUCTION
Fusarium head blight (FHB), also known as scab and mainly caused by Fusarium graminearum Schwabe [teleomorph: Gibberella zeae (Schwein.) Petch] (Bai and Shaner, 1994;McMullen et al., 1997), is a devastating fungal disease of smallgrain cereals including durum and common wheat and barley, resulting in severe yield and quality losses (Gilbert and Tekauz, 2000;McMullen et al., 2012). Moreover, as food for humans and feed for animals, FHB infected grain also creates health risks due to contamination with mycotoxins. This is a particular concern for durum wheat, as its main purpose is for human consumption (Bai and Shaner, 2004;Zhao et al., 2018;Haile et al., 2019;He et al., 2019). Canada is the largest producer and exporter of durum wheat supplying more than a half of the world's total exported durum (International Grains Council, 2020). Since the early 1990s, FHB has become the major disease threatening durum production in Canada and has caused major economic losses for producers (Gilbert and Tekauz, 2000). In 2016, a severe FHB epidemic caused 65% of the common wheat and 36% of the durum wheat to be downgraded in Saskatchewan, Canada, with an estimated economic loss of $1 billion (Canadian Grain Commission, 2017). It is therefore a priority to develop durum wheat with desirable FHB resistance to protect it from losses.
Currently, the combination of agronomic and chemical control along with genetic resistance is the most effective means to manage FHB (Gilbert and Haber, 2013;Prat et al., 2014). Genetic resistance is preferred due to its lower cost, higher efficacy, and environmental benefit (Prat et al., 2014). Genetic resistance to FHB in wheat is quantitative in expression due to control by multiple minor genes. FHB resistance is also significantly affected by environment (Bai and Shaner, 2004;Buerstmayr et al., 2009Buerstmayr et al., , 2019, thus having lower to moderate heritability (Van Sanford et al., 2001). Therefore, when visual assessment of FHB is performed in the field, lines must be tested in multiple independent environments with intensive phenotyping to reliably identify QTL for resistance.
Fusarium head blight resistance can be categorized into three main types or components: (1) type I -resistance to initial infection measured by the incidence of disease in the presence of natural or augmented artificial inoculum (e.g., spray inoculation); (2) type II -resistance to fungal spread measured by the severity of disease; and (3) type III -resistance to the accumulation of the toxin deoxynivalenol (DON) in infected spikes (Miller et al., 1985;Mesterhazy, 1995;Bai and Shaner, 2004). Till now, more than 556 QTL contributing to FHB resistance have been identified on all 21 chromosomes of hexaploid wheat (Buerstmayr et al., 2009;Liu et al., 2009;Löffler et al., 2009;Venske et al., 2019). These QTL can be refined largely into 56 clusters by meta-QTL analysis (Venske et al., 2019). In spite of the relatively large number of identified QTL for FHB, only three QTL, Fhb1 on chromosome arm 3BS Liu et al., 2006), Qfhs.ifa-5A on 5AS (Fhb5) (Buerstmayr et al., 2002;Somers et al., 2003;Steiner et al., 2019a) and Fhb2 on 6BS Cuthbert et al., 2007) have been validated. All of these resistance loci originate from the Chinese cultivar Sumai 3, which displays among the highest levels of FHB resistance observed (Buerstmayr et al., 2009). Fhb1 is the best validated, and most frequently studied and deployed resistance QTL (Buerstmayr et al., 2019). It is currently the only resistance QTL confirmed to be present in several new FHB North American and European varieties with strong resistance (Hao et al., 2019). Fhb1 is reported primarily as conferring strong Type II resistance, and accounting for 20-60% of phenotypic variation in breeding populations (Miedaner and Korzun, 2012). Fhb1 was recently claimed to be cloned by three research groups as two different candidate genes (Rawat et al., 2016;Li et al., 2019;Su et al., 2019) with conflicting interpretations, leaving room for independent validation.
As hexaploid wheat has significantly more sources of FHB resistance, introgression of resistance from hexaploid into durum wheat is one possible way to expand the durum resistance gene pool. Previous attempts to introgress FHB resistance from Sumai 3 into durum were largely unsuccessful (Prat et al., 2017). However, several recent successes have been reported with Fhb1 from Sumai 3 (Giancaspro et al., 2016;Prat et al., 2017) as well as a non-Sumai 3-related FHB resistance sources (Chu et al., 2011;Zhao et al., 2018). Despite these partial successes, no commercial durum cultivars with QTL from these non-adapted sources have been released due to the lengthy breeding process, linkage drag or suppression of resistance in durum backgrounds. Because of these challenges, utilizing the FHB resistance already present in durum cultivars is gaining favor as a promising approach to bring durum wheat cultivars with improved resistance to market more quickly. Durum cultivars with an improved level of FHB resistance have been developed and released by the North Dakota durum breeding program using this strategy (Zhang et al., 2014). With the same approach, recent durum cultivars, including Brigade (Clarke et al., 2009) and Transcend (Singh et al., 2012) with a better level of FHB resistance have also been successfully developed and released by Canadian durum breeding programs selecting for reduced symptoms in FHB nurseries. Regardless of this initial success, there is still a need to know and identify additional native sources of resistance as well as more exotic sources. Understanding the association of FHB resistance with developmental traits, flowering time and plant height is also important for recommending which resistance loci may be most relevant to a particular breeding program.
Genome wide association studies (GWAS) are a promising way to detect FHB resistance QTL present in diverse genetic sources. Only a few GWAS have been conducted on FHB resistance, including winter wheat (Wang et al., 2017), elite Chinese wheat (Zhu et al., 2020), durum breeding panels (Steiner et al., 2019b) and type II FHB resistance durum diversity panels (Ghavami et al., 2011). In this study, we aimed to use GWAS to explore FHB resistance of domestic durum cultivars and breeding material as well as exotic sources of resistance, including Sumai 3 and emmer wheat introgression lines. With GWAS in multiple environments, we aimed to: 1) explore and characterize FHB resistance QTL in durum wheat from the domestic as well as exotic sources, and 2) identify resistance QTL that colocalize with flowering time and plant height.

Plant Materials
In total, 186 diverse durum wheat (Triticum turgidum L. ssp. Durum (Desf.) Husn.) lines were selected to constitute a durum association mapping (AM) panel targeted to improve FHB resistance in durum wheat. This panel was primarily composed of durum from Canada, including elite Canadian cultivars, advanced breeding lines, recently developed germplasm from Canadian breeding programs and from research projects (Supplementary Table 1

Phenotyping
Lines of the durum AM panel were evaluated for FHB infection in Morden and Brandon, MB, Canada in 2015 to 2017 with artificial inoculation and Indian Head, SK, Canada in 2015 and 2016 with natural infection. At both Morden and Brandon, FHB nurseries, corn spawn inoculum of Fusarium graminearum was used. Corn spawn consisted of grains that were inoculated with a mixture of two F. graminearum isolates, a 3-acetyl-deoxynivalenol (3ADON, M9-07-01) and a 15-acetyl-deoxynivalenol (15ADON, M1-07-02) isolate, after which colonized kernels were air dried. In Morden, approximately 2-3 weeks prior to heading, the corn spawn inoculum was spread at 8 g per single meter row with two applications at weekly intervals. Plots were irrigated three times per week using Cadman Irrigation travelers with Briggs booms. At Brandon, the corn spawn inoculum was applied between the rows at a rate of 40 g/m 6 weeks after planting, with a second application performed at the same rate 2 weeks after the first. Plots were irrigated three times per week with a mist irrigation system to create favorable conditions for F. graminearum infection. In Indian Head, FHB was achieved solely by natural disease infection. FHB incidence (INC, percentage of spikes showing symptoms) and severity (SEV, average percentage of spike with visual symptoms of infection) were estimated with visual assessment. FHB index (IND) was calculated with the formula: (INC × SEV)/100. Plant height (HT) and days to anthesis (DTA) were also recorded for Morden plots.

Genotyping
Genomic DNA of the durum AM panel was extracted from freeze-dried fresh leaf tissue of seedlings with a CTAB based protocol carried out on an automated AutoGen DNA isolation system (AutoGen, Holliston, MA). DNA was quantified with a Quant-iT TM PicoGreen R dsDNA Assay Kit (Thermo Fisher Scientific Inc., Bartlesville, OK, United States) and diluted to 50 ng/µL for SNP array genotyping. Genotyping of DNA was performed with the Illumina iSelect 90K SNP array (Wang et al., 2014) according to the manufacturer protocol (Illumina). SNP arrays were scanned with an Illumina HiScan. Raw intensity files from the HiScan were imported into GenomeStudio Version 2013 (polyploid clustering module v1.0.0, Illumina). SNP calling was performed with the method described by Wang et al. (2014) with 3 cluster steps of the cluster algorithm DBSCAN then OPTICS. All SNPs were subsequently visually checked, and incorrectly clustered SNPs or SNPs with more than 4 clusters were manually removed. Finally, SNPs with minor allele frequency (MAF) below 0.05, and missing genotypes higher than 15% were filtered out. This resulted in a total of 6900 high quality polymorphic SNPs of which 5933 markers were anchored to the wheat consensus map (Wang et al., 2014) for the downstream genome wide association analysis.

Statistical Analysis
Statistical analysis was performed with R 3.4.2 (R Core Team, 2017) with the lme4 package (Bates et al., 2015). Phenotypic traits from each disease nursery site across multiple years were fitted with the linear mixed model (Bates et al., 2015). The model is implemented as: P iy = µ + G i + E y + (G i XE y ) + E iy , where, P iy are the values of the tested phenotypic trait, µ is the population mean, G i is the effect of genotypes, E y is the effect of environments (here, by Year), E iy is the residual, where i is the genotype, y is the year. The restricted maximum likelihood (REML) method within lme4 (Bates et al., 2015) was used to estimate the variance components of each trait. The broad sense heritability (H 2 ) was estimated with the equation across multiple years in each disease nursery site, where: δ 2 G is the genotypic variance, δ 2 GXE is the variance of interaction between genotype and year, δ 2 e is the error variance, y is the number years, and p is the total number of replications in all tested years. The least squares means were used for trait correlation and association mapping analysis. The correlation coefficients of disease response, plant height and days to anthesis across multiple years and multiple sites were calculated with the Pearson correlation test and visualized with the R package "corrplot" (Wei and Simko, 2017).

Linkage Disequilibrium, Population Structure, and Kinship Analysis
Linkage disequilibrium (LD) was estimated by correlation coefficient analysis and used the squared correlation coefficients (r 2 ) for all 5933 anchored SNP markers implemented in Tassel v.5.5.0 (Bradbury et al., 2007). The r 2 values of unlinked genetic markers (defined as genetic distance > 30 cM) were square-root transformed into a normal distribution. The baseline (or critical) r 2 value, a value that suggested LD was likely caused by genetic linkage, was determined by taking the 95% percentile of this distribution (Breseghello and Sorrells, 2006). The scatter plot of r 2 versus genetic distance (cM) was fitted using a non-linear model described by Remington et al. (2001) that was implemented in software PopLDdecay .
Population structure of the durum AM panel was determined with STRUCTURE v2.3.4 (Pritchard et al., 2000) with a pruned SNP marker dataset that was generated with the LD (linkage disequilibrium)-based pruning approach implemented in the software PLINK (Purcell et al., 2007). A total of 2306 pruned markers with LD (r 2 ) ≤ 0.2 were used for population structure analysis. STRUCTURE analysis was performed with a 50000 burn-in length and 100000 Markov chain Monte Carlo (MCMC) iterations from K = 2 to K = 12 (K, specialized clusters of the AM panel). Fifteen independent STRUCTURE runs were conducted for each specialized K. The optimal cluster (K) was determined by the K method (Evanno et al., 2005), implemented in the software Structure Harvester (Earl and vonHoldt, 2012). Independent runs of the optimal K were summarized using CLUMMP v1.1.2 software (Jakobsson and Rosenberg, 2007). The CLUMMP generated Q matrix was used to graph the population structure using Structure Plot software (Ramasamy et al., 2014) and perform downstream GWAS analysis. A phylogenetic tree was built with the neighbor-joining (NJ) method in MEGA6 (Tamura et al., 2013) and visualized with Figtree v1.4.4 1 . Principal component analysis (PCA) was performed with R package Genome Association and Prediction Integrated Tools (GAPIT) (Lipka et al., 2012;Tang et al., 2016).

Genome-Wide Association Study
Association mapping was performed on the durum association panel using the phenotypic data collected from the multiple nurseries in multiple years, including HT, DTA, INC, SEV and IND. Association mapping was performed using 5933 mapped SNPs that had a MAF > 0.05 using both Tassel v5.5.50 (Bradbury et al., 2007) and GAPIT (Lipka et al., 2012;Tang et al., 2016). Different association models were tested in both software packages, and QQ-plots generated from all the models were compared to select the model that best controls false positives and negatives. All the data presented here were generated in TASSEL using a mixed linear model (MLM) incorporating the STRUCTURE (Q) matrix as a fixed factor and the kinship (K) matrix as a random factor (Q + K MLM). To be considered a QTL in this dataset, we selected SNPs that were significant 1 http://tree.bio.ed.ac.uk/software/figtree/ (p < 0.05, marker-wise) in at least four of the tested environments for FHB resistance or two for plant height and flower time, and with at least one environment with a highly significant response (p < 0.001). Significant SNPs on the same linkage group were grouped into a QTL region if markers were linked with LD > 0.2.

RESULTS
Population Structure and Linkage Disequilibrium (LD) Analysis STRUCTURE analysis, principal component analysis (PCA) and NJ-phylogenetic tree analysis were all used to determine clustering of lines within the durum AM panel, and two subpopulations were consistently indicated, as shown by different colors in Figure 1. Subpopulation 1 (shown as green in Figure 1 and Supplementary Table 1) contained 124 lines, and consisted of a large proportion of Canadian cultivars and inbred lines including the older cultivar Kyle, more recent cultivars Strongfield and currently most popular cultivars as Brigade, Transcend and CDC Credence. Subpopulation 2 included 62 lines (shown in red in Figure 1 and Supplementary Table 1), consisting of the founder landrace Pelissier and the majority of lines from Austria. All of the inbreeding lines derived from introgression of FHB resistance genes from Sumai 3 into European durum wheat cultivars were contained in subpopulation 2, as were the majority of T. dicoccoides introgression lines. The baseline critical threshold r 2 value of LD was identified as 0.2, corresponding to a genetic distance around 3.0 cM from the whole genome analysis (Supplementary Figure 1).

Phenotypic Analysis
Mean values (across years) of FHB INC, SEV, IND, DTA and HT of lines from the durum AM panel at Brandon, Morden, and Indian Head, were summarized in Supplementary  Table 1. Across environments, FHB INC tended to be higher than SEV (Figure 2) which is reflected in the overall means ( Table 1). The lowest INC was observed at Indian Head in 2016, the location with the lowest severities in both 2015 and 2016. Moderate SEV were observed at Brandon in 2016 and 2017. Generally, a large differential in FHB INC and SEV was observed as indicated by the range for each environment in Table 1, except Indian Head where the maximum severity of disease was less than 100%. Plant height showed a larger range with the average shortest 55 cm and the highest 148 cm while DTA was observed in a range of 13 days in 2017 and 20 days in 2015 (Table 1 and Supplementary  Figure 2). For both INC and SEV, moderate to high broad sense heritability was observed with the two sites under artificial inoculation showing lower heritability than the natural infection site (Table 1). HT showed the highest heritability, while DTA had the lowest heritability (Table 1). For FHB INC and SEV, moderate to high correlations were observed in all tested environments (years and sites). Generally, both HT and DTA had very significant negative correlations with INC and SEV (Figure 3). Analysis of variance (ANOVA) revealed that genotypic effects were significant for all phenotypic traits (P < 0.001, Supplementary Table 2).

GWAS Analysis of FHB Resistance, HT and DTA
With GWAS analysis, 31 genomic regions were significantly associated with FHB resistance traits (Figures 4, 5). The quantilequantile (QQ) plots (Supplementary Figure 3) showed that, for the majority of traits, an appropriate model was fitted for the GWAS test. The GWAS results were summarized in Table 2 and  Supplementary Table 3. SNPs located within the same region were grouped into QTL, and Table 2 shows the QTL names and physical location of the associated SNPs based on their location on the IWGSC Chinese Spring (CS) reference 1.0 (CS Ref 1.0; International Wheat Genome Sequencing Consortium [IWGSC], 2018). For each significant QTL, the lowest -log10 (pvalue) is shown for each environment and trait tested whenever the p-value is less than p = 0.05. As shown in Table 2, there was significant variation in detection of QTL across all of the environments, and more detection of INC than SEV across the environments. The majority of the FHB resistance QTL colocalized with DTA and/or HT.
A major QTL, 1B.1, was found between 544 and 580 Mb on 1B ( Figure 5 and Table 2). It was significant for INC, SEV and IND, and explained as much as 20% of the phenotypic variation ( Table 2 and Supplementary Table 3). The QTL 1A.3 was located in the syntenic region of 1B.1, between 503 and 580 Mb ( Figure 5 and Table 2), and it was also significant for IND and INC though present in fewer environments and with lower significance than 1B.1 ( Table 2). 1B.1 colocated with significant HT and DTA QTL, while 1A.3 was significant for DTA.
Another major QTL was at 30-31 Mb on 2AS, termed 2A.1 (Figure 5 and   Table 2). It was detected at a low level for HT in one environment. 5B.2 was located from 577 to 691 Mb on 5BL ( Table 2). It explained up to 9.6% of phenotypic variation, and was most stable for INC in Brandon and Morden. This QTL was also associated with DTA, and minor effects were observed on IND and SEV, including at Indian Head (Figure 4 and Table 2). Three QTL were identified on chromosome 3B (Figure 5). The 3B.1 QTL was located around 3.7 Mb. It was identified in significant levels for INC, IND and SEV, and explained as much as 20.8% of phenotypic variation ( Table 2). The QTL also affected HT with a very large effect on DTA. A stable QTL, designated 3B.3, was located on chromosome 3B at 141-233 Mb ( Table 2). This QTL affected up to 9% of phenotypic variation, and also conferred a very stable effect on HT and smaller effect on DTA ( Table 2). The third 3B QTL, 3B.2, was located around 9.8 Mb ( Table 2), approximately 1 Mb from Fhb1 in common wheat (Rawat et al., 2016;Li et al., 2019;Su et al., 2019). It had no observable effect on HT or DTA, but also had a quite minor effect, explaining at most 7.7% of phenotypic variation ( Table 2). Though this QTL was less stable, because of the location of 3B.2 in the region of Fhb1 and because of the importance of this gene to FHB resistance in common wheat, we chose to further characterize the QTL in the durum AM panel. Pedigree information and genotypes of 3B.2 identified three different haplotypes for the significant marker, BS00079522_51, which were defined as tSumai3, tNative and tEmmer types (Supplementary Table 4). The tSumai 3 haplotype was derived from the introgression of Fhb1 from Sumai 3 into durum wheat (Supplementary Table 4). All Canadian cultivars shared the tNative haplotype, and the tEmmer haplotype was found in durum wheat introgressed from Td161 and a few durum wheat experimental lines from Austria (Supplementary Table 4). Allele effect analysis identified that the tEmmer type of 3B.2 conferred an effect that increased disease susceptibility (Figure 6).
There were a small number of QTL that did not co-locate with DTA or HT QTL. These include 1A.1, 1A.2, 6A.1, and 7A.3. The QTL 1A.1, located near the distal end of the short arm of chromosome 1A within a region from 13 to 20 Mb, was only significant for INC. Also on 1A was 1A.2, which mapped to 366 Mb on chromosome 1AL. It was detected in seven of the eight different environments, though not consistently across INC, IND and SEV, and a minor association with HT was also identified in one environment at this locus (Figure 4, Table 2, and Supplementary Table 3). 6A.1 was positioned at 12-23 Mb on 6AS. It had a significant effect on FHB, explaining up to 16% of phenotypic variation. This QTL also had a very minor effect for both HT and DTA with each only observed in a single environment. The 7A.3 QTL located to the distal region of 7A, around 671 Mb, had an effect on INC, SEV and IND, with no QTL for height or DTA found in this region (Figure 4 and Table 2). This QTL was detected only in Brandon and Morden field sites, and explained up to 9.6% phenotypic variation ( Table 2).

Phenotypic Data Analysis
The moderate to high heritability observed for FHB resistant traits in multiple environments in the durum AM panel indicated a large part of the phenotypic variation was contributed by genetic variation. The positive correlation between plant height and days to anthesis indicated that the genetic control of plant height and flowering time was partially shared ( Table 1; Langer et al., 2014). The high proportion of disease susceptibility we observed in the field tests supports literature emphasizing the limited tetraploid wheat resources with a high level of FHB resistance (Oliver et al., 2008). The observed significantly negative correlations between FHB resistance and plant height and days to anthesis also agreed with previous findings summarized by Prat et al. (2014) and Steiner et al. (2017). Because the significant negative correlations between both DTA and HT and FHB traits ranged from −0.24 * * * to −0.60 * * * and −0.18 * to −0.42 * * * , respectively, there is considerable scope to shift this negative relationships (i.e., to have DTA more consistently around −0.24 and the correlation with HT toward −0.18). By adopting strategies to stratify experimental genotypes into groups by both days to anthesis and plant height, it may be possible to recombine earlier to flower and shorter plants with reduced FHB symptoms. The correlation will not be broken but it can be shifted so that earlier maturing and shorter genotypes can be recombined with reduced FHB symptoms. Using this strategy, the negative relationship between plant height and FHB traits has been shifted by recombining semi-dwarf stature with a moderate level of resistance in hexaploid wheat cultivars such as Carberry (DePauw et al., 2011) and AAC Brandon (Cuthbert et al., 2017), both of which became widely adopted by producers. Adopting this strategy in durum wheat genetic enhancement could prove equally effective.

Genetic Architecture of FHB Resistance in the Durum AM Panel and Its Association With Flower Time and Plant Height
Compared to common wheat, durum wheat has limited genetic variation, and less effort has been committed to improve durum resistance to FHB (Buerstmayr et al., 2009(Buerstmayr et al., , 2019Prat et al., 2014Prat et al., , 2017. Within the current study, we identified a large number of QTL associated with FHB resistance with GWAS analysis from multiple environments and sources, broadening the resistance gene pool in durum wheat. The minor effect of these multiple QTL reinforces what is already known about the polygenic nature of FHB resistance, but also reveals the necessity of combining genes from multiple sources (Buerstmayr et al., 2009(Buerstmayr et al., , 2019Liu et al., 2009).
The major and most consistent FHB QTL found in previous studies is the hexaploid wheat Sumai 3 derived Fhb1, located on 3BS around 7.6-13.9 Mb Liu et al., 2006).   Introgression of Fhb1 into durum wheat has been challenging, with one possible reason being the unstable expression in a durum genetic background (Zhao et al., 2018). Recently, Prat et al. (2017) successfully introgressed Fhb1 into durum wheat, and some of those introgression lines are part of this AM panel. A QTL was found in the same region as Fhb1 in this study, designated 3B.2. This QTL was detected in limited environments with a minor effect. QTL 3B.2 had three distinct haplotypes (Supplementary Table 4), and compared to haplotypes of Sumai 3 (tSumai 3) and Canadian cultivars (tNative), the haplotype from the experimental lines derived from emmer wheat Td161 (tEmmer) conferred disease susceptibility (Supplementary Table 4). This finding confirms previous findings that the Fhb1 region from Td161 contributed to disease susceptibility when compared to the susceptible durum wheat Floradur (Buerstmayr et al., 2012). The resistance haplotype found in the GWAS study by Steiner et al. (2019b) corresponds to the tNative haplotype presented in this study. The tNative haplotype is the only haplotype found in the Canadian and American cultivars presented in both studies, while both the tNative and tEmmer haplotypes exist in durum wheat from Austria, CIMMYT, ICARDA, Italy and Morocco (Steiner et al., 2019b). Altogether, these findings indicate that one of the two non-Sumai 3 Fhb1 region haplotypes found in tetraploid wheat contributed to disease susceptibility when compared to the other. Further characterizing the region with additional markers is needed to help resolve the source of the alleles and further understand the effects of the three haplotypes identified in this study. Two additional 3B QTL were found significantly associated with all of the traits, 3B.1 in the telomeric region of 3BS, and 3B.3 in the centromeric region of the short arm (3BSc). Recently, Wu et al. (2019) reported a QTL positioned at 2.0 Mb on the reference sequence from elite Chinese common wheat germplasm, almost the same region as the 3B.1 identified in this durum AM panel study. The 3B.3 QTL was one of the most stable QTL identified, with a larger effect on FHB resistance than other QTL in this AM panel. Notably, the resistant 3BSc haplotypes were identified in the durum wheat lines that also had Fhb1 introgressed from Sumai 3 by Prat et al. (2017). The location of 3B.3 corresponds to the 3BSc region QTL previously reported as important to FHB resistance, particularly in Canadian elite germplasm, where 3BSc conferred a larger effect than Fhb1 (McCartney et al., 2007). Also in agreement with findings from McCartney et al. (2007), the 3BSc QTL conferred a large effect on both plant height and DTA in elite Canadian wheat. Further research is needed to explore effects of Fhb1, 3B.1 and 3BSc in durum wheat.

QTL With No or Weak Association With Flowering Time and Height
The common association between plant height, flowering time and FHB resistance was illustrated in this study. Of the 31 FHB QTL regions identified, all but five also had strong associations with plant height and/or flowering time. The relatively small effects of these QTL compared to other QTLs detected in this study may be related to the strong influence of flowering time on FHB resistance, potentially overinflating the effects of the QTL for FHB resistance due to the timing of flowering. Due to the progression of the FHB symptoms over time, the correlation between days to anthesis and disease development are confounded by the length of time for disease development. Due to cost constraints, disease rating was not evaluated over a time course to control for this effect, and thus we cannot exclude the observed correlation between FHB resistance and DTA may be caused by these confounding effects.
Fusarium head blight resistance QTL that are not associated with height or flowering time are much more appealing targets, as the negative influence of taller plants and complicated relationship with flowering time can be avoided. The targeted breeding of these QTL for resistance that do not carry extra undesirable traits will have the most likely success. The most favorable of these QTL may be 3B.2, but the QTL 1A.1, 1A.2, 6A.1, and 7A.3 with no association or weak association with DTA and HT are also desirable candidates. The 1A.1 QTL was located in the same region as the major QTL previously reported on the distal part of 1AS (summarized by Buerstmayr et al., 2009;Liu et al., 2009;Venske et al., 2019). Jiang et al. (2007a,b) located an FHB SEV QTL from the Chinese wheat line CJ9306 to position 27.2 Mb, and GWAS by Zhu et al. (2020) similarly identified an FHB QTL for IND from Chinese elite germplasm in the same region. A recent study by Sari et al. (2018) (Wu et al., 2019) and for FHB resistance based on point inoculation in CIMMYT line C615 (Yi et al., 2018). In our study, we found this QTL was also associated with FHB incidence, index and severity. Within the AM panel of our study, although the resistance allele of 1A.1 was not found in Canadian cultivars, the 1A.2 occurred in several current Canadian cultivars with improved FHB resistance, including CDC Precision (Pozniak and Clarke, 2017b) and Brigade (Clarke et al., 2009;Supplementary Table 3).
The 6A.1 QTL's large effect on FHB resistance makes it appealing despite a small undesirable influence on DTA and HT. No major QTL clusters have been reported in a similar region as 6A.1, though Yi et al. (2018) reported a minor QTL in this region detected from a susceptible wheat line in one environment, and Lu et al. (2013) identified a minor QTL in the proximal 6A region for both FHB resistance and plant height. Because the 6A.1 resistance haplotype is present in a large number of Canadian durum wheat cultivars, including Brigade (Clarke et al., 2009), Transcend (Singh et al., 2012), CDC Credence (Sari et al., 2018) and CDC Precision (Pozniak and Clarke, 2017b;Supplementary Table 3), it should be possible for Canadian breeding programs to build on this resistance, though the effect of the QTL in Canadian elite durum cultivars remains to be validated. The 7A.3 QTL, located at 671 Mb, with its relatively large effects on all FHB resistant traits without being associated with plant height or flowering time also make it another good target for breeding FHB resistance. Previous research identified a major QTL for type II resistance based on point inoculation in the vicinity of 7A.3 through the physical mapping of the SSRs gwm276 and gwm262 to positions of 642.9 and 681.4 Mb (Semagn et al., 2007;Buerstmayr et al., 2009). Wu et al. (2019) also reported a QTL affecting DON accumulation in the same region of elite Chinese germplasm, while Sari et al. (2018) reported QTL for SEV and IND in the same region from the durum wheat inbred line DT696.
From the durum AM panel in our study, 2A.2 located in the same region as a native durum FHB resistance QTL in previous research in cultivars Ben by Zhang et al. (2014) and Joppa by Zhao et al. (2018). In addition, the QTL 2A.2 was also found consistently associated with DTA, suggesting it plays a role in controlling flowering. In this durum AM panel, the resistance haplotype of 2A.2 was found in DT696 (Sari et al., 2018), an adapted source of FHB resistance in durum wheat, as well as several Canadian cultivars with improved FHB resistance derived from this line, including Brigade (Clarke et al., 2009), Transcend (Singh et al., 2012) CDC Credence (Sari et al., 2018), and CDC Precision (Pozniak and Clarke, 2017b;Supplementary Table 3). Despite its association with DTA, the effectiveness of the 2A.2 in native durum cultivars from Canada and United States make it another good target to breed durum wheat with improved FHB resistance.

QTL Co-located With Flowering Genes
The majority of the QTL identified from this AM panel were found associated with flowering time and/or plant height. As mentioned previously, the Notably, three QTL pairs, including 1A.3 and 1B.1, 2A.1 and 2B.1, and 5A.1 and 5B.2, were found in syntenic regions of the A/B genome that harbor known orthologous gene pairs controlling flower time. 1A.3 was in a similar region of a major QTL found in United States winter wheat cultivar NC-Neuse (Petersen et al., 2016(Petersen et al., , 2017. The FLOWERING LOCUS T3-A1 (TaFT3-A1) gene that promotes flowering was found physically mapped around 528.1 Mb of 1A in CS Ref 1.0 (Zikhali et al., 2017;International Wheat Genome Sequencing Consortium [IWGSC], 2018), which is close to the region of 1A.3. The major QTL 1B.1 located to the region coinciding with a QTL of FHB resistance from the European winter wheat Arina (Semagn et al., 2007;Buerstmayr et al., 2009;Liu et al., 2009), as well as loci controlling DTA identified in the recent durum wheat GWAS by Steiner et al. (2019b). This QTL conferred a stable and large effect for INC, SEV, HT and DTA. Recently, the photoperiod gene FLOWERING LOCUS T3-B1 (TaFT3-B1) that promotes flowering time, was physically identified at position 581 Mb of 1B (Zikhali et al., 2017), the same region as 1B.1. The 1B.1 and 1A.3 QTL occur in syntenic region of the genome, indicating the orthologous gene pair, TaFT3-B1 and TaFT3-A1, as candidate genes underling the QTL effect in these regions.
The 2A.1 QTL conferred main effects for INC, IND, DTA and HT, physically positioned to around 27-31 Mb on chromosome 2A. This location is very near to the photoperiod gene Ppd1A, which has an important role in controlling flowering time and height, indicating 2A.1 as candidate gene controlling the QTL. Giancaspro et al. (2016) found a similar QTL positioned at 10 Mb on 2AS for FHB resistance in durum wheat, derived from the introgression of FHB resistance from Sumai 3, but with no report on its association with plant height. Gadaleta et al. (2019) identified a wall-associated receptor-like kinase (WAK2) in this region as the candidate gene for FHB resistance. Our study found a 2B QTL, designated 2B.1 that colocalizes with Ppd1B located in a syntenic region of 2A.1. This QTL contributed to INC, SEV, IND, DTA and HT. Thus, our findings support the Ppd loci on 2AS and 2BS as candidate genes responsible for the observed effects, although further studies with well stratified plant height and FHB rating DTA are required in order to explore the factors underlying these QTL.
Both the QTL 5A.1 on 5AL and 5B.2 on 5BL occur in syntenic regions that harbor orthologs of the well-known vernalization genes VRNA1 (at 585.1 Mb) and VRNB1 (at 613.0 Mb). 5A.1 and 5B.2 both conferred a stable effect for INC and IND, and while 5B.2 also had a large effect of on DTA, 5A.1 had no effect on DTA and only a minor effect on HT in one environment. Sari et al. (2018) reported a major FHB resistance QTL from the Canadian durum wheat line DT696 in the same region as 5A.1, also finding no DTA or HT QTL in this region. Xu et al. (2020) found QTL located in the same regions as 5A.1 and 5B.2 in common wheat that controlled anther extrusion, heading time and FHB resistance. There is potential that these vernalization genes are responsible for the FHB resistance coming from these regions, and that the VRNA1 gene has just a minor effect on flowering time in durum wheat. The resistance haplotype of 5A.1 was found in Canadian durum cultivars including Brigade (Clarke et al., 2009) and CDC Alloy (Pozniak and Clarke, 2017a;Supplementary Table 3). Because of the presence of the resistant haplotype in current durum cultivars, and the minor effect on flowering time, we believe the VRNA1 region QTL from this study and Sari et al. (2018) is a good target to improve FHB resistance in durum wheat. However, there is still need for further research to explore the mechanism of colocalization between the vernalization genes and FHB resistance and their effect on flowering in durum.

CONCLUSION
With genome wide association analysis we identified 31 QTL for FHB resistance. This confirms the quantitative nature and polygenic control of the FHB resistance and also signifies that this durum AM panel contains a large amount of genetic variation for FHB resistance loci. These QTL capture a large amount of the major QTL reported for hexaploid and tetraploid wheat which should facilitate improving FHB resistance in durum wheat. Five QTL found primarily for FHB resistance, including 1A.1, 1A.2, 5A.1, 6A.1, and 7A.3, could be used as initial targets to improve resistance in durum wheat without detrimental effects. Although 2A.2 is associated with DTA, the resistant haplotype exists in several Canadian and United States cultivars with improved FHB resistance, and we think that due to its adaption to durum cultivars in North America it is also a good target. The majority of these QTL identified were associated with plant height and/or flowering time, indicating that phenology, flowering and height genes formed a complex network affecting FHB resistance in durum wheat. Prior knowledge of the haplotypes of these genes in breeding materials will provide an informed approach to stack these genes and give breeders the ability to design a better strategy to use these sources to improve FHB resistance. However, more research is needed to identify the mechanism of the trait associations, and truly determine whether pleiotropic effects of same gene, linkage drag of resistant genes, and/or disease escape due to flowering time and plant height are in effect. Only by completely understanding these relationships, can a better strategy, from genetic, genomics and breeding perspectives be developed to significantly increase FHB resistance in durum wheat. Finally, considering the attributes of QTL identified in this study, including the large number of minor effects, the varied expression across environments, and the complex interaction with flowering time and height, we suggest intercrossing the multiple sources of resistance. Then the progeny should be selected using a multi-trait based, high-throughput marker assisted selection approach that incorporates resistance, flowering time and height loci, in combination with intensive phenotyping, with the genotypes grouped by days to flower and plant height, across multiple target environments, as the most promising approach to develop durum wheat with a better level of resistance.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
YR, RD, RK, PF, and WZ conceived and designed this study. YR and RK constructed the durum association mapping (AM) population. HB, RD, RK, RC, and YR contributed to the durum germplasm development in the AM population. YR, RK, HC, and SB contributed to seed increase of this AM population and the design of field trials. MH, AB, and SK conducted field trials and disease evaluation at FHB nurseries located in Morden and Brandon, Manitoba. YR, SB, RK, RC, and HC performed field trials and disease evaluation at FHB nursery in Indian Head, Saskatchewan. KB and BP contributed to the 90K SNP genotyping and marker identification. KB, WZ, and YR analyzed the data and interpreted results. KB, WZ, RR, and YR contributed to data validation. WZ, KB, and YR wrote the original draft. KB, RK, WZ, YR, RD, HB, HC, and PF reviewed and edited the manuscript. YR was the principal investigator and supervised the project. All authors contributed to the article and approved the submitted version.