The diversity panel employed in this study comprises 293 Ethiopian durum wheat landraces and 25 Ethiopian durum wheat improved varieties (Table S1). The panel was assembled to represent Ethiopian durum wheat diversity. Landraces were obtained from the Ethiopian Biodiversity Institute (www.ebi.gov.et) in form of seeds. During the prior seed amplification, each accession was cleaned to obtain a reference genotype that was used for genomic DNA extraction and for field experiments downstream. Full details about the diversity panel assembly procedure can be found in Mengistu et al. (2016). Genomic DNA was extracted using the GenElute Plant Genomic DNA Miniprep Kit (Sigma-Aldrich, St Louis, MO) in Mekelle University, Mekelle, Tigray region, Ethiopia. After the check for quality and quantity, the DNA was typed with the 90K SNP wheat array (Wang S. et al., 2014) at TraitGenetics GmbH (Gatersleben, Germany). Samples not providing genotypic data of sufficient quality were discarded from the molecular analyses. SNPs were filtered for failure rate below 20% and for minor allele frequency above 5%. Details about the genotyping procedure can be found in Mengistu et al. (2016). Complete genotyping data for the varieties included in this study is reported in Table S9.

The experiment was conducted under natural infestation during 2012 and 2013 growing seasons at Geregera, Meket district (Wollo, Amhara region, 11°4′N/38°52′E; WGS 84 coordinates), at an altitude of 2,876 m above sea level. The experiment was conducted in rainfed conditions. Sowing was performed after the onset of the main rainy season, on July 5, 2012, and July 9, 2013, respectively. The field was laid out in a partially balanced lattice design, with replications of block size 20 rows by 20 columns, block length 49.5 m. The genotypes sown but not used for this study were considered as fillers. This design was used to allow the adjustment of treatment means for block effect as well as to provide effective control within replicate variability. The seed rate used was 85 kg ha⁻¹, with seeds drilled evenly in rows. Recommended rate of fertilizer for wheat for Geregera area were used, with 64 kg N (Urea) and 46 kg P₂O₅ (DAP) per hectare. P₂O₅ was entirely applied at planting time, while N was applied 1/3 at planting and 2/3 at tiller initiation. Weeding was done three times by hand. In 2012, the plot size was 2 × 0.6 m, and each plot had three rows with 0.2 m spacing. In 2013, the plot size used was 2 × 0.8 m, and each plot had four rows with 0.2 m spacing. The spacing between rows and replications were 0.5 and 1.5 m, respectively. Five plants in the middle rows were randomly selected and tagged to be used consistently to record STB infection. The same field was used to collect agronomic and phenologic traits.

The Septoria disease severity (SDS) was scored visually, according to a double-digit scale (00–99) modified from Saari and Prescott (1975) for wheat foliar diseases. The assessment was taken in five randomly selected productive plants. The evaluation was carried out individually at two time points: the heading stage and the maturity stage. For each score, the percentage of disease severity was estimated based on Equation (1), following the formula used by Sharma and Duveiller (2007): (1)SDS=[(D19)(D29)]100 where D1 represents the vertical disease progress as deriving from the average relative height reached by the disease as recorded from five random plants (0–9). D2 represents the severity of the disease, measured as the average relative coverage of the diseased leaf area (i.e., the necrotic leaf area) recorded from the upper four alive leaves in the same five plants. The SDS index is thus composed by a first digit representing the blotch development up the plant height (e.g., 5 if the disease reached the mid-point of the plant or 50%, 8 if it reached the flag leaf, 9 if it reached the spike), and a second digit representing severity per se (e.g., 1 for 10% to 9 for 90%). SDS values range from 0 to 100, where 0 would indicate complete resistance, and 100 would indicate complete susceptibility. The SDS trait was sided by the Septoria progress coefficient (SPC) trait, calculated as in Equation (2) following the formula from Eyal and Ziv (1974): (2)SPC=(SDHPH) where SDH (Septoria disease height) is the maximum height from the ground where pycnidia of the pathogen are found on the plant, in cm. PH is the height of the accessions, averaged over the same plants used to derive SDH. The coefficient indicates the position of pycnidia relative to plant height, regardless of pycnidial coverage, and allows the comparison of infection placement on cultivars with different plant stature. SPC is thus a ratio indicating the height reached by the disease. SPC of 0.0 would mean no disease at all, while SPC of 1.0 would indicate that the disease is covering the entire plant.

Ten traits accounting for phenology, yield, and yield components were investigated. Days to flowering and days to maturity were recorded for whole plots when 50% of the plants reached the corresponding Zadock's growth stage. The full plot was used to measure grain yield (as grams of grain produced per plot, then projected to t/ha), biomass (as the dry weight of the above-ground harvested biomass per plot, then projected to t/ha) and 1,000 grain weight (as the weight of 1,000 kernels, in grams). The moisture content of grains was measured using a digital grain moisture meter and used to adjust grain yield and biomass to 12.5% moisture content. Three randomly selected plants per plot were used to measure tillering (as the number of effective/productive tillers per plant), plant height (in cm), spike length (distance between the pedicule base and the tip of the spike excluding awns, in cm), number seeds per spike, and the number of spikelets per spike. For further details about the agronomic data collection, see Mengistu et al. (2016).

All phenotypic data was normalized using arcsine square root transformation. The phenotypes were analyzed using a linear mixed model (LMM) including genotypes as fixed effects. The effects of year, block and the interaction of genotypes by year were treated as random effects. The parameters of the model were estimated by the method of restricted maximum likelihood (REML) using SAS statistical software (SAS Institute, Cary, NC). BLUE values were calculated from replicates of each line for each trait for the different trials. Test of X² goodness of fit was applied to investigate if there was normality performance of each trait at each trial. BLUE values are referred to as “combined phenotypes” throughout the manuscript text. Phenotype means were transformed back to percentages for discussion. The software R (R Core Team, 2013) was used to conduct the analyses and produce plots with custom script available upon request. Histograms were plotted using R/ggplot2 (Wickham, 2009), correlations were plotted using R/corrplot (Wei, 2013).

Association mapping was conducted using R/GAPIT (Lipka et al., 2012). The initial 81,584 markers available on the array have been filtered to produce the working SNP set used in downstream analyses. Markers failing in over 20% of the samples were removed. Markers having a minor allele frequency below 5% were also filtered. Map position of markers was derived from the consensus durum wheat genetic map (Maccaferri et al., 2015). Pairwise LD measures (r²) were obtained for all markers in each linkage group using R/LDheatmap (Shin et al., 2006). To calculate LD decay, only markers pairs within 50 cM were considered. Mean LD for pairs was calculated for markers at increasing genetic distances with a sliding window and plotted against the genetic distance with custom R scripts available upon request.

The working set of markers was input together with STB resistance phenotypes combined over the 2 years into a mixed linear model (MLM) accounting for uneven relatedness among samples. This method effectively controls population structure to lower type I errors. Kinship was calculated in R/GAPIT following VanRaden's method (VanRaden, 2008). Principal components (PC) describing the genetic structure of the panel were calculated with R/GAPIT and iteratively added to the fixed part of the model, from PC1 to PC10. The SUPER compression model was used Wang Q. et al. (2014). The best fit of the model was visually evaluated on quantile–quantile plots. Marker-trait associations (MTA) were deemed highly significant when surpassing a threshold calculated using the Bonferroni method and accounting for multiple testing at a nominal p-value of 0.15. GWA scan providing significant MTAs according to this threshold are the sole reported. The false discovery rate (FDR) was computed using Storey's method (Bass et al., 2015) and was reported with GWA results as an alternative, less stringent multiple testing correction method. In order to provide a better description of the genomic locations associated with SDS and SPC, the same GWA model was run on data collected each year separately and on agronomic trait values combined over the 2 years. Allelic states at significant markers were converted to the arbitrary numeric values (−1, homozygous for the most frequent allele; 0, heterozygous; +1, homozygous for the less frequent allele). R² values for significant marker tests were produced regressing numeric allele states at the markers over the corresponding phenotypes with a linear model. The regressions coefficient was then used to derive positive alleles at each significant MTA.

Significant markers without map positions were traced back on their genetic positions by using their numerical encoding as phenotypes in a GWA scan using a compressed linear model in R/GAPIT. The deriving MTA denote genomic regions in LD with the unmapped marker used as phenotype, and thus indicate the likely genetic position of the latter. Custom R scripts available upon request were used. A putative QTL is defined as a genomic location identified by one or more MTAs.

KASP markers were designed from SNP markers underlying putative QTL for STB. When possible, multiple markers per QTL were developed. KASP markers were designed and tested on a set of 46 Ethiopian samples belonging to the durum wheat reference collection (DWRC) at LGC (http://www.lgcgroup.com/). Working KASP markers are reported in Table S8.

S. tritici blotch (STB) natural infestation was slightly different over the 2 years. An ANOVA showed that the genotype effect was significant in each year, and that the year effect was significant when considering combined data, except for SDS at maturity (Table 1). The interaction of year by genotype was significant in all the disease traits analyzed. The distribution of Septoria disease resistance was similar among landraces and MVs for the three traits considered (Figure 1). The distribution of the traits is pseudo-normal, confirming the quantitative nature of Septoria resistance in the Ethiopian material. SDS was expectedly higher at maturity than at heading stage. Ethiopian durum wheat landraces are extremely diverse in their resistance to STB (Table S1). SDS at heading ranges from 3.38 to 43.97% in landraces, whilst in MVs range from 8.45 to 28.22%. A hundred and thirty-five landraces (46.1%) and 13 MVs (52%) show a disease severity below 15% at heading. At this stage, however, only two MV show SDS below 10%, while 51 landraces do. At maturity, the range for SDS ranges from 7.64 to 55.95% in landraces and from 5.65 to 48.96% in MVs. Among landraces, 135 (40.3%) fall within 20% SDS, while five MVs (50%) do. Only one MV but 43 landraces showed SDS at maturity below 15%. The range of SPC is broader for landraces (0.43–0.81) than for MVs (0.45–0.72). In 2012, the MV with the lower SDS at heading was Mangudo (8.45), while the MV with lower SDS at maturity was Selam (5.65) (Table S1). Selam was the MV having the lower SDS at heading and maturity in 2013, and together with Mangudo was the most resistance MV across the 2 years' data in disease severity at maturity and heading, respectively. Tossa was the most susceptible line for SDS at both heading and maturity across the 2 years. In many cases, the same varieties had different affection severities across the 2 subsequent years. This can be seen comparing SDS and SPC values for the same genotype across the two seasons (Figure S1). Although the pressure of disease at heading and maturity was generally lower in 2013 than in 2012, several varieties showed contrasting trends resulting in higher SDS values in 2013 than in 2012, especially at the maturation stage. Conversely, SPC shown a generalized increase in 2013 as compared to 2012, with most of the varieties showing increased affection values over the 2 years (Figure S1).

Disease resistance traits were correlating with some of the agronomic traits measured in the same experimental field (Table S2, Figure 2A). The phenology of the plant is linked with the progression of the disease, as reported by a significant negative correlation between SDS and days to flowering and days to maturity. Grain yield and biomass, conversely, show low correlations with all disease traits, and spike length is not significantly correlated to any of those. This result shows that in our panel those genotypes with the highest productivity also suffer the highest STB affection. The number of seeds per spike and the grain size (reported as thousand grain weight) are inversely correlated with SDS (Table S2). A principal component analysis (PCA) was used to visualize the relation between disease and agronomic traits in the collection (Figure 2B). The first two PCs cumulatively explained 46.6% of the total variance. Genotypes having a shorter flowering time also have a higher infection, although the causal direction of this relationship cannot be derived from the data. This could be contributed by disease escape mechanisms, i.e., late flowering varieties escaping the disease spread and appearing more resistant as a consequence (Arraiano et al., 2009). Measures of yield and yield components are only poorly correlated to SPC and SDS.

The Ethiopian durum wheat panel was genotyped with the wheat 90K array (Wang S. et al., 2014). After filtering for sample quality, 300 accession were retained for the diversity analysis and the GWA. After filtering for data quality, 16,223 polymorphic SNP markers were retained and used to evaluate the genetic structure of the panel and to conduct a GWA study on the measured phenotypes. The PCA and kinship analysis conducted on molecular data showed that the durum wheat cultivated in Ethiopia is clearly divided into two groups. The group of MVs is highly similar within itself and different from that of landraces (Figure 3A). Landraces, on the other hand, are highly admixed and uniformly differentiated from MVs (Figure 3B). If not controlled for, this structure could inflate GWA statistics.

The GWA was conducted on STB resistance data combined over the 2 years, then on STB data collected each year separately. The GWA scan was performed also on agronomic data combined over the 2 years. We report and discuss marker-trait associations (MTA) surpassing the Bonferroni threshold. Individual marker p-values and FDR corrected q values for GWA scan on STB are reported in Table S3. Quantile-quantile plots for selected models are reported in Figure S2 for STB resistance data and in Figure S3 for agronomic data. Overall, 35 MTAs were identified with 32 unique SNP markers, of which 13 did not have a chromosomal position on the durum wheat genetic map (Maccaferri et al., 2015). The position of unmapped markers providing significant associations was derived exploiting their linkage disequilibrium (LD) with markers having a map position: for disease data, they clustered with MTAs already reported on chromosome (Chr) 1A and Chr 3A. A previously unmapped MTA for SDS at heading recorded in 2012 was mapped to Chr 5B (Table S4). The R² for MTAs varied substantially among traits and markers, the highest being 10.6% for an SPC MTA on Chr 3A. SDS at heading was the trait providing the second highest R², at 8.7% for and MTA on Chr 2B (Table S5). The MTAs were combined according to their genetic position to identify putative QTL for STB resistance. In order to identify a window of genetic distance in which to combine MTAs, a pairwise LD analysis was conducted using the SNP markers produced on the panel. The LD in the Ethiopian durum wheat was dropping below r² = 0.2 (representing null LD) within 1.5–6 cM in all chromosomes (Table S6, Figure S4), denoting a relatively high mapping definition. MTAs falling on the same linkage group within the genetic distance for LD decay specific for that chromosome were assigned to the same putative QTL, giving a total of five putative QTL identified (Table 2, Figure 4). A putative QTL deriving from the MTA with estimated genetic position on Chr 5B is discussed separately. Diagnostic KASP markers for qSTB.1, qSTB.3, qSTB.4, qSTB.5, and for the putative QTL on Chr 5B are reported in Table S8.

Combined measures for SDS at heading identified one MTA pointing to a putative QTL at 67.7–69.2 cM on Chr 1A (qSTB.1) and one MTA at 85.8 cM on Chr 2B (qSTB.2) (Table 2, Figure 4, Figure S5). One additional unmapped MTA for SDS at heading was mapped to qSTB.1 (Table S4). The combined measures of SDS at maturity did not provide any MTA at the stringent significant threshold utilized. However, the manhattan plot for SDS at maturity clearly show a significance peak barely missing the threshold at ~4 cM on Chr 4B (Table S3, Figure S5). Although not significant in absolute terms, the shape of this this peak is suggestive of the presence of multiple markers in LD with a genetic element influencing the phenotype. The combined measure of SPC, representing the relative height of the plant reached by the disease, identified two major putative QTL at 71.6–72.5 cM on Chr 3A (qSTB.3) and at 167.5 cM on Chr 4A (qSTB.4), supported by nine and one MTAs, respectively (Table 2, Figure 4, Figure S5). Combined SPC identified also 11 unmapped MTAs that were traced to qSTB.3 (Table S4).

In order to provide a finer dissection of the molecular basis of STB resistance in Ethiopian durum wheat landraces, the GWA scan was repeated on STB data specific to 2012 (Figure S6) and 2013 (Figure S7). Although single-year GWA data allows a lesser degree of generalization, it may account for yearly specificity in STB infections. The SDS at heading in 2012 also identified qSTB.1, reported by the combined measures over the 2 years as well (Table 2, Figure 4). In this case, the putative QTL was supported by six MTAs. SDS at heading in 2012 identified one MTA without genetic position that we mapped at 43.4 cM on Chr 5B (Table S4). The same phenotype measured in 2013 did not provide associations significant at the Bonferroni threshold. SDS measured at maturity in 2012 and 2013 and independently used in a GWA scan did not identify any significant MTAs. The SPC measured in 2012 provided a strong MTA supporting qSTB.3, also identified by the combined measures of SPC (Table 2, Figure 4). The MTAs without map positions also resulted overlapping qSTB.3 (Table S4). SPC data collected in 2013 identified one putative QTL at 135.2 cM on Chr 5A (qSTB.5) (Table 2, Figure 4).

Significant associations deriving from agronomic traits reported R² as high as 26.6% in the case of days to maturity, and above 26% in the case of several markers associated to 1,000 grain weight and number of spikelets per spike (Table S5). Days to heading provided three strong signals on Chr 1B at 81 cM, Chr 2B at 80 cM, and Chr 4B between 3 and 6 cM (Table S7, Figure S8). Days to maturity identified two putative QTL on Chr 1B at 35 cM and on Chr 4A at 51 cM, respectively. A highly significant putative QTL was identified from the GWA scan for plant height on Chr 3B at 166 cM. Spike traits are contributed by several major putative QTL (Table S7, Figure S9). Seeds per spike identified significant MTAs on Chr 3A at 141 cM, on Chr 4A at 170 cM, Chr 4B at 52 and 87 cM, Chr 5A at 14 cM and 142 cM, and on Chr 5B at 43 cM. Spike length identified a unique MTA on Chr 4B at 91 cM. Four putative QTL were identified by measures of 1,000 grain weight on Chr 1A between 60 and 70 cM, on Chr 5A at 44 cM, on Chr 5B between 38 and 41 cM, and on Chr 6A at 85 cM. The number of spikelets per spike provided suggestive peaks on Chr 1A and 4A, but no highly significant putative QTL.

In few cases, the putative QTL provided by disease data overlap those deriving from agronomic measures. The putative QTL for seeds per spike on Chr 4A is proximate to qSTB.4, identified by combined measures of SPC (Table 2), and that on Chr 5B co-maps with the previously unmapped MTA for SDS at heading in 2012 (Table S4). The putative QTL mapped on Chr 1A for 1,000 grain weight co-maps with qSTB.1, and that on Chr 5B overlaps that for seeds per spike and the previously unmapped MTA for SDS at heading in 2012 (Table S4). Putative QTL from phenology and plant height do not co-map with any STB MTA.

Several unmapped MTAs for agronomic data were traced to their most plausible genetic position, many of which overlapped putative QTL identified by mapped markers (Table S4). Novel genomic loci of agronomic significance reported by MTAs with estimated map positions include several putative QTL for seeds per spike, at ~108 cM on Chr 2A, 162 cM on Chr 4A, 0, 18 and 60 cM on Chr 4B, and 95 cM on Chr 7A. An MTA mapping at ~40 cM on Chr 4B is jointly identified by 1,000 grain weight and seeds per spike.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.