Genetic structure and genome-wide association study of a genomic panel of two-row, spring barley (Hordeum vulgare L.) with differential reaction to Fusarium head blight (Fusarium graminearum Schwabe) and deoxynivalenol production

Abstract

Fusarium head blight (FHB), primarily incited by Fusarium graminearum Schwabe, is the most devastating disease of barley (Hordeum vulgare L.) in Canada. Contaminated grains are unsuitable for use as livestock feed or by the malting industry, primarily due to contamination by mycotoxins such as deoxynivalenol (DON). Reducing DON content has been a long-term goal of barley breeders; however, the complex genetics of resistance and laborious testing required have made it difficult to develop resistant cultivars. An Illumina 50 K SNP beadchip was used to genotype a diverse collection (n = 400) of two-row barley genotypes, selected primarily from North American-breeding programs. This genomic panel was phenotyped in three environments in Manitoba, Canada, over two growing seasons. Genotypes were evaluated for heading date, height, FHB score, and DON content. A genome-wide association study (GWAS) was carried out using a mixed linear model (MLM: Q + K) accounting for population structure (Q) and kinship (K) as covariates. This study characterized the population structure of two-row barley germplasm important to North American-breeding programs. Significant marker-trait associations were identified on all chromosomes for FHB, and on all chromosomesexcept 1 H and 6 H for DON content. Individual marker effects were small, explaining only up to 5% of the phenotypic variation. Genes within genomic regions that were associated with FHB-related traits displayed multiple functions linked with disease resistance. While numerous minor, marker-trait associations were identified, phenotypic evaluations are still the best predictors of these traits, despite their time-consuming nature and dependenceon environmental conditions.

Résumé

La brûlure de l’épi causée par le fusarium (BEF), provoquée principalement par Fusarium graminearum Schwabe, est la maladie de l’orge (Hordeum vulgare L.) la plus dévastatrice au Canada. Les grains contaminés ne peuvent être utilisés pour l’alimentation du bétail ni par l’industrie brassicole, notamment à cause de la contamination par des mycotoxines telles que le désoxynivalénol (DON). La réduction du contenu en DON a été un objectif à long terme des sélectionneurs d’orge. Toutefois, la génétique complexe de la résistance et les essais laborieux requis ont rendu difficile le développement de cultivars résistants. Une puce Illumina 50 K SNP a été utilisée pour génotyper une collection variée (n = 400) de génotypes d’orge à deux rangs, issus principalement des programmes nord-américains de sélection. Ce panel génomique a été phénotypé dans trois environnements manitobains, au Canada, durant deux saisons de croissance. Les génotypes ont été évalués en fonction de la date d’épiaison, de la hauteur, du taux de BEF et du contenu en DON. Une étude d’association pangénomique (EAP) a été menée avec un modèle mixte linéaire (MML: Q + K) qui prenait en compte la structure de la population (Q) et la parenté (K) comme covariables. Cette étude a caractérisé la structure de la population du germoplasme de l’orge à deux rangs, essentiel pour les programmes de sélection nord-américains. Des associations importantes de marqueurs et de caractères ont été détectées sur tous les chromosomes pour la BEF, et sur tous les chromosomes sauf 1 H et 6 H pour le contenu en DON. Les effets des marqueurs individuels étaient minimes, expliquant seulement jusqu’à hauteur de 5 % la variation phénotypique. Des gènes des régions génomiques qui étaient associés aux caractères de la BEF ont affiché de nombreuses fonctions liées à la résistance à la maladie. Bien que de nombreuses associations mineures de marqueurs et de caractères aient été répertoriées, les évaluations phénotypiques demeurent toujours les meilleurs indicateurs de ces caractères, et ce, malgré leur nature chronovore ainsi que leur dépendance aux conditions environnementales.

Keywords:

• barley
• deoxynivalenol
• disease resistance
• Fusarium head blight
• GWAS
• linkage mapping
• quantitative

Mots clés:

• Orge, désoxynivalénol
• résistance à la maladie
• brûlure de l’épi causée par le fusarium
• EAP
• cartographie de liaison
• quantitatif

Introduction

Fusarium head blight (FHB) is a devastating disease of barley (Hordeum vulgare L.) and other cereals, caused primarily by Fusarium graminearum Schwabe (teleomorph: Gibberella zeae (Schwein.) Petch), which can significantly reduce yield and quality of grain. Fungal infections are characterized by stained and/or shriveled kernels which are sometimes covered with pink mycelium. Of particular importance are B-type trichothecene mycotoxins such as deoxynivalenol (DON) which are commonly associated with this disease (Laitila 2015). Contaminated barley grains are not acceptable for use as livestock feed or malt production. Deoxynivalenol is thermostable and it and its conjugated form (deoxynivalenol-3-glucoside, D3G) can be transmitted to beer made from infected barley (Lancova et al. 2008). Due to food safety issues, strict maximum limits have been set for DON by the malting industry for commercial use (< 0.5 mg kg⁻¹) (He et al. 2015). Since the mid-1990s, FHB has become the most damaging threat to barley production in western Canada (Clear et al. 1996; Tekauz et al. 2000).

The development of genetically resistant cultivars remains the most sustainable approach to FHB management. However, resistance sources are limited in barley (Bai and Shaner 2004) and provide only partial resistance. Fusarium head blight resistance is complex and under polygenic control. Numerous bi-parental quantitative trait loci (QTL) mapping studies have been conducted with various FHBresistance sources (de la Peña et al. 1999; Zhu et al. 1999; Ma et al. 2000; Dahleen et al. 2003; Mesfin et al. 2003; Hori et al. 2005, 2006; Horsley et al. 2006; Yu et al. 2010). The FHB-related QTL identified were commonly associated with plant height and heading date, and were often environmentally specific. Near-isogenic lines (NILs) and transcriptomics have been used to study functions of these FHB-related QTL in six-row barley located on chromosome 2 H Bins 8 & 10; 3 H Bin 6; 6 H Bin 7 (Jia et al. 2011; Huang et al. 2016). Linkmeyer et al. (2013) applied marker-assisted selection (MAS) to incorporate a resistant allele for Qrgz-2 H-8, but concluded that the introgressions did not attain resistance as high as observed for the Swiss-landrace ‘Chevron’, which is one of the best sources of partial resistance. Therefore, MAS has had limited application as a breeding method for developing FHB resistance in barley.

Modern technologies based on single nucleotide polymorphism (SNP) markers have emerged to facilitate the application of whole-genome scans in barley (Waugh et al. 2009). The Illumina 9 K iSelect SNP chip (7842 markers) consists of 2832 original Illumina GoldenGate BOPA SNPs (Close et al. 2009) plus over an additional 5 K of Next-Generation sequence-based SNPs (Comadran et al. 2012), whichwere applied in QTL-mapping of a recombinant inbred line (RIL) population to identify genomic regions involved in FHB resistance. Huang et al. (2018) mapped FHB and DON QTL in a RIL population involving ‘Rasmusson’, a six-row barley cultivar with ‘Chevron’ (Moderately resistant – MR) in ancestry. Several QTL were identified on multiple chromosomes for FHB (2 H, 5 H, 6 H, and 7 H) and DON content (2 H, 3 H, and 7 H). Ogrodowicz et al. (2020) also used the 9 K SNP chip for QTL mapping in a two-row RIL population and identified regions associated with FHB index on several chromosomes (2 H, 3 H, 5 H, and 7 H), but found no association for DON content. However, several of the QTL associated with resistance in these studies were also negatively associated with important agronomic traits.

Genome-wide association study (GWAS) is a powerful tool that uses whole-genome markers, and offers an alternative approach to QTL mapping which is restricted to bi-parental populations. Genome-wide association study has the benefits of higher resolution by relying on more-extensive, historical recombination events (Bush and Moore 2012). Its power,however, is dependent on sample size and the number of molecular markers used. Using the Illumina GoldenGate or the iSelect 9 K SNP chip, GWAS has identified numerous QTL for FHB and DON within six-row contemporary breeding germplasm from American barley-breeding programs (Massman et al. 2011; Sallam et al. 2015). Mamo and Steffenson (2015) applied GWAS within a mixed panel of 2-row/6-row landraces from Ethiopia and Eritrea, where they identified associations in the 2-row sub-group for FHB (1 H, 4 H, and 7 H) and DON content (2 H, 4 H, and 5 H). Bedawy et al. (2018) used a combination of 9 K SNP markers and diversity array (DArT) markers to identify QTL associated with a F. graminearum leaf assay and spike inoculations in a greenhouse (1 H and 5 H).The majority of past research contributions in genetics have been focused on six-row resistance development, while two-row barley is by far the dominant type produced in Canada.

Over the past two decades, FHB resistance has been enhanced through an accumulation of minor alleles in lines tested in large screening nurseries (Legge et al. 2004). Progress has led to the release of moderately resistant cultivars with lowered DON content, developed through exotic introgressions, e.g., ‘AAC Connect’ and ‘AAC Goldman’ (Legge et al. 2017, 2018), or identification of transgressive-segregants,e.g., ‘CDC Mindon’ (Rossnagel et al. 2008) and ‘Lowe’ (Juskiw et al. 2019).In this study, a large genomic panel of contemporary two-row barley (n = 400) along with the more recent, larger barley 50 K SNP microarray chip (containing 44 040 molecular marker assays) (Bayer et al. 2017) was used to assess the population structure of the North American assembled germplasm collection and perform GWAS to identify genomic regions associated with quantitative resistance to FHB and low DON accumulation under western Canadian field conditions. This study provides valuable information that could support the barley breeders in their quest to develop resistance to this destructive disease.

Materials and methods

The plant materials used in this study consisted of 400 two-row, hulled barley genotypes with spring-habit, primarily adapted to production within North America. Germplasm was selected based on the previous reaction to FHB and DON content, as tested primarily in a FHB nursery at Brandon, MB. Genotypes originating from elite by elite-breeding crosses and/or exotic introgressions were incorporated. Pedigrees were taken into consideration to increase diversity and exclude half-sibs.

The majority of the panel (> 75%) consisted of cultivars and advanced breeding lines from the three western Canadian barley breeding centres: Brandon Research and Development Centre, Agriculture and Agri-Food Canada, Brandon, MB (AAFC-BRDC); Crop Development Centre, University of Saskatchewan, Saskatoon, SK (CDC); Field Crop Development Centre, Olds College, Lacombe, AB (FCDC) (Fig. 1). In addition to the material from these breeding institutions, genotypes were also selected from breeding programs belonging to other organizations such as Busch Agricultural Resources Inc. (BARI); WestBred/Highland Specialty Grains (WB/ HSG); North Dakota State University (NDSU), Ottawa Research and Development Centre, Agriculture and Agri-Food Canada (AAFC-ORDC), Lethbridge Research and Development Centre, Agriculture and Agri-Food Canada (AAFC-LeRDC), Molson Coors Brewing Co (COORS) and Agricore United (AU). Additionally, a small portion (< 5%) of the panel included international cultivars or landrace accessions. Furthermore, the panel consisted of genotypes that could be classified as malting (n = 171) and general-purpose (n = 229, feed/forage) barley. For the majority of the barley genotypes, all of the seed used in the study was received directly from the breeders except for the landraces which were obtainedfrom a gene bank (Agriculture and Agri-Food Canada, Plant Gene Resources of Canada).

Fig. 1. Pie chart of genetic panel composition by the breeding institutions (Agriculture and Agri-Food Canada, Brandon Research and Development Centre [AAFC-BRDC]; Crop Development Centre [CDC]; Field Crop Development Centre [FCDC]; D = Busch Agricultural Resources, Inc. [Bari]; WestBred LLC/ Highland Specialty Grains [WB/HSG]; North Dakota State University [NDSU]; Agriculture and Agri-Food Canada, Lethbridge Research and Development Centre [AAFC-LeRDC]; Agriculture and Agri-Food Canada, Ottawa Research and Development Centre [AAFC-ORDC]; Molson Coors Beverage Company [COORS]; Agricore United [AU]; Accessions are landrace or exotic cultivars).

Phenotype

Field experiments. In 2014 and 2015, the experiments were conducted in field FHB nurseries at three locations in Manitoba, Canada: Brandon (49°51ʹ56.1”N 99°58ʹ57.9”W), Carman (49°30ʹ05.1”N 98°01ʹ43.5”W) and Carberry (49°52ʹ24.8”N 99°24ʹ30.8”W). Standard agronomic practices were applied within locations, with the following seeding dates: Brandon (May 28, 2014; May 20, 2015), Carberry (May 19, 2014; May 20, 2015), Carman (May 29, 2014; June 1, 2015). All 400 two-row barley genotypes of the panel were grown together in a randomized complete block design (RCBD), at a seeding rate of 30–40 seeds per meter. All sites were irrigated with various methods: Brandon (PVC risers with NaanDanJain 502-H sprinkler head); Carberry (center-pivot overhead irrigation); Carman (mist irrigation). While inoculation methods differed by site, all nurseries used mixtures of two isolates each of 3ADON and 15ADON chemotypes of F. graminearum. Before inoculation, isolates were mixed in equal proportions by volume. Grain spawn (maize kernels colonized by F. graminearum) was spread on the soil surface (5 g m⁻²) at Brandon and Carberry, whereas Carman plots were sprayed at 75% spike emergence with a macroconidial suspension (5 x 10⁴ spores ml⁻¹) and then at three day-intervals thereafter. Following inoculation at Carman, mist was applied at ten-minute intervals for ten hours each time to ensure successful infection.

Plots were rated for FHB3.5 weeks following heading, during the soft dough growth stage (Zadoks – Z85, Zadoks et al. 1974) using a 0–5 visual scale (A. Tekauz personal communication, 0 = no infection; 1 = incidence low, up to 5% of spikes infected; severity low, 1 or 2 kernels per spike affected [up to 7% of spike]; 2 = incidence low to moderate, 5 to 15% of spikes infected; severity low to moderate, 1 to 4 kernels per spike affected [up to 15% of spike]; 3 = incidence moderate, 15 to 30% of spikes infected; severity moderate, 2 to 8 kernels per spike affected [up to 25% of spike]; 4 = incidence moderate to high, 30 to 50% of spikes infected; severity moderate to high, 4 to 12 kernels per spike affected [up to 40% of spike] and 5 = incidence high, 50% or more spikes infected; severity high, 5 to 15+ kernels per spike affected [up to 50%+ of spike diseased]). Heading (50% spike emergence) (HD50) was recorded only at Brandon and Carberry and used in the calculation of days to heading by subtracting the seeding date. Height (HT) was measured as the distance from the soil surface to the top of the spike excluding awns and was recorded as an average of three measurements at all three locations.

Grains were harvested at maturity with a stationary, research combine set at low wind speed to retain thin kernels. Grains were cleaned using an SLN3 sample cleaner (Pfeuffer GmbH, Kitzingen, Germany). A 20 g subsample that was free of debris and/or chaff was then ground using a Perten 3610 lab mill with a fine particle disc set (Perten Elmer Inc. Waltham, MA). Deoxynivalenol content was analyzed on a 2-gram sub-sample by enzyme-linked immunosorbent assay (ELISA) technique using Veratox ® 5/5 (Neogen Corporation, Lansing, MI) as per the manufacturers protocol (limit of detection = 0.1 mg kg⁻¹). Sample extracts were analyzed in pairs, where analysis was repeated for deviations above 10%.

The least square (LS) means of all traits were calculated for genotypes using SAS JMP Genomics 9.1 software (SAS Institute Inc., 1989–2019) in a model overall years with variables as fixed effects. LS means were used for descriptive statistics and inputs for modelling genome-wide association study (GWAS). Pearson correlation coefficients were calculated in R statistic environment (R Core Team 2021) using the ‘cor’ function and heat map graphics were produced with the ‘ggcorrplot’ package.

DNA extraction

Two barley seeds were surface-sterilized in 6% vol/vol sodium hypochlorite and rinsed in distilled water. Seeds were germinated on moist cotton balls and grown under natural light until the two-leaf stage. Leaves were randomly cut from a single plant and flash frozen in liquid-N then freeze-dried in a lyophilizer (Labconco Corporation, Kansas City, MO). Total genomic DNA was extracted from 100 mg of lyophilized tissue per genotype using Qiagen, DNeasy 96 Plant Kit (Qiagen, Canada). The quality and concentration of DNA were evaluated by agarose gels and a Nanodrop^TM 1000 spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE). Samples were normalized to 50 ng µL⁻¹ (20 µL) of DNA using micro-pure DNase-free water.

High throughput genotyping

DNA samples were assayed on an Illumina iScan (Illumina, San Diego, CA) using a custom iSelect – 50 K SNP microarray chip consisting of 44 040 molecular markers (Bayer et al. 2017) at AAFC, Morden Research and Development Centre, Morden, MB. A custom cluster file (Trait Genetics, GmbH, Germany) was used to call SNP alleles using Illumina Genome Studio v2.0.5 software. Allelic data in Plink format were imported into TASSEL v5.2.58 then filtered for minor allele frequency (MAF) > 5% and for missing data <10% (total of 25 068 SNPs). Missing values were numerically imputed by the mean value method using five nearest neighbours. Data were imported into SAS JMP Genomics 9.1 (SAS Institute Inc., 1989–2019) and converted into numerical values (0 = major allele; 1 = heterozygous; 2 = minor allele).

Linkage disequilibrium and population structure analysis

Resolution of trait mapping is affected by linkage disequilibrium (LD). LD decay was analyzed separately for each chromosome. Squared allele frequency correlations (r²) were calculated pairwise for all filtered molecular markers using SAS JMP Genomics version 9.1 (SAS Institute Inc 2019). Intra-chromosomal graphs were produced by plotting r² values of marker pairs against physical genomic map position (base pair) over each chromosome with a smoothing line function (α = 0.005).

STRUCTURE v2.3.4 software (Pritchard et al. 2000) was used to perform a Bayesian clustering analysis of population structure. A sub-set of > 2000 polymorphic SNP markers of equidistant spacing overall seven chromosomes was used with a burn-in of 100 000 iterations and 100 000 Markov Chain Monte Carlo (MCMC) iterations. Population sub-populations from K1 (no structure) to K5 were tested using an admixture model. Results were analyzed using Structure Harvester v0.6.94 (Earl and vonHoldt 2012) where mostlikely sub-populations were identified via the Evanno method (Evanno et al. 2005). A ΔK index was calculated and plotted against the K value to observe and identify the peak value and determine the mostlikely number of sub-populations. A Q matrix (sub-population coefficients) was extracted from the population structure analysis based on two sub-populations. Visual analysis of the Q matrix output by STRUCTURE was conducted using Structure Plot v2.0 (Ramasamy et al. 2014). Alternatively, population structure was evaluated by TASSEL v5.2.58 through principal component analysis (PCA) along with the production of a neighbour-jointing tree visualized by Archaeopteryx v0.9909 to evaluate patterns of relatedness. A matrix of pairwise kinship values (K matrix) was also calculated within SAS JMP Genomics version 9.1 (SAS Institute Inc 2019) via anidentity-by-descent method. The K-matrix was used in subsequent GWAS analysis as a random factor.

Association mapping was conducted using all 400 genotypes for traits via various methods (with Q-matrix based on PCA or STRUCTURE results; general linear model (GLM + Q) and mixed linear model (MLM + Q + K) were evaluated by quantile-quantile (Q-Q) plots for best fit (Supplementary Figure S1). For simplicity, results herein are attributed to the Q (STRUCTURE) + K mixed model. Association analysis was conducted on overall means and within a location. SNP markers -Log₁₀ [P-value] ≥ 3 were considered statistically significant, a commonly accepted limit for the GWAS study (Chan et al. 2010). Manhattan plots were generated in SAS JMP Genomics 9.1 (SAS Institute Inc 2019) to visually identify peaks of trait associations by physical chromosome position.

Candidate gene prediction

Genes were determined by the physical position of the SNP allele on the Morex V3(MorexV3_pseudomolecules_assembly, accession GCA_904849725) which is the most recent released version of the barley reference genome (Mascher et al. 2021).Associated genes of interest were identified through evaluation of nearest proximity to genomic position (± 50 Kbp) and often adjacent to the SNP marker using the Barleymap web tool (URL floresta.eead.csic.es/barleymapaccessed Apr. 12, 2022)and the GrainGenes genomic web browser (URL wheat.pw.usda.gov/GG3/genome_browser accessed on Apr. 12, 2022. High-confidence gene annotations were determined through BARLEX genome explorer (URL apex.ipk-gatersleben.de/apex/f? p = 284:10 accessed Apr.12, 2022). Gene functions were further evaluated by associated gene ontology (GO) terms (Ashburner et al. 2000, http://geneontology.org/); Pfam domains (Mistry et al. 2020, URL pfam.xfam.org/ accessed Apr. 12, 2022); InterPro terms (Blum et al. 2021, URL www.ebi.ac.uk/interpro/ accessed Apr. 12, 2022).

Results

Population structure analysis

STRUCTURE software was applied to evaluate sub-structures within the panel of 400 genotypes. The optimal number of sub-populations (K) determined by STRUCTURE analysis, identified K = 2 sub-populations (Supplementary Figure S2). Individuals from Cluster 1 (> 90% sub-population membership) included 104 lines primarily contributed by the AAFC-BRDC-breeding program, or other programs where an AAFC-BRDC parent was involved in a cross. Representative individuals of Cluster 1 included cultivars and elite breeding materials such as: ‘AAC Connect’, ‘AAC Goldman’, ‘AAC Synergy’, ‘AC Oxbow’, ‘CDC Fraser’, ‘Cerveza’, ‘CDC Kindersley’, ‘Ellice’, ‘Major’ ‘Newdale’. Cluster 2, comprised all remaining entries, included cultivars and lines from various programs and accessions (Fig. 2). Genotypes in extreme opposition to Cluster 1, included feed cultivars such as ‘Brahma’, ‘McLeod’, ‘Ursa’, and ‘Xena’. This structure was also reflected in the phylogenetic tree constructed by the neighbor-joining method (Fig. 3).

Fig. 2. STRUCTURE Sub-population membership for 400 two-row barley grouped by breeding institution. A = Agriculture and Agri-Food Canada, Brandon Research and Development Centre (AAFC-BRDC); B = Crop Development Centre (CDC); C = Field Crop Development Centre (FCDC); D = Busch Agricultural Resources, Inc. (Bari); E = WestBred LLC/ Highland Specialty Grains (WB/HSG); F = North Dakota State University (NDSU); G = Agriculture and Agri-Food Canada, Lethbridge Research and Development Centre (AAFC-LeRDC); H = Agriculture and Agri-Food Canada, Ottawa Research and Development Centre (AAFC-ORDC); I = Molson Coors Beverage Company (COORS); J = Agricore United (AU); K = Accession. Yellow = Cluster 1 and Orange = Cluster 2. The Y axis represents the proportion of cluster (K = 2) membership.

Fig. 3. Phylogenetic tree constructed by neighbor-joining (NJ) method based on single-nucleotide polymorphism (SNP) genetic markers (25 K) displaying selected cultivars and germplasm from the various institutions (Agriculture and Agri-Food Canada, Brandon Research and Development Centre [red]; Crop Development Centre [blue]; Field Crop Development Centre [green]; D = Busch Agricultural Resources, Inc. [purple]; WestBred LLC/ Highland Specialty Grains [Orange]; North Dakota State University [brown].

Principal component analysis (PCA) was also conducted to analyze the panel structure. Analysis of principal components demonstrated clustering within barley type (malt or feed) and separation on the primary principal component axis. Results of the PCA are presented by type (Fig. 4A) and by breeding program origin (Fig. 4B). Patterns of genetic relatedness from various analyses were confirmed by pedigree information (breeders’ personal communication).

Fig. 4. Plots of the first two principal components of genotype information (single-nucleotide polymorphism genetic markers) of the genomic panel with colouration highlighting A) malting or feed barley types; B) breeding institution origin (Agriculture and Agri-Food Canada, Brandon Research and Development Centre [AAFC-BRDC]; Crop Development Centre [CDC]; Field Crop Development Centre [FCDC]; D = Busch Agricultural Resources, Inc. [Bari]; WestBred LLC/ Highland Specialty Grains [WB/HSG]; North Dakota State University [NDSU]; Agriculture and Agri-Food Canada, Lethbridge Research and Development Centre [AAFC-LeRDC]; Agriculture and Agri-Food Canada, Ottawa Research and Development Centre [AAFC-ORDC]; Molson Coors Beverage Company [COORS]; Agricore United [AU]; Accessions are landrace or exotic cultivars).

Linkage disequilibrium and polymorphic SNP coverage

The genomic panel was primarily composed of adapted breeding materials, where it is expected to show some level of LD decay as a result of selective breeding. The degree of population structuring is also related to the level of LD decay. The plot of squared allele frequency correlations (r²) vs. genomic position (Mbp) demonstrated a decay of LD values with the genetic distance in the genomic panel (Supplementary Figure S3). The mean distance at which r² fell below 0.2 was approximately 5 Mbp. Of the 44 K SNP markers evaluated, 57% (> 25 K) were found to be polymorphic and useful for GWAS. The SNP markers displayed adequate genomic coverage, however, some gaps were observed over certain chromosomes (1 H, 2 H, and 4 H), particularly in centromere regions (Supplementary Figure S4).

Phenotype

Overall means for all traits measured are presented by year for each location in Table 1. Fusarium head blight occurred in all locations in all years. Deoxynivalenol was detected in all of the harvested grains and displayed a continuous distribution (Supplementary Figure S5). Disease symptoms were similar for Brandon and Carman, yet considerably lower at Carberry. Deoxynivalenol content varied by location and year. While the lowest FHB infection level was observed at Carberry and particularly in 2014, this location produced the highest DON content for that year (Table 1). Carman displayed the highest consistency in overall DON content.

Table 1. Means and standard errors of height (HT, cm), fusarium head blight (FHB, 0–5) and days to fifty percent heading (HD50, days) and deoxynivalenol (DON, mg kg⁻¹) content over years and locations.

Location	Year	HD50^a	Min-Max	HT^b	Min-Max	FHB^c	Min-Max	DON^d	Min-Max
Brandon	2014	57.4 ± 0.1	49.0–62.5	83.5 ± 0.4	55.5–121.5	2.0 ± 0.1	0.0–5.0	7.1 ± 0.2	0.7–30.3
	2015	53.6 ± 0.1	39.0–64.0	82.2 ± 0.4	53.0–113.0	2.4 ± 0.1	0.0–5.0	19.7 ± 0.5	4.6–90.0
Carberry	2014	71.6 ± 0.1	57.0–76.0	75.9 ± 0.4	53.0–110.0	0.4 ± 0.0	0.0–5.0	27.2 ± 0.5	2.7–69.8
	2015	61.4 ± 0.2	43.0–67.5	67.6 ± 0.3	45.0–89.5	0.7 ± 0.0	0.0–4.5	6.0 ± 0.1	0.8–19.4
Carman	2014	-	-	89.3 ± 0.4	60.0–126.5	2.6 ± 0.1	0.0–5.0	12.2 ± 0.3	1.7–57.0
	2015	-	-	81.4 ± 0.4	55.5–103.0	2.4 ± 0.1	0.5–5.0	11.7 ± 0.4	2.3–32.6

^aHD50 = days to 50% of plot headed ^bHT = height from ground to top of spike excluding awns

^cFHB = fusarium head blight score

^dDON = deoxynivalenol content

Min = minimum

Max = maximum

Heading notes were not recorded at Carman.

Relationships between traits over locations were analyzed by comparison of Pearson correlation coefficients (Supplementary Figure S6). The most robust within-trait relationship between locations was found for days to heading and plant height. Correlation of FHB and DON accumulation between locations was moderately high and moderate, respectively. With regards to between trait correlation, FHB was moderately associated with DON at Brandon (r² = 0.25, P-value < 0.0001) and Carman (r² = 0.19, P-value < 0.0001), but at Carberry the relationship was non-existent (r² = 0.00, P-value > 0.05). Height displayed a negative relationship with FHB (r² = −0.36, P-value < 0.0001), but showed no association with DON (r² = −0.04, P-value < 0.0001).

FHB reaction and DON content

Overall, there was a wide range of mean FHB reactions among genotypes, where some genotypes displayed minimal infection and others had near-maximum scores (FHB, 0.3–4.8). Likewise, the mean DON content among genotypes ranged from 4.9 to 36.9 mg kg⁻¹. Genotypes previously identified as resistant at a single location (Brandon), were confirmed in this study by testing conducted in additional environments and using alternative evaluation methods (conidial spray). A list of genotypes with consistently low FHB and/or DON content across locations is presented in Table 2. The group of landrace exotics generally held both the lowest FHB and DON content. In this group, germplasm such as CIho 4196 (Horsley et al. 2006) and Harbin (Hori et al. 2006) have been used extensively as resistance sources for breeding. However, several breeding lines with adaptive characteristicsfor production in North America demonstrated resistance to DON accumulation near those exotic accessions including genotypes such as: NDSU cultivars ‘Bowman’ and ‘Conlon’ and several Canadian breeding lines incorporating resistance from: ‘Svansota’, ‘MNBrite’, ‘Shyri’ and ‘Gobernadora’. In some cases high FHB did not result in high DON content as exemplified by HDE84194-622-1, which produced high FHB but low DON content values at all locations.

Table 2. Selection of germplasm with low Fusarium head blight (FHB) score or deoxynivalenol (DON) content. Numbers represent grand average over site years, with standard errors in brackets.

Name	Type	Mean FHB (0–5)	Mean DON (mg kg^−1)
CN 5317 (396)	exotic	0.5 (0.2)	5.7 (0.9)
CIho 14 360	exotic	0.3 (0.1)	6.0 (1.2)
Kutahya	exotic cultivar	0.3 (0.1)	6.6 (1.4)
CIho 4196	exotic	0.4 (0.1)	8.0 (1.7)
CI 2208	exotic	0.4 (0.1)	9.3 (2.1)
Harbin	exotic cultivar	0.3 (0.1)	9.7 (1.2)
HDE84194-622-1	exotic	1.1 (0.3)	6.5 (1.3)
Bowman	NDSU cultivar	4.7 (0.2)	5.2 (1.3)
Conlon	NDSU cultivar	1.4 (0.3)	6.9 (1.4)
2ND27421	NDSU line	3.0 (0.5)	4.9 (0.9)
2ND24263	NDSU line	1.9 (0.5)	7.1 (1.5)
TR06676	FCDC line	2.1 (0.5)	7.0 (1.5)
TR10692	FCDC line	1.0 (0.3)	5.8 (1.4)
SB071676	CDC line	1.7 (0.3)	6.7 (1.5)
BM0332D-280	AAFC line	0.5 (0.2)	9.5 (1.9)
BM0124-347-1	AAFC line	0.5 (0.2)	10.3 (2.1)
BM9637-53	AAFC line	0.4 (0.1)	11.1 (2.5)
CDC Bold^a	CDC cultivar	3.9 (0.4)	26.3 (5.6)
BM9925-192^a	AAFC line	2.6 (0.4)	36.9 (5.5)

^asusceptible lines included for contrast.

AAFC = Agriculture and Agri-Food Canada; CDC = Crop Development Centre; FCDC = Field Crop Development Centre; NDSU = North Dakota State University.

Association mapping

Genome-wide association study (GWAS) identified a total of 128 SNP markers (Fig. 5) associated with HD50, where approximately half (52%) of these markers were located on chromosome 2 H (Supplementary Figure S7), with several peak regions (55.0, 360.0, 501.5, 538.7, 549.6, and 635.3 Mbp). Chromosome 3 H held the second-highest number of significant markers with peaks at 28.5, 32.1, and 594.8 Mbp. A total of 14 marker associations were observed on chromosome 7 H with a peak at 55.2 Mbp. Significant marker associations were observed on all other chromosomes including peaks on chromosomes 4 H: 559.2 Mbp and 5 H: 579.8 Mbp.

Fig. 5. Manhattan plots ofall site years of 400 genotypes for days to fifty-percent heading (HD50), height (HT), FHB (Fusarium head blight score) and deoxynivalenol content (DON) over the chromosomal position. The horizontal red dashed line represents a genome-wide significance threshold (-Log₁₀ P-value ≥ 3.0). SNP markers are displayed cumulatively by chromosome (Chr) which are represented by colour.

A total of 60 SNP markers were identified associated with HT, which were identified on all chromosomes (Fig. 5). As seen for HD50, the largest number of markers (40%) were located on chromosome 2 H (Supplementary Figure S7), with association peaks (25.4, 40.5, 382.1, 425.6 Mbp). Chromosome 5 H held the second-highest number of significant markers with peaks 11.0, 499.3, and 584.1 Mbp).

Fusarium head blight (FHB) was associated with a total of 75 SNP markers (Fig. 5), located across all chromosomes (Table 3). Chromosome 5 H held the highest number of markers (23%), where peaks were observed at 25.2 Mbp (SCRI_RS_7368) and 499.3 Mbp (JHI-Hv50k-2016-321 425), where the latter coincided with a height-associated region. A region near the end of 4 HS was associated with FHB, with a peak at 2.8 Mbp (JHI-Hv50k-2016-227 104). Likewise, three peak regions were identified on 7 H including 28.2 Mbp (SCRI_RS_219709), 484.1 Mbp (JHI-Hv50k-2016-488 582), and 584.3 Mbp (JHI-Hv50k-2016-498 181). A region for FHB on 2 H at peak 538.7 Mbp (JHI-Hv50k-2016-101 173), coincided with a region associated with days to heading. Other regions associated with FHB were identified on chromosomes 1 H: 11.8 Mbp (JHI-Hv50k-2016-10 994), 3 H: 547.9 Mbp (JHI-Hv50k-2016-201 074), and 6 H: 50.9 Mbp (JHI-Hv50k-2016-386 764).

Table 3. Significant molecular markers (-Log₁₀ [P value] ≥3) associated with Fusarium head blight (FHB) averaged over all site years (n = 6).

SNP marker ID ^a	Chr ^b	Allele ^c	Position ^d	Estimate	F value	-Log10 (P value)	R^2 (%)
JHI-Hv50k-2016-10 994	1H	A/G	11796532	0.3	16.9	4.3	4.1
JHI-Hv50k-2016-55 687	1H	A/G	510489416	−0.1	11.0	3.0	2.7
JHI-Hv50k-2016-90 430	2H	A/C	357421942	0.3	11.3	3.1	2.8
JHI-Hv50k-2016-92 837	2H	A/G	355524119	0.3	10.9	3.0	2.7
JHI-Hv50k-2016-93 148	2H	A/G	382064575	0.3	12.7	3.4	3.1
JHI-Hv50k-2016-94 189	2H	A/G	416370735	0.3	12.7	3.4	3.1
JHI-Hv50k-2016-94 503	2H	A/C	425608540	0.3	12.7	3.4	3.1
JHI-Hv50k-2016-101 173	2H	G/A	538658014	0.3	15.5	4.0	3.8
JHI-Hv50k-2016-163 997	3H	G/A	35959325	−0.3	11.1	3.0	2.7
JHI-Hv50k-2016-200 089	3H	A/G	542418217	−0.2	12.7	3.4	3.1
JHI-Hv50k-2016-201 074	3H	G/A	547887089	0.2	13.5	3.6	3.3
JHI-Hv50k-2016-201 914	3H	A/G	552148262	0.2	10.8	3.0	2.7
JHI-Hv50k-2016-226 348	4H	G/A	200043	0.3	11.2	3.0	2.7
JHI-Hv50k-2016-226 367	4H	A/G	210117	0.3	11.2	3.0	2.7
JHI-Hv50k-2016-225 695	4H	A/G	2006309	0.3	11.0	3.0	2.7
JHI-Hv50k-2016-227 048	4H	A/G	2600202	0.2	12.8	3.4	3.1
JHI-Hv50k-2016-227 104	4H	G/A	2772321	0.3	20.9	5.2	5.0
JHI-Hv50k-2016-227 095	4H	G/C	2772483	0.3	19.5	4.9	4.7
JHI-Hv50k-2016-227 292	4H	C/G	3371574	−0.3	11.9	3.2	2.9
JHI-Hv50k-2016-228 556	4H	G/A	7766003	0.3	13.1	3.5	3.2
JHI-Hv50k-2016-271 609	4H	C/G	600783796	−0.2	11.4	3.1	2.8
JHI-Hv50k-2016-278 585	5 H	T/A	4185576	0.3	11.4	3.1	2.8
JHI-Hv50k-2016-278 627	5H	C/G	4290748	0.3	11.4	3.1	2.8
JHI-Hv50k-2016-282 091	5H	G/C	10988013	0.3	12.5	3.3	3.1
JHI-Hv50k-2016-282 179	5H	G/A	10993387	0.3	12.5	3.3	3.1
JHI-Hv50k-2016-283 470	5H	C/G	15019718	0.2	10.9	3.0	2.7
JHI-Hv50k-2016-284 034	5H	C/A	17366809	0.2	12.1	3.3	3.0
JHI-Hv50k-2016-284 439	5H	A/T	18246675	0.3	11.9	3.2	2.9
JHI-Hv50k-2016-286 464	5H	G/A	25035335	0.2	15.0	3.9	3.6
SCRI_RS_7368	5H	A/G	25203928	0.3	21.0	5.2	5.0
SCRI_RS_196566	5H	G/A	474102141	0.3	12.2	3.3	3.0
JHI-Hv50k-2016-312 973	5H	G/A	474233116	0.3	12.2	3.3	3.0
BOPA1_6833-658	5H	C/A	474574103	0.3	11.0	3.0	2.7
JHI-Hv50k-2016-313 362	5H	C/A	474867454	0.3	12.2	3.3	3.0
JHI-Hv50k-2016-321 425	5H	A/G	499317939	0.4	16.7	4.3	4.0
JHI-Hv50k-2016-321 503	5H	A/G	499524364	0.3	12.4	3.3	3.0
JHI-Hv50k-2016-324 079	5H	G/C	507354940	0.3	11.3	3.1	2.8
JHI-Hv50k-2016-332 519	5H	A/G	522216590	0.2	12.4	3.3	3.0
JHI-Hv50k-2016-372 853	6 H	G/A	9839643	−0.2	10.9	3.0	2.7
JHI-Hv50k-2016-386 716	6 H	A/G	50211492	0.3	10.8	3.0	2.6
JHI-Hv50k-2016-386 764	6H	G/A	50905912	0.3	12.2	3.3	3.0
SCRI_RS_219709	7 H	G/A	28232544	0.1	11.7	3.2	2.9
JHI-Hv50k-2016-468 635	7 H	G/A	71317253	0.3	11.3	3.1	2.8
JHI-Hv50k-2016-488 566	7 H	G/C	483488573	0.2	11.0	3.0	2.7
JHI-Hv50k-2016-488 570	7 H	A/G	484010690	0.2	12.3	3.3	3.0
JHI-Hv50k-2016-488 582	7 H	T/A	484148180	0.2	12.3	3.3	3.0
BOPA1_3232-201	7 H	G/C	488294965*	0.2	10.8	3.0	2.6
JHI-Hv50k-2016-488 814	7 H	G/A	489070312	0.2	11.6	3.1	2.8
JHI-Hv50k-2016-488 919	7 H	G/A	490959746	0.2	11.6	3.1	2.8
JHI-Hv50k-2016-498 095	7 H	G/A	584016335	0.3	11.1	3.0	2.7
JHI-Hv50k-2016-498 181	7 H	A/G	584342965	0.3	12.2	3.3	3.0

^aMolecular marker of 50 K SNP array (Bayer et al. 2017)

^bChr = Chromosome

^cAllele (Major/Minor)

^dGenomic position in base pairs (bp), reference genome ‘Morex v3’ gene models (Mascher et al. 2021).

* = 100% BLASTn alignment.

Deoxynivalenol content (DON) was associated with a total of 42 SNP markers (Fig. 5; Table 4), which mostly occurred on chromosome 5 H (40%). Associations were observed over multiple regions of chromosome 5 H with peaks occurring at 4.8 Mbp (SCRI_RS_114102), 34.0 Mbp (JHI-Hv50k-2016-288 841), 48.3 Mbp (JHI-Hv50k-2016-290 535), 60.8 Mbp (JHI-Hv50k-2016-291 313), and 530.9 Mbp (JHI-Hv50k-2016-337 389). Likewise, peaks were observed on chromosome 4 H at 7.7 Mbp (JHI-Hv50k-2016-228 552) and 571.1 Mbp (JHI-Hv50k-2016-261 905). Peaks were also observed on chromosome 3 H at 403.2 Mbp (JHI-Hv50k-2016-180 623) and 461.9 Mbp (JHI-Hv50k-2016-186 010). Multiple regions were also identified on chromosomes 2 H (11.4 and 18.1–19.1 Mbp) and 7 H (3.0 and 615.7 Mbp), while no significant markers were identified in association with DON on chromosomes 1 H nor 6 H.

Table 4. Significant molecular markers (-Log₁₀ [P value] ≥3) associated with deoxynivalenol (DON) content averaged over all site years (n = 6).

SNP marker ID ^a	Chr ^b	Allele ^c	Position ^d	Estimate	F value	-Log10 (P value)	R^2 (%)
JHI-Hv50k-2016-65 236	2H	G/A	11377865	1.3	13.7	3.6	3.3
JHI-Hv50k-2016-69 936	2H	A/G	18121966	1.3	10.8	3.0	2.6
JHI-Hv50k-2016-70 485	2H	C/A	19058700	1.4	11.6	3.1	2.8
JHI-Hv50k-2016-180 623	3H	G/C	403198090	−1.8	12.9	3.4	3.2
JHI-Hv50k-2016-185 349	3H	A/G	457715004	1.5	11.8	3.2	2.9
JHI-Hv50k-2016-185 889	3H	A/G	461002904	1.5	11.3	3.1	2.8
JHI-Hv50k-2016-186 010	3H	A/G	461910907	1.6	12.5	3.3	3.0
JHI-Hv50k-2016-216 594	3H	G/A	595248892	1.5	11.2	3.0	2.7
JHI-Hv50k-2016-228 552	4H	C/G	7701844	1.7	14.9	3.9	3.6
JHI-Hv50k-2016-261 905	4H	A/G	571137453	1.4	14.5	3.8	3.5
JHI-Hv50k-2016-262 043	4H	A/C	571707652	1.3	10.9	3.0	2.7
JHI-Hv50k-2016-270 157	4H	C/G	597834082	−1.4	12.4	3.3	3.0
SCRI_RS_179365	5H	G/A	4178434	1.5	11.8	3.2	2.9
JHI-Hv50k-2016-278 585	5H	T/A	4185576	1.8	14.7	3.8	3.6
JHI-Hv50k-2016-278 627	5H	C/G	4290748	1.8	14.7	3.8	3.6
SCRI_RS_114102	5H	C/A	4849576	1.7	15.0	3.9	3.6
SCRI_RS_155555	5H	G/A	4880406	1.7	13.1	3.5	3.2
JHI-Hv50k-2016-279 292	5H	C/A	5494556	1.8	13.4	3.5	3.3
JHI-Hv50k-2016-288 841	5H	G/A	34036714	2.3	13.4	3.5	3.3
JHI-Hv50k-2016-288 993	5H	G/A	35094836	2.3	13.4	3.5	3.3
JHI-Hv50k-2016-290 535	5H	G/C	48340379	1.6	12.4	3.3	3.0
JHI-Hv50k-2016-290 781	5H	A/G	55083578	2.0	11.9	3.2	2.9
JHI-Hv50k-2016-290 990	5H	C/A	56917275	2.0	11.9	3.2	2.9
JHI-Hv50k-2016-291 313	5H	A/G	60783173	2.1	12.3	3.3	3.0
JHI-Hv50k-2016-337 389	5H	A/T	530890859	1.7	20.6	5.1	4.9
JHI-Hv50k-2016-337 660	5H	G/A	531560545	1.6	12.3	3.3	3.0
JHI-Hv50k-2016-437 412	7H	G/A	2991267	−1.5	12.0	3.2	2.9
JHI-Hv50k-2016-512 584	7H	G/A	615711190	0.9	12.3	3.3	3.0
BOPA1_8758-564	7H	A/G	615732324*	0.9	12.3	3.3	3.0

^aMolecular marker of 50 K SNP array (Bayer et al. 2017)

^bChr = Chromosome

^cAllele (Major/Minor)

^dGenomic position in base pairs (bp), reference genome ‘Morex v3’ gene models (Mascher et al. 2021).

* = 100% BLASTn alignment.

Candidate genes

I) Days to heading (HD50) and height (HT). Marker-traitassociations were often coincidental for HD50 and HT. Plant height (HT) was associated on chromosome 2 H (25.3 Mbp) and borderline (-Log₁₀ P-value = 2.9) for HD50. Photoperiod response regulator gene (HORVU.MOREX.r3.2 HG0107710, Ppd-H1),which determines the sensitivity to long days is located at this position (Turner et al. 2005). Similar association patterns over locations corroborate the influence of chromosome 2 H onthese traits.

Plant height (HT) was associated with chromosome 1 H at 479.1 Mbp in proximity to HORVU.MOREX.r3.1 HG0077240(477.8 Mbp), Protein FLOWERING LOCUS T, HvFT3, Ppd-H2. Days to heading (HD50) showed a marker-trait-association inchromosome 2 H at 538.7 and 635.3 Mbp which are respectively associated with HORVU.MOREX.r3.2 HG0177270 (Hv FTL5, Flowering locus T) and HORVU.MOREX.r3.2 HG0204770 (HvAP2, ethylene-responsive transcription factor). Both HD50 and HT showed a marker-trait-association on 2 H: 403–407 Mbp region, which are respectively associated with HORVU.MOREX.r3.2 HG0153260 (early-responsive to dehydration stress protein, HvERD4) and HORVU.MOREX.r3.2 HG0153990 (light-dependant short hypocotyls protein). Days to heading (HD50) was associated with positions on chromosome 3 H: 28.5and 32.3 Mbp. Cytochrome P450 superfamily protein (HORVU.MOREX.r3.3 HG0232150) is located at 28.5 Mbp, whereas Six-rowed spike4 [Vrs4, syn. Intermedium-E], HORVU.MOREX.r3.3 HG0233930 (LOB domain-containing protein) is located at 33.7 Mbp. Vrs4 is a major morphometric gene that controls the two-row spike morphology through Vrs1 (HORVU.MOREX.r3.2 HG0184740, 570.8 Mbp) and its suppression of lateral spikelet fertility (Koppolu et al. 2013).

Plant height (HT) also showed association with several nearby regions on chromosome 5 H, including 499.3 Mbp nearby HORVU.MOREX.r3.5 HG0498800 (Protein STAY-GREEN, chloroplastic) andHORVU.MOREX.r3.5 HG0498830 (Pseudo response regulator, PRR95). Both traits also demonstrated association on chromosome 5 H at 585.8 Mbp in the proximity (583.9 Mbp) to gibberellin 20 oxidase (HORVU.MOREX.r3.5 HG0536610, HvGA20ox1).

II) Fusarium head blight. Genome-wide association study (GWAS) of FHB identified nucleotide-binding site leucine-rich repeat (NBS-LRR) class disease resistance proteins association in two regions of chromosome 1 H: 11.6 (HORVU.MOREX.r3.1 HG0005640) and 510.4 Mbp (HORVU.MOREX.r3.1 HG0093120). Multiple locations on chromosome 2 H were associated with multiple protein kinase proteins: HORVU.MOREX.r3.2 HG0153670 (357.4 Mbp); HORVU.MOREX.r3.2 HG0156970 (382.1 Mbp),and HORVU.MOREX.r3.2 HG0177510 (538.2 Mbp). Another association was observed on chromosome 2 H at 425.6 Mbp in proximity to HORVU.MOREX.r3.2 HG0161420 (MLO-like protein). On chromosome 3 H, FHB at 359.6 Mbp adjacent to HORVU.MOREX.r3.3 HG0234370 (NBS-LRR resistance protein) and at542.4 Mbp was associated with HORVU.MOREX.r3.3 HG0300850 (Cytochrome P450).

On chromosome 4 H, FHB was associated at 2.1 Mbp close to HORVU.MOREX.r3.4 HG0331940 (Copper ion-binding protein); HORVU.MOREX.r3.4 HG0331250 (cellulose synthase-like protein), andHORVU.MOREX.r3.4 HG0331260 (ABC transporter ATP-binding protein ARB1). Association was also seen with the adjacent region 2.8 Mbp in the position of genes HORVU.MOREX.r3.4 HG0332220 (GDSL esterase/lipase)and HORVU.MOREX.r3.4 HG0332240 (Fatty acyl-CoA reductase 3).

On chromosome 5 H, FHB-association at 25.2 Mbp wasin proximity to Chaperone heat shock protein DnaJ (HORVU.MOREX.r3.5 HG042960). Also associated with FHB were regions: 474.1 Mbp, HORVU.MOREX.r3.5 HG0489730 (UDP-Glycosyltransferase superfamily protein), and HORVU.MOREX.r3.5 HG0489810 (caffeoyl-CoA 3-O-methyltransferase); 499.3 Mbp, HORVU.MOREX.r3.5 HG0498800 (Protein STAY-GREEN, chloroplastic),and 507.4 Mbp HORVU.MOREX.r3.5 HG0501610 (Glutathione S-transferase). Other detected associations with genes were chromosome 6 H at50.9 Mbp, HORVU.MOREX.r3.6 HG0556730 (LRR protein kinase), and chromosome 7 H at484.0 Mbp HORVU.MOREX.r3.7 HG0712930 (NBS_LRR resistance protein), and 584.3 Mbp, HORVU.MOREX.r3.7 HG0732460 (Protease inhibitor/seed storage/lipid transfer family protein).

III) Deoxynivalenol (DON) content. Multiple DON associations were observed on chromosome 2 H at11.4 Mbp, HORVU.MOREX.r3.2 HG0100300 (Peroxidase);18.1 Mpb HORVU.MOREX.r3.2 HG0104570 (LRR protein kinase); and 19.1 MbpHORVU.MOREX.r3.2HG0104860 (Expansin). On chromosome 3 H associations occurred at 403.2HORVU.MOREX.r3.3 HG0276410 (Trichome birefringence-like protein);457.7 Mbp HORVU.MOREX.r3.3 HG0284050(IQ-domain 6, a calmodulin-binding protein);461.9 Mbp HORVU.MOREX.r3.3 HG0284810 (WRKY DNA-binding protein 50), and 595.2 Mbp HORVU.MOREX.r3.3 HG0319600 (Serine/Threonine-kinase).

DON was associated with several genomic regions of chromosome 5 H, 4.2–5.5 Mbp which were coincident for association with FHB including genes: HORVU.MOREX.r3.5 HG0420980 (Early-responsive to dehydration stress protein, ERD4), HORVU.MOREX.r3.5 HG0421020 (Pathogenesis-related protein 1, PR1), HORVU.MOREX.r3.5 HG0421310(Leucine-rich receptor-like protein kinase family protein) and HORVU.MOREX.r3.5 HG0421350 (Neutral ceramidase).Several associations were seen between 34.0–60.8 Mbp were in the region of CBS domain-containing protein (HORVU.MOREX.r3.5 HG0431770);HORVU.MOREX.r3.5 HG0433910(WRKY DNA-binding protein 74); HORVU.MOREX.r3.5 HG0436650 (DA1, ubiquitin-activated limiter of seed size);; HORVU.MOREX.r3.5 HG0437600 (Myb family transcription factor-like protein). Mapping to 530.9 Mbp, a SNP (JHI-Hv50k-2016-337 389) explained 4.9% of the phenotypic variance (-Log₁₀ P-value = 5.1), associated with HORVU.MOREX.r3.5 HG0512320 (Cytochrome P450 superfamily protein).

Discussion

The assembled panel of genotypes displayed a wide range of FHB infections in this study, which resulted in varying levels of DON accumulation. Deoxynivalenol (DON) production was positively associated with Fusarium infection, and as seen in other studies the relationship is only moderate (Buerstmayr et al. 2004; Choo et al. 2004; He et al. 2015; Tucker et al. 2019). Later infections which do not produce high FHB levels, may result in very high mycotoxin accumulation. FHB infection at Carberry in 2014, which had the lowest level of infection, yet some of the highest DON levels of the study. This demonstrates the complexities in the relationship between FHB and DON, both of which are highly subjective to environmental conditions at their respective, developmental stages. This also illustrates the essential need to conduct DON testing in barley, which can occur in seemingly uninfected kernels.

Patterns of relatedness were analyzed for germplasm of the genomic panel, where clustering of materials occurred mainly by barley end-use i.e. malting or livestock feed. The genetic variance was narrower within the malting class of barley. This is to be expected given the restrictions posed by malting quality standards. As a result, use of elite crosses is required to maintain Canadian malting quality. In addition to type, the genotypes also demonstrated clustering by institution. However, a substantial level of genetic transfer was apparent through the documented admixture within breeding pools of the Canadian breeding programs. The STRUCTURE analysis identified two main, genetically distinct sub-populations. Modern malting barley cultivars such as ‘AAC Connect’, ‘AAC Synergy’, and related breeding lines from the AAFC-BRDC program, clustered distinctly from the remaining entries in the genomic panel.

Fusarium head blight (FHB)-marker association was observed on chromosome 1 H at 11.5 Mbp. A QTL (Qlds.1 H) has been previously detected by Bedawy et al. (2018) on chromosome 1 H 22.1 Mbp and associated with SCRI_RS_239784 through greenhouse spike inoculation and F. graminearum leaf assay. Likewise on chromosome 1H, an association was observed at 510.5 Mbp, which coincided with a previously identified QTL for FHB (BOPA1_2572-986, 506.6 Mbp) in two-row barley (Mamo and Steffenson 2015). An association with FHB was observed at 357.4 Mbp close to leucine-rich receptor (LRR)-like protein kinase, HORVU.MOREX.r3.2HG0153670. Aprevious association with FHB was reported between this region and a QTL, FHB2H.50–56 (bin 6–7, Massman et al. 2011) along withSNP marker SCRI_RS_130072 (356.9 Mbp) for DON content (Mamo and Steffenson 2015). Interestingly, FHB symptoms were not associated with the Ppd-H1 region in this study, as previously reported in bi-parental QTL studies (Huang et al. 2018; Ogrodowicz et al. 2020, QTL termed QFHB.IPG-2H_1).

Genome-wide association study (GWAS) showed several regions of chromosome 2H between 355.5–382.1 Mbp associated with FHB, where several markers were previously associated with DON content (SRCI_RS_230984, 356.2 Mbp; BOPA1_1062-478, 359.5 Mbp; SCRI_RS_197997, 360.7 Mbp) (Mamo and Steffenson 2015). In this study, this region (355.4–360.0 Mbp) was also associated with HD50 and HT. Mildew Locus O (MLO) proteins like HORVU.MOREX.r3.2HG0161420 (425.6 Mbp) are susceptibility factors. Mutant forms of these genes (mlo) confer resistance to powdery mildew fungus but enhances susceptibility to hemi-biotrophic organisms. On chromosome 4 H, FHB was associated with multiple SNPs between 0.2–7.8Mbp. A QTL DON-QTL-4 H-12_11485 associated with DON resistance was previously identified in proximity at 6.4Mbp (Mamo and Steffenson 2015) and associated with QTL DON4H.03 (bin 1) (Massman et al. 2011). Several genes with potential resistance contributions in this region including HORVU.MOREX.r3.4 HG0331250 (cellulose synthase-like) HORVU.MOREX.r3.4 HG0331260 (ABC transporter) and HORVU.MOREX.r3.4 HG0332240 (Fatty acyl-CoA reductase 3). A GDSL esterase/lipase gene, HORVU.MOREX.r3.4 HG0332220 associated with JHI-Hv50k-2016-227 104 was highly significant (-Log₁₀ P-value = 5.2).

Resistance factors can be induced or possess a constitutive basis (Kushalappa et al. 2016). Lipids are important regulators of plant defence, where the well-known defence hormone, jasmonic acid (JA) is a derivative of lipid-metabolism. Genes contributing to lipid catabolism such as GDSL esterase/lipases and non-specific lipid-transfer proteins (nsLTP, HORVU.MOREX.r3.7 HG0732460), have been implicated with FHB resistance in transcriptional studies of cereals (Foroud et al. 2012; Gottwald et al. 2012; Huang et al. 2016; Tucker et al. 2021). McLaughlin et al. (2021) demonstrated the importance of nsLTP genes, through the use of transgenics AtLTP4.4 from Arabidopsis in wheat, which suppressed DON production in the field, where they hypothesized it reduced DON-induced ROS. Fatty acyl-CoA reductase 3 (HORVU.MOREX.r3.4 HG0332240) which is involved in cuticular wax production, may be an important barrier for barley resistance to F. graminearum (Kumar et al. 2016). Caffeoyl-CoA 3-O-methyltransferase (HORVU.MOREX.r3.5 HG0489810), is involved in the production of feruloyl-CoA lignin which contributes to the fortification of the cell wall. Structural-related proteins were identified which can aid fungal penetration such as expansin (HORVU.MOREX.r3.2HG0104860) which causes cell loosening or cellulose synthase (HORVU.MOREX.r3.4 HG0331250) and UDP-glucose 6-dehydrogenase (HORVU.MOREX.r3.5 HG0489730) leading to hemicellulose and pectin.

On chromosome 5 H, FHB association was seen at 25.2 Mbp (SCRI_RS_7368, -Log₁₀ P-value = 5.2). Fusarium head blight (FHB) association was seen over the range 474.1–474.9 Mbp. A major QTL (Qsds.5 Hc) explaining 39% of genetic variance was identified in a greenhouse study in association with SNP marker SCRI_RS_174710 located on 5 H at 523.9 Mbp (Bedawy et al. 2018). This association was in proximity to markerJHI-Hv50k-2016-332 519 (522.2 Mbp). Also associated was JHI-Hv50k-2016-321 425 (-Log₁₀ P-value = 4.3) at 499.3 Mbp beside HORVU.MOREX.r3.5 HG0498800HORVU5Hr1G08121 (Protein STAY-GREEN) and 507.4 Mbp beside HORVU.MOREX.r3.5 HG0501610 (Glutathione S-transferase). FHB association was seen on chromosome 6 H at 50.2–50.9 Mbp in the same region of QTL for FHB, FHB6H.42–61,and DON content DON6H.42–67 (bins 5–7, Massman et al. 2011)and corresponding gene HORVU.MOREX.r3.6 HG0556730 (LRR protein kinase).

Regions of association with DON did not always coincide with those of FHB and showed associations on chromosomes 2H- 5 H and 7 H. Deoxynivalenol (DON) associations on 3 H occurred at numerous regions independently of FHB, which corresponded to genes of interest: 403.2 Mbp (52.0 cM) – HORVU.MOREX.r3.3 HG0276410 (Trichome birefringence-like protein); 457.7 Mbp (57.1 cM) HORVU.MOREX.r3.3 HG0284050 (IQ-domain 6, a calmodulin-binding protein) and 461.9 Mbp (57.4 cM) – HORVU.MOREX.r3.3HG0284810 (WRKY DNA-binding protein 50) in a region with QTL DON3H.52–65 identified with DON content (Massman et al. 2011). A DON-association was also observed at 595.2 Mbp which was coincident with the association for FHB and the number of days to heading. Deoxynivalenol (DON) association for overall analysis was identified in the tip of chromosome 4 HS at 7.7 Mbp which corresponds with the FHB association at the region (0.21–7.8Mbp) and aforementioned DON QTL. Deoxynivalenol (DON) also showed an association on chromosome 4 H at 571.1 Mbp in the region of HORVU.MOREX.r3.4 HG0405710(Ethylene-responsive transcription factor 1, ERF1). On chromosome 5 H, multiple regions were associated with DON content. Several associations were seen between 4.2–5.5 Mbp, with a peak at 4.9Mbp (-Log₁₀ P-value = 3.9). These regions overlapped those of FHB association and include several receptor genes along with HORVU.MOREX.r3.5 HG0421020 (Pathogenesis related 1, PR1).Increased expression of salicylic acid (SA) and SA-inducible PR1 protein correlated with lower FHB severity in wheat (Makandar et al. 2012). A large association was observed at 530.9 Mbp (-Log₁₀ P-value = 5.1) in the position of HORVU.MOREX.r3.5 HG0512320 (Cytochrome P450 superfamily protein). A QTL for FHB severity was previously reported at this location, QFHB.IPG-5 H close to marker SCRI_RS_184066 (529.5 Mbp) (Ogrodowicz et al. 20 020).

Several leucine-rich receptor-like protein kinase (LRR-RLK) family proteins were identified in association with FHB or DON (HORVU.MOREX.r3.2HG0104570, HORVU.MOREX.r3.2HG0156970, HORVU.MOREX.r3.5 HG0421310, HORVU.MOREX.r3.6 HG0556730). LRR-RLK proteins are broad-spectrum, pattern-recognition receptor (PRR) proteins localized to cell surfaces and are an important component of the surveillance system for the detection of microbial signatures playing an important part in the initial stage of plant defence (Tang et al. 2017). These play a role in pattern-triggered immunity (PTI) signaling by recognition of pathogen-associated molecular patterns (PAMPs). Following PAMP recognition, PTI response activates MAP kinases and reactive oxidative species (ROS) production. To circumventPTI resistance response, pathogens introduce effector molecules to the host. Plants have evolved a secondary layer of defence, called effector-triggered immunity (ETI) which involves specialized, resistance (R)-proteins (Jones and Dangl 2006).Nucleotide-binding site leucine-rich repeat (NBS-LRR) are R-proteins that interact (directly or indirectly) with pathogen effector molecules (alternatively called avirulence [Avr] proteins) whichtrigger anETI defence response, associated with an oxidative burst, hypersensitive response, and induction of pathogenesis-related genes. Several NBS-LRR genes were identified as being associated with FHB resistance (HORVU.MOREX.r3.1HG0005640, HORVU.MOREX.r3.1HG0093120, HORVU.MOREX.r3.3HG0234370).The number of marker-trait associations nearby genes involved in pathogen-sensing aspects suggests they may be important factors in barley-Fusarium interactions.

Transcription factors are important factors in modulating transcriptional programming in immune response (Tsuda and Somssich 2015). Genes of interest included those carrying myeloblastosis (MYB) domain (HORVU.MOREX.r3.5 HG0437600), which are key regulators of phenylpropanoid metabolism (Liu et al. 2015). Also identified by this study as a gene of interest was ethylene-responsive transcription factor 1 (ERF1, HORVU.MOREX.r3.4HG0405710) where ethylene signaling is known to be highly involved in FHB outcomes albeit complex (Foroud et al. 2019). Barley cultivars have demonstrated differential response in WRKY transcription factors (Huang et al. 2016; Tucker et al. 2021), with some contributing to resistance while others to susceptibility to F. culmorum crown rot. (Uluhan et al. 2019). Karre et al. (2019) demonstrated HvWRKY23 as a major component of FHB through regulation of flavonoid production. HORVU.MOREX.r3.3HG0284810 (WRKY DNA-binding protein 50) and HORVU.MOREX.r3.5 HG0433910 (WRKY DNA-binding protein 74) are both genes of interest identified in the study. Arabidopsis gene AtWRKY50 is associated with salicylic acid (SA) defense response through binding to the PR1 promoter region and activating its expression (Hussain et al. 2018).

While barley genotypes differed in their response to FHB and DON accumulation, as expected, none showed immunity, or the ability to prevent DON production entirely. Defence strategies may involve limiting DON creation, but also involve chemical modification toxins, which may reduce their toxicity, but also may increase mobility for storage in vacuoles or cell walls (Boutigny et al. 2008). Several gene classes have been commonly connected to DON-specific resistance and/ or detoxification (Gunupuru et al. 2017) including those in this study: Cytochrome P450 (HORVU.MOREX.r3.5 HG0512320), ABC transporters (HORVU.MOREX.r3.4HG0331260), GSTs (HORVU.MOREX.r3.5 HG0501610) and UGTs (HORVU.MOREX.r3.5 HG0489730).

Numerous marker-trait associations for both FHB (51 SNPs over 22 regions) and DON (29 SNPs over 13 regions) were identified in this study, however, no regions could be defined as a major QTL. Low correlations for marker-trait associations are typical of GWAS. This study employed the 50 K SNP chip which provided sufficient marker coverage at > 25 K polymorphic markers. This finding, along with the observed continuous distribution of FHB-related traits, implies that a large number of minor genes likely contribute to the overall phenotype of two-row barley supporting the polygenic nature of FHB in barley. Several regions also coincided with HD50 or HT, implying that correlations between traits may be the result of pleiotropy, however, gene functions may imply tight linkages. Due to small additive effects and/or percentage of total variation explained by associated regions, genomic regions of interest may have limited breeding value in MAS, but could be monitored. Thus a more holistic approach should be employed when breeding FHB-resistant two-row barley, such as genomic selection (Meuwissen et al. 2001) which considers the joint merit of all molecular markers to predict genomic estimated breeding values (GEBVs) of breeding lines. Molecular markers identified in this study could provide additional information for developing prediction models. However, singular SNP molecular markers could not adequately predict FHB and DON, where phenotypic evaluations are still the best predictors of these traits even though they are time consuming and dependent on environmental conditions.

Acknowledgements

We thank Aaron Beattie (Crop Development Centre) and Patricia Juskiw (Field Crop Development Centre) for providing seeds. We thank Andy Tekauz and Jeannie Gilbert (Agriculture and Agri-Food Canada (AAFC) Cereal Research Centre) for providing F. graminearum single-spore isolates for inoculum production.We acknowledge the technical assistance of Nathan Henderson, Kevin Moore, Mira Popovic and the wheat-breeding group at AAFC-BRDC.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online here: https://doi.org/10.1080/07060661.2022.2086925

Data availability statement

The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus (Edgar 2002) and are accessible through GEO Series accession number GSE188791 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE188791). LSMeans of phenotypic data are also available (Supplementary Table S1).

Additional information

Funding

This research was funded by the Brewing and Malting Barley Research Institutes (BMBRI), Manitoba Crop Alliance, Growing Innovation: Agri-Food Research and Development Initiative (ARDI) project to D. Fernando and A. Badea. Support was also delivered through part of the Barley Cluster project led by Alberta Barley with funding from the Western Grains Research Foundationand AAFC under the Growing Forward 2 program. Funding for the University of Manitoba FHB nursery came from the WGRF core agreement.