Multi-Location Evaluation of Global Wheat Lines Reveal Multiple QTL for Adult Plant Resistance to Septoria Nodorum Blotch (SNB) Detected in Specific Environments and in Response to Different Isolates

The slow rate of genetic gain for improving resistance to Septoria nodorum blotch (SNB) is due to the inherent complex interactions between host, isolates, and environments. Breeding for improved SNB resistance requires evaluation and selection of wheat genotypes consistently expressing low SNB response in different target production environments. The study focused on evaluating 232 genotypes from global origins for resistance to SNB in the flag leaf expressed in different Western Australian environments. The aim was to identify resistant donor germplasm against historical and contemporary pathogen isolates and enhance our knowledge of the genetic basis of genotype-by-environment interactions for SNB response. Australian wheat varieties, inbred lines from Centro Internacional de Mejoramiento de Maiz y Trigo (CIMMYT), and International Center for Agricultural Research in the Dry Areas (ICARDA), and landraces from discrete regions of the world showed low to moderate phenotypic correlation for disease response amongst genotypes when evaluated with historical and contemporary isolates at two locations across 3 years in Western Australia (WA). Significant (P < 0.001) genotype-by-environment interactions were detected regardless of same or different isolates used as an inoculum source. Joint regression analysis identified 19 genotypes that consistently expressed low disease severity under infection with different isolates in multi-locations. The CIMMYT inbred lines, 30ZJN09 and ZJN12 Qno25, were particularly pertinent as they had low SNB response and highest trait stability at two locations across 3 years. Genome wide association studies detected 20 QTL associated with SNB resistance on chromosomes 1A, 1B, 4B, 5A, 5B, 6A, 7A, 7B, and 7D. QTL on chromosomes 1B and 5B were previously reported in similar genomic regions. Multiple QTL were identified on 1B, 5B, 6A, and 5A and detected in response to SNB infection against different isolates and specific environments. Known SnTox-Snn interactions were either not evident or variable across WA environments and SNB response may involve other multiple complex biological mechanisms.


INTRODUCTION
Septoria nodorum blotch (SNB) caused by the fungal pathogen Parastagonospora (syn. ana, Stagonospora; teleo, Phaeosphaeria) nodorum (Berk.) Quaedvlieg, Verkley, and Crous is a significant necrotrophic disease on the leaf and glume of bread wheat. Disease epidemics cause significant yield losses in wheat growing regions of the world, including annual yield losses of 12.9% in Western Australia (WA; Murray and Brennan, 2009). SNB resistant varieties through breeding efforts are estimated to contribute 30% of pathogen control whereas cultural practices and fungicide application contributes the remaining 70% in WA (Murray and Brennan, 2009). Therefore, deploying new improved SNB resistant varieties whilst integrating on-farm management practices as a dual approach provide significant benefits to control the pathogen and disease.
Parastagonospora nodorum have asexual pycnidiospores and sexual ascospores where the latter is the primary source of inoculum (McDonald et al., 1994). There is a high degree of genetic diversity within and between isolate populations from major wheat growing continents with annual cycles of sexual reproduction contributing to population structure (Stukenbrock et al., 2006). Allele assortment, therefore, significantly influences isolate aggressiveness on the host plant in specific environments (Engle et al., 2006;Ali and Adhikari, 2008) whereby breeding for adult plant resistance (APR) would benefit from selection via natural infection or using a mixture of the most genetically diverse isolates for artificial inoculation in a particular environment (Cowger and Silva-Rojas, 2006). It is unlikely that any single isolate or a specific environment would provide sufficiently robust evaluation for host APR in the field as pathogen aggressiveness amongst isolates and host responses vary between wheat genotypes (Scharen and Eyal, 1983;Scharen et al., 1985;Eyal, 1999;Ali and Adhikari, 2008). Consequently, selection of isolates from diverse geographic origins and their effect on wheat cultivars in different environments are important considerations when evaluating disease symptoms and understanding the genetic control of resistance to SNB. Moreover, variation in heading date and height are important determinants in breeding and selection whereby wheat lines ranked as SNB resistant variation could be misclassified particularly for late maturing genotypes that escaped disease infection (Francki, 2013). In some cases, genetic factors controlling height and heading date were, indeed, linked to genes controlling SNB response rather than pleiotropic effects of agronomic characteristics affecting disease evaluation (Shankar et al., 2008;Francki et al., 2011).
Genetic analysis has shown seedling resistance to be under the control of minor genes and independent to APR genes when a bi-parental mapping population was evaluated for SNB response specifically in WA environments (Shankar et al., 2008) with general disease susceptibility often increasing toward physiological maturity of the plant (Develey-Rivière and Galiana, 2007). Infection at or beyond ear emergence is the most damaging, reducing photosynthetic capacity and grain yield (Mullaney et al., 1983;Spadafora et al., 1987;Wainshilbaum and Lipps, 1991). Given disease symptoms proliferate in the warmer spring temperatures in WA as the crop matures, APR is therefore, a high priority in breeding to maintain photosynthetic capability in the leaf and glume when SNB epidemics are at the greatest inseason risk. APR is under polygenic control and genes for leaf and glume resistance have previously been reported to be minor with additive effects (Rosielle and Brown, 1980;Fried and Meister, 1987;Bostwick et al., 1993). A number of genetic studies based on bi-parental mapping populations identified quantitative trait loci (QTL) specifically for SNB resistance on the flag leaf and glume on chromosomes 1B, 3B, 2A, 2D, 4B, and 5B using either mixed isolates as inoculum or natural infection in field environments (reviewed in Francki, 2013). In some instances, QTL for SNB resistance were not detected in all environments (Shankar et al., 2008;Francki et al., 2011;Lu and Lillemo, 2014) indicating that host resistance is controlled by genes expressed in response to different isolates, environments or isolate-by-environment interactions. Moreover, a number of the QTL identified using moderate density molecular marker maps appeared to be colocated on the same regions, such as those on chromosomes 1B, 2A, 2D, and 5B and genes at these loci controlled SNB resistance for different isolates and environments (Shankar et al., 2008;Francki et al., 2011). However, higher resolution genetic maps using the iSelect Infinium 90K single nucleotide polymorphic (SNP) genotyping array (Wang et al., 2014) allowed for and enabled the discrimination of co-locating QTL for SNB resistance into discrete but linked loci (Francki et al., 2018). Multiple loci on the same chromosome, therefore, responded independently to different isolates and environments.
Although genetic maps with high density SNP markers have been beneficial in discriminating genetic loci in QTL mapping, trait variation in bi-parental mapping populations are constrained to alleles from either parent. Genome wide association studies (GWAS) enhances resolution in markertrait associations whilst simultaneously evaluating a broader gene pool and alleles for desired trait variation. Alleles and marker polymorphisms in a GWAS panel are associated with traits of interest through historical recombination events and linkage disequilibrium (LD). There are increasing reports on the application of GWAS for genetic analysis of traits, including disease resistance. For example, marker-trait associations linked to resistance against stem and stripe rust (Bajgain et al., 2015;Bulli et al., 2016;Gao et al., 2016;Maccaferri et al., 2016) and Fusarium head blight (Arruda et al., 2016) by GWAS further unraveled the genetic complexity of fungal resistance and identified known and discovered novel loci with an expanded repertoire of molecular markers linked to genes controlling disease resistance. Therefore, the application of GWAS allows the evaluation of wider accessions for disease response whilst simultaneously identifying new marker-trait associations with finer mapping resolution based on LD and discovery of new alleles not necessarily represented in parents of bi-parental mapping populations. There are reports for identifying the genetic control of SNB in wheat using GWAS. Included are those for analyzing seedling resistance in a wheat panel of 528 lines that were screened using a single P. nodorum isolate in a single controlled environment using lower density DArT and the Illumina iSelect beadchip 9K SNP assay (Adhikari et al., 2011;Gurung et al., 2014). Considering plant response to disease changes through physiological development from seedling to maturity (Develey-Rivière and Galiana, 2007) and environmental influences on fungal isolate aggressiveness (Pariaud et al., 2009), GWAS analysis using the iSelect Infinium 90K genotyping array (Wang et al., 2014) provides opportunities for a detailed investigation on discrete QTL controlling adult plant response to SNB in multi-environments based on evaluation of a wider gene pool. In recent GWAS studies a number of QTL were detected in only one field environment with some detected in two or more environments and only one QTL on chromosome 2D and 2A detected in all environments (Ruud et al., 2019;Lin et al., 2020). It appears, therefore, that host response to SNB is largely variable and complex where several loci are detected in specific environments indicating considerable genotype-environment interactions.
The aim of this study was to evaluate SNB response of flag leaf in 232 diverse bread wheat accessions from global breeding programs and discover the genetic interaction of host resistance to different P. nodorum isolates in SNB trials by GWAS using the iSelect Infinium 90K genotyping array. The evaluation for SNB response was done in two locations annually for 3 years by inoculating field trials with a mixture of historical and contemporary P. nodorum isolates collected from diverse geographical regions in WA and the potential role of NE-host interactions elucidated by comparing QTL for SNB response with known Tsn and Snn loci. The outcome of the study will refine our knowledge on host genes responding to SNB and the interplay between isolate diversity and different environments to improve breeding resistance in bread wheat.

Plant Material
The GWAS accession panel consisted of 71 wheat lines from Australian origin, 72 inbred and commercial lines from Centro Internacional de Mejoramiento de Maiz y Trigo (CIMMYT), 78 inbred lines from International Center for Agricultural Research in the Dry Areas (ICARDA), and 11 landraces from various origins. Inbred lines were accessed through the CIMMYT Australia ICARDA Germplasm Evaluation (CAIGE) project 1 . 1 http://www.caigeproject.org.au/ Pedigrees and origins of 232 wheat accessions are provided in Supplementary Table S1.

Isolates and Preparation of Inoculum
Parastagonospora nodorum isolates from various wheat growing regions of WA were sourced from the culture collection at the Department of Primary Industries and Regional Development (DPIRD), formerly the Department of Agriculture and Food WA. Inoculum used for evaluation of SNB in field trials each year consisted of a mixture of 16-20 isolates where at least 35% of the inoculum represented contemporary isolates collected from the previous year of the trials. Mixed inoculum used in each year included at least four of the same isolates represented in successive years. Supplementary Table S2 provides details of the isolates used in inoculum each year, including the origin and year of collection. Fungal cultures were prepared by growing the isolates on sterile wheat grain for 3 months to induce formation of pycnidia and pycnidiospores (Fried, 1989). Fungal cultures were air-dried and ground to a coarse powder, equal amounts of each isolate was mixed and stored at 4 • C. The mixed fungal isolates were rehydrated to a suspension of 10 6 spores/ml with 0.5% gelatine prior to use as inoculum in field trials.

Field Trials
Trials in 2016 and 2017 were sown at Northam and Katanning, and Northam and Manjimup in 2018. Locations were selected based on varying mean annual temperatures, rainfall and maximum relative humidity (Northam: 18.1 • C, 432 mm, 71%; Katanning: 15.4 • C, 470 mm, 70%; and Manjimup: 14.7 • C, 981 mm, 76%). Trials were arranged in completely randomized block design with 3 replications for each accession with 0.6 m spacing between plots. Each plot consisted of a single 1.9 m length row with 0.4 m row spacing in 2016 and 2017 trials and two rows of 1.9 m length with 0.2 m row spacing in 2018. A spreader 2-row plot of the susceptible early maturing variety "Amery" was adjacent to each treatment plot. Susceptible check cultivars (6 reps) comprising cultivars "Amery" (early maturing), "Arrino" (early-mid maturity), "Millewa" (mid-late maturity), "EGA2248" (mid-late maturity), and "Scout" (late maturity) were included in all trials. Trials were rainfed or irrigated with above ground sprinklers to simulate a rainfall event and induce disease when required.
Inoculum (10 6 spores/ml with 0.5% gelatine) was applied using a motorized mister commencing at Feekes growth stage 3-4 with a further three inoculations at 10-14 day intervals. Inoculum was applied at a rate of 28.5 m 2 /L prior to a rainfall event. Sprinkler irrigation was applied to trials in the absence of rainfall within 2 days post inoculation. SNB response on the flag leaf was visually assessed using percent leaf area diseased (PLAD) scale on susceptible check cultivars where 0% and 100% were scored as highly resistant and highly susceptible, respectively. All accessions were scored in each trial when at least 2 susceptible check varieties had PLAD scores >70%. Disease ratings were scored for 5 random individual plants from the middle of each plot and mean plot scores were used in statistical analysis.
Measurements for plant height (cm) were taken as the distance from the top of the soil to the top of the heads (excluding awns). Three measurements were taken in each plot and mean plot scores used for analysis. Heading date was measured as days from date of sowing to when 50% of the ears were fully emerged for each plot.

Analysis of SNB Disease Response
All statistical analyses for phenotypic evaluation were done using Genstat, 19th edition 2 . Disease and agronomic data were tested for normality and error variance homogeneity across environments. Residuals plotted against fitted values revealed a random distribution (data not shown) indicating there was no need for data transformation. Generalized linear models and linear mixed models were used in phenotypic analysis of trait data. Treatment factors and co-variates were fitted to fixed models to estimate main effects and interactions. Finlay-Wilkinson joint regression analysis was used to compare genotypes for SNB response and agronomic traits at two locations across 3 years. Broad-sense heritability estimates were calculated using the formula H 2 = σ 2 g /σ 2 g + σ 2 e /r, where σ 2 g , and σ 2 e are the genotypic and error variance, respectively, and r is the number of replications.
Genotyping DNA samples from wheat accessions were assayed using the 90K Infinium SNP chip array. Raw intensity data was analyzed in GenomeStudio, and NormTheta and NormR values for each SNP were extracted. A custom perl script was used to cluster samples and assign genotypes to known polymorphisms. Homozygous genotypes were called when a sample was located at a previously identified cluster from a bi-parental mapping population, whether they had been genetically mapped or not. Clusters that had been genetically mapped in bi-parental mapping populations were assigned to chromosomes and clusters reported in the wheat SNP Consensus map 90K Array (Wang et al., 2014) were also assigned a map position, with some SNP having multiple loci. Physical location of genetic markers were extracted from the International Wheat Genome Sequencing Consortium (IWGSC) RefSeq v1.0 physical map via https://wheat-urgi.versailles.inra.fr/ (Alaux et al., 2018).
The Tsn1 locus was genotyped across lines of the GWAS panel using the linked marker fcp620 as previously described (Zhang et al., 2009). The effect of sensitive versus insensitive alleles was then determined for each environment (p-value threshold p < 0.05) to establish whether Tsn1 had a significant effect on PLAD scores. The genotypes for fcp620 were also compared to the marker-trait associations (MTA) co-locating with Tsn1, IWB14942, and IWB43679.

Population Structure, Genome-Wide Association and Statistical Analysis
Monomorphic markers and markers with less than 80% call rate were removed from the dataset. Markers with minor allele 2 https://genstat.kb.vsni.co.uk frequencies (MAF) less than 5% were removed prior to analysis, reducing the total number of SNPs to 20,563. The data was then split into three separate analyses, one set containing the filtered data set of 20,563 SNPs, a set where the markers were further filtered to include SNPs with no less than 90% call rate and pruning the data based on a LD of r 2 ≤ 0.2 using the PLINK 1.9 command "-indep-pairwise 50 2 0.2, " reducing the data set to 2,941 SNPs, and a third set where the markers were filtered to include SNPs with no less than 90% call rate and pruning the data based on a LD of r 2 ≤ 0.1 using the PLINK 1.9 command "-indep-pairwise 50 2 0.1, " reducing the data set to 1,142 SNPs.
To generate a Q matrix the LD pruned data (r 2 ≤ 0.1) was imported into STRUCTURE version 2.3.4. (Pritchard et al., 2000;Falush et al., 2003Falush et al., , 2007Hubisz et al., 2009) and applied to full set, pruned r 2 ≤ 0.2 and pruned r 2 ≤ 0.1 marker datasets for MTA. A burn-in/MCMC of 20,000/40,000 was applied, selecting the Admixture model with correlated allele frequencies from K = 1-10 with 10 independent runs each. The output was imported into Structure Harvester (Earl and von Holdt, 2012) and determined that K = 2. CLUMPP version 1.1.2 (Jakobsson and Rosenberg, 2007) was implemented to summarize the output as a Q matrix for use in association analysis in TASSEL v.5.2.52 (Bradbury et al., 2007). Within TASSEL a genotypic kinship matrix (K) was estimated by selecting the "Centered_IBS" method. General linear model (GLM; Q), GLM (PCA), mixed linear model (MLM; Q + K), and MLM (PCA + K) were all initially explored with visual assessment of quantilequantile plots (Q-Q plots). The suitable number of PCs for each trait was determined by testing one through 15 PCs. The option "P3D" was not selected during the MLM analysis with the variance component re-estimated after each marker. After assessing Q-Q plots it was determined that MLM (PCA + K) was the most suitable method, accounting for both population structure and cryptic relatedness. Significant associations were only detected in the full data set with no significant associations detected in either pruned data sets. The pruned data sets were removed from further analysis. The R programs "qqman" and "Rcolorbrewer" were used to draw Manhattan plots (Turner, 2017;R Core Team, 2018). Two-dimensional displays of the top PCs were drawn in R.
A genome-wide significance threshold for MTAs was set at p < 2.43 × 10 −6 [−log 10 (p) > 5.61] using Bonferroni correction with α = 0.05. Bonferroni correction is highly conservative and reduces type I errors, however, is most practical when all tests are independent (Bland and Altman, 1995;Abdi, 2007). To estimate the number of independent tests the tagger function in Haploview was implemented as described in Maccaferri et al. (2016) with a r 2 of 0.1. This returned a genome-wide threshold significance of p < 7.65 × 10 −5 [−log 10 (p) > 4.12] and was considered a moderate level of significance compared to the Bonferroni corrected threshold. As reported in a number of related studies a threshold of significance of p < 1 × 10 −3 [−log 10 (p) > 3.00] was included as a suggestive level of significance (Maccaferri et al., 2015;Alomari et al., 2017;Muqaddasi et al., 2019).
Linkage disequilibrium was assessed to identify MTAs that occurred as single entities or clusters of SNPs in strong LD. The "-block" function in PLINK 1.9 was utilized to identify intrachromosomal SNPs in strong LD as previously defined (Gabriel et al., 2002) with the bottom of the 90% D-prime confidence interval greater than 0.70 and the top of the confidence interval at least 0.98 (Purcell et al., 2007). LD decay was estimated to gage the approximate size of QTL intervals. Marker pairwise r 2 values were calculated in PLINK 1.9 with a sliding window of 50 and then LD decay curves fitted by non-linear regression for each subgenome (A, B, and D) as previously described (Marroni et al., 2011) with decay of r 2 against distance. LD decay plots were drawn in R with a critical threshold of r 2 = 0.2 which represented the LD half decay point (R Core Team, 2018). QTL were defined as having moderate to highly significant MTA for a single SNP locus. Linkage decay values were used as estimates for a QTL interval when multiple moderate to highly significant MTA were identified within a similar genomic region.
The most significant SNP marker was chosen as a representative for all QTL and their individual and collective effect of allele stacking on PLAD scores for each environment was assessed and P-values calculated in R (R Core Team, 2018).

Phenotypic Analysis of SNB Response
The population mean values for SNB response measured as PLAD ranged from 25.0 to 53.0 ( Table 1) indicating variable disease pressure across environments and years. Analysis of variance indicated significant differences between genotypes for PLAD, heading date and height in each environment ( Table 1). High broad sense heritability was observed for PLAD in each of the two locations across three years (H 2 = 0.64-0.88) and similarly for heading date and plant height ( Table 1).
Pearson's correlation co-efficient for PLAD scores within and between years was low to moderate but highly significant (r = 0.272 to 0.568, P < 0.001; Table 2). It appears, therefore, that SNB response of many individual genotypes is variable even when inoculated with the same isolates across different locations in the same year indicating host-by-isolate-byenvironment interactions. The exception was between Northam and Manjimup in 2018 where correlation for PLAD was high (r = 0.797, P < 0.001; Table 2) indicating that, in some instances, the relative response of genotypes to SNB between environments in the same year was similar when inoculated with the same isolates.
The contribution of genotype, environment and their interactions in the same year was further analyzed by fitting genotype, environment and their interactions as terms in linear mixed models. Among the three sources of variation, the largest proportion of SNB response was significantly (P < 0.01) accounted by genotypes (72.94-79.71%) followed by genotype-by-environment interactions (10.22-22.63%) and environment (4.43-10.07%) for each year (Table 3). Moreover, a highly significant (P < 0.001) but smaller proportion of variation was accounted by genotypes at two locations across 3 years (44.52%) whereby genotype-by-environment interactions (34.52%), and environment (20.96%) was higher compared to  their respective sources of variation across environments within years ( Table 3). The genotype-by-environment interactions within and between years indicated that either varying host response to disease, genetically distinct isolates, diverse isolate virulence, or aggressiveness or a combination of these factors interact to significantly influence the expression of SNB response across genotypes in different environments.

Genotype Performance Using Joint Regression Analysis
Despite highly significant genotype-by-environmental interactions for PLAD within and between years ( Table 3), means of each genotype were fitted with average environmental means using a Finlay and Wilkinson joint regression model to identify genotypes with low mean PLAD scores at two locations across 3 years and ranked in ascending order based on sensitivity to SNB response ( Table 4). Genotypes with PLAD scores <30.0 and with similar heading date and plant height to the susceptible check varieties were determined as consistently expressing SNB resistance. A total of 19 genotypes had mean PLAD values across all locations ranging from 15.33 to 28.94 compared to control susceptible varieties with similar heading dates and plant height ( Table 4). The CIMMYT inbred lines ZJN12 Qno 25 and 30ZJN09, in particular, were identified as having low SNB severity with high stability and lowest mean square deviation indicating greater predictable SNB disease response in any given environment (Table 4). Interestingly, ZJN12 Qno 25, and 30ZJN09 had comparable or improved resistance to SNB but with higher stability and predictability than other resistant lines, 6HRWSN125, and EGA Blanco ( Table 4) that have been deployed as donor parents in doubled haploid populations for QTL studies (Shankar et al., 2008;Francki et al., 2011). Moreover, ZJN12 Qno 25, 30ZJN09, 6HRWSN125, and EGA Blanco have distinctly different pedigrees (Supplementary Table S1) so it is likely that the source of SNB resistance is from different parental origins.

Fixed Effects of Morphological Traits on SNB Response
Pleiotropic effects of morphological characteristics can have significant implications when interpreting the genetic control of SNB response and often exacerbated in QTL studies when populations have significant differences in trait measurements. There was a significant and moderate to high negative correlation between heading date and PLAD scores (r = -0.51 to -0.86, P < 0.001) with similar negative correlation between plant height and PLAD scores (r = -0.52 to -0.79, P < 0.001) in all environments supporting potential pleiotropic effects of morphological characteristics on disease evaluation. Since heading date and plant height were significantly different (P < 0.001) between genotypes in each environment (Table 1), they were fitted as variates in a linear mixed model to estimate the significance of any main fixed effects and their interactions on PLAD. Heading date had highly significant (P < 0.001) main effects on SNB response at two locations across 3 years when adding to or sequentially dropping terms from the fixed model whereas plant height had significant main effects (P < 0.05) in most environments ( Table 5). There was no significant heading date-by-height interactions (P > 0.05) in any of the environments ( Table 5). Heading date and plant height were, therefore, fitted as co-variates in general linear model and adjusted mean PLAD scores were used for subsequent genetic analysis to reduce spurious MTA with agronomic characteristics and improve the accuracy of association of SNP markers with genes controlling PLAD response in GWAS. a Genotypes ranked in ascending order from the most stable assessed by sensitivity of PLAD response to environmental effects. b Predictability of response to SNB assessed by mean square deviation.

Population Structure and Linkage Disequilibrium
A total of 20,563 markers filtered from the 90K Infinium SNP chip array was used to investigate the relatedness of the GWAS panel based on principal component analysis (PCA). PCA described 15.6% of the genetic variance between PC 1, 2, and 3 (6.2, 4.7, and 4.7%, respectively) and regarded as having low population structure (Figure 1). Despite wheat lines sourced from different breeding programs and continents and selected based on different pedigrees (Supplementary Table S1), demarcation for clustering into distinct sub-populations was not apparent (Figure 1) indicative of shared genetic relatedness of distant ancestors in inbred lines and varieties. A comparison between the entire 20,563 SNP marker dataset and the LD pruned datasets (2,491 SNP with no less than 90% call rate and LD of r 2 ≤ 0.2 SNP, and 1,142 SNP with no less than 90% call rate and LD of r 2 ≤ 0.1) showed no difference in principal component analysis and population structure, with the top PCs accounting for slightly less of the total variation in the LD pruned datasets (Figure 1 and Supplementary Figure S1). The extent of LD of the GWAS panel was estimated based on the pairwise squared correlation coefficients (R 2 ) for the 20,563 SNP loci. LD decay was estimated by non-linear regression at 6.97 cM and 6.14 cM for the A and B subgenomes, respectively, whereas the D subgenome had a significantly higher LD of 24.62 cM for threshold R 2 = 0.2 (Figure 2).

Marker-Trait Associations for SNB Response
A mixed linear model (MLM) was applied to account for any confounding effects of population structure and cryptic relatedness and reduce rate of false positives or spurious SNP marker associations linked to genes controlling SNB resistance. Adjusted mean PLAD values from six individual environments in 2016-2018 were used to reduce any confounding effects of heading date and height and analyzed using 20,563 filtered SNP markers. Analysis of quantile-quantile (Q-Q) plots for adjusted mean PLAD scores showed deviations of the observed association compared to the association statistics expected under the null hypothesis of no association (Figure 3), indicating SNP markers were associated with trait variation at two locations across 3 years. Association tests were analyzed independently for each of the environments with thresholds of p < 2.43 × 10 −6 [−log 10 (p) > 5.61], p < 7.65 × 10 −5 [−log 10 (p) > 4.12], and p < 1 × 10 −3 [−log 10 (p) > 3.00] ranking as a high, moderate and suggestive levels of significance, respectively, in Manhattan plots (Figure 3).
A total of 47 SNP loci associated with PLAD with moderate to high significance were detected on nine chromosomes at two locations across 3 years including 1A, 1B, 4B, 5A, 5B, 6A, 7A, 7B, and 7D ( Table 6). None of the SNP markers detected for SNB were associated with known loci controlling variation for heading date and plant height (Supplementary Tables S3, S4 and Supplementary Figures S2, S3) using moderate to high threshold values of p < 7.65 × 10 −5 [−log 10 (p) > 4.12] so it appears that markers are associated with disease response and not a pleiotropic effect of agronomic characteristics. The number of QTL detected for each environment was variable. A minimum of two QTL were detected for Katanning 2017 and Northam in 2018 with a maximum of seven QTL detected for Katanning in 2016 (Figure 3 and Table 6). The majority of MTA were detected on the A and B subgenomes with fewer genes controlling PLAD on the D subgenome when the GWAS panel was evaluated in all environments (Figure 3). Although SNP were evenly spaced across the genome, the average distance between SNP ranged from 0.52 cM in the full dataset to 9.06 cM in the smaller pruned dataset of 1,142 SNP. The full and pruned dataset had large gaps of 79.53 cM and 135.67 cM, respectively, both on the D subgenome. SNP per chromosome ranged from 83 SNP to 1907 SNP for the full dataset and 19 SNP to 111 SNP for the pruned dataset of 1,142 (Supplementary Table S5).
Since genotype-by-environment interaction was significant ( Table 3) and the disease correlation showed low to moderate Pearson's coefficients between environments (Table 2), it was expected that the same QTL would not necessarily be detected across environments within and between years. A total of 20 QTL were identified at two locations across 3 years (  (Figure 3 and Table 6). Interestingly, TABLE 6 | Summary of SNP marker associations with PLAD scores from two locations across 3 years in 2016-2018. Gray shading represent SNP markers that are in strong LD as defined by Gabriel et al. (2002).

Environment
Chromosome    Table 6). Phenotypic variance accounted by each SNP ranged from 7 to 10% and allele effects on reducing disease score ranged from 7.47 to 24.45% (Table 6). Therefore, multiple loci on the same chromosome control SNB resistance with varying effects of genes and alleles on reducing disease severity within and between environments and isolates. Individual SNP alleles representative of each locus consistently showed significantly (P < 0.05) reduced PLAD across all environments when stacking of more than one allele (Supplementary Figure S4). Marker alleles for QTL were different for some SNP in each of the four resistant wheat genotypes (  Table S1) with consistent low PLAD at two locations across 3 years ( Table 4). The resistant wheat genotypes probably inherited different genes and alleles for SNB response. Single nucleotide polymorphic markers from the 90K Infinium SNP chip array used to genotype the GWAS panel enabled a direct comparison with QTL based on the physical map marker position (Alaux et al., 2018) and the wheat consensus map (Wang et al., 2014). Loci detected in this study, therefore, were compared to those previously reported for SNB response in other WA environments on chromosome 1B and 5B from a bi-parental doubled haploid (DH) mapping population derived from the parents EGA Blanco and Millewa (Francki et al., 2011;Francki et al., 2018). All of the 22 MTA identified for four QTL on chromosome 1B from GWAS analysis ( Table 6) did not co-locate with the QTL, QSnl07.daw-1B, and QSnl08.daw-1B (4.3 Mbp to 6.9 Mbp) in the EGA Blanco/Millewa population (Figure 4). Therefore, it appears that regions on chromosome 1B contain multiple genes that respond to SNB in different environments and isolates. Similarly, the SNP marker IWA5670 (at 514.9 Mbp) associated with SNB resistance in Katanning in 2016 was identified flanking the QTL QSnl07.daw-5B in the EGA Blanco/Millewa mapping population (Figure 4) and likely represents an alternative locus on chromosome 5B responding to SNB in specific environments and against diverse isolates. SNP markers 1WB43679 (at 539.5 Mbp) and IWB14942 (at 546.7 Mbp) on chromosome 5B associated with SNB resistance at Northam and Manjimup in 2018, respectively, aligned to the QTL, QSnl07.daw-5B in the EGA Blanco/Millewa population (Figure 4) indicating similar regions harbor QTL detected in more than one environment.

Comparison of MTA for SNB Resistance With Known SnTox-Snn Loci
The distribution of MTA across nine chromosomes were compared to locations for known Snn loci. SNP markers identified for Tsn1 (Ruud et al., 2017;Ruud et al., 2019) Snn1 and Snn3-B1 (Ruud et al., 2017;Downie et al., 2018) loci were assigned directly to the genetic consensus map (Wang et al., 2014). The closest SSR markers linked to Snn4 (Abeysekara et al., 2009), Snn5 (Friesen et al., 2012), and Snn6 (Gao et al., 2015) were crossreferenced with genetic maps consisting of SSR and SNP markers (Francki et al., 2018) to identify closely linked SNP as anchoring markers to assign Snn4, Snn5, and Snn6 on the genetic consensus map (Wang et al., 2014). A total of six from eight known NE-host loci were mapped to nine chromosomes containing QTL detected across WA environments (Figure 5). The Tsn1 locus represented by marker fcp620 identified on 5B in the EGA Blanco/Millewa population was co-located with two SNP markers, IWB14942 and IWB43679, at 252.96 cM on the consensus map (physical positions 546.70 Mbp and 539.46 Mbp, respectively) associated with SNB response in Northam and Manjimup 2018 (Figure 4). A QTL on 6A detected from the Northam 2018 environment only was within 6 cM of Snn6 whereas the co-location of QTL on remaining chromosomes was not apparent with other known Snn loci (Figure 5).
The marker for Tsn1 (fcp620) was used to genotype the GWAS panel for sensitivity or insensitivity to SnToxA with 100% genotype identity with SNP marker IWB14942 located at 546.70 Mbp (Figure 4). However, an additional SNP marker FIGURE 6 | Allelic effects at fcp620 locus at two locations across 3 years (2016)(2017)(2018). P-values for significance are shown.
IWB43679 located at 536.46 Mbp had only 69% genotype identity and, therefore, did not co-segregate with Tsn1. Box plot analysis confirmed the allelic effects of fcp620 linked to Tsn1 significantly (P < 0.05) reducing PLAD in only two environments, Northam 2018 and Manjimup 2018 (Figure 6) which corroborated similar allelic effects of IWB14942 at these sites (Supplementary Figure S5).

DISCUSSION
The analysis of global wheat accessions for their response to APR against mixed isolates of SNB in WA field environments and cultivar-by-isolate-by-environment interactions are poorly understood. As environmental influences may either affect host gene response, virulence or aggressiveness of genetically diverse isolates, this study evaluated wheat accessions in two environments for three successive years against a mixture of isolates to expand our knowledge on SNB disease response and capture genotype-isolate-environment interactions. The outcomes of this study provided further insight to the complexity of these interactions and the underlying genetic control of SNB response for Australian wheat cultivars, inbred lines from CIMMYT, ICARDA, and some landraces.
Parastagonospora nodorum isolates have displayed a high degree of genetic diversity and the ability to evade plant host resistance (McDonald and Linde, 2002). Reports have shown that a wide range of P. nodorum isolates collected from different geographical regions of the United States can differ significantly in their aggressiveness and vary independently of the cultivar tested (Rufty et al., 1981;Scharen et al., 1985;Krupinsky, 1997), so it is plausible that isolates collected in diverse geographical regions of WA were genetically distinct with varying effects on disease development in regions other than the environment from where they were collected. Moreover, a number of environmental factors including climatic conditions for optimal spore germination, lesion development and sporulation or indeed the genetic basis influencing relative fitness and differential adaption of isolates to host cultivars can influence aggressiveness in different environments (Pariaud et al., 2009). SNB evaluation in this study detracted from using a single isolate as inoculum in field trials as it may neither be prominent in any season and environment, a strain assortment representative of a contemporary isolate, nor represent the breadth of diversity needed for identifying effective and robust SNB resistance. Instead, historical and contemporary SNB isolates presumed to be genetically distinct were used in mixed inoculum for evaluation and comparison of SNB response across wheat genotypes in different production environments. Highly significant genotype-by-environment interactions was consistent within and between years in this study and provided evidence that either alternative host genes, isolate virulence or aggressiveness, environment influences or a combination of these factors contribute to varying SNB response of wheat genotypes across environments. We anticipated consistent high correlation across genotypes in each year when inoculated with the same isolates under similar conditions in different trial locations but this was not evident particularly in 2016 and 2017 and across all years. Instead, both low to moderate correlation and high genotype-by-environment interactions were apparent. It is, therefore, reasonable to assume that different environments affect expression of alternative host genes against the same and different isolates. Conditions favoring isolate virulence or aggressiveness in one environment may also be different in another with a profound and variable effect on different host genes. Further studies to identify genetic factors controlling variability in isolate virulence and their interactions with alternative host genes would continue to unravel the biological and genetic complexity of genotypeby-environment-by-isolate effects on varying SNB disease response in wheat.
Despite varying SNB response in genotypes across environments, this study exploited the genetic diversity of breeding programs in Australia, Mexico, Middle East and landraces from discrete regions of the world to evaluate SNB response in WA environments. There were a total of 19 genotypes having robust resistance at two locations across 3 years against a total of 42 different isolates. Interestingly, EGA Blanco and 6HRWSN125, previously used as donor parents in QTL analysis (Shankar et al., 2008;Francki et al., 2011) showed consistent low PLAD scores indicating these lines have sustained expression of resistance particularly against contemporary isolates. Similarly, inbred lines from CIMMYT including 30ZJN09 and ZJN12 Qno 25, exhibited consistent stable low mean PLAD scores against historical and contemporary isolates in each year. Since EGA Blanco, 6HRWSN125, 30ZJN09, and ZJN12 Qno 25 have no common pedigree it appears, therefore, that SNB resistance could be derived from different genes or alleles with consistent SNB resistance in multiple environments. The different SNP alleles for QTL on chromosomes 1B, 4B, 5A, 5B, 7A, and 7D partly supports the inheritance of alternative genes or alleles while shared favorable alleles were found on 1B, 5B, and 6A. The four wheat genotypes with stable SNB resistance across multi-environments evaluated in this study are ideal donor parents for SNB resistance breeding.
A comparison of population structure using the full set of 20,563 SNP markers compared to the LD pruned datasets in this study showed no difference in population structure. The use of LD pruned datasets in genome-wide association studies were assumed to be favorable for avoiding deleterious effects caused by overrepresentation of markers for some regions of the genome or by artificially creating structure in the population caused by LD between markers rather than by true ancestry (supplementary note in Price et al., 2006). However, in practise, this is often not found to be the case and it is usually recommended to include a full set of cleaned marker data to accurately model more distant relationships and complex ancestry (Price et al., 2006;Eu-ahsunthornwattana et al., 2014). The PCA using the full set of 20,563 SNP markers confirmed that population structure was low and 15.6% of the total genetic variances was contributed by the first three principal components. Reduced population structure with a similar low total genetic variance was previously reported where wheat genotypes were predominantly sourced from one breeding program having extensive and reciprocal germplasm exchange with other programs resulting in a genetic bottleneck (Arruda et al., 2016). Similarly, CIMMYT wheat germplasm is widely distributed globally (Ortiz et al., 2008) so it is reasonable to conclude that reduced population structure of the panel in this study was a consequence of historical and frequent exchange of germplasm between breeding programs in Australia, CIMMYT and ICARDA. Given that population structure had low genetic variance, long-range LD blocks would be assumed for subgenomes of wheat. Although LD decay for the A and B subgenomes were similar to those reported for spring wheat populations (Chao et al., 2010;Bajgain et al., 2015), LD for the D subgenome was considerably higher in this study. The long-range LD decay to 50% of its original value at 26 cM in this study was similar to the D subgenome of hexaploid wheat reported to be 22 cM in a GWAS panel of germplasm sourced from CIMMYT and breeding programs in South America (Mora et al., 2015). Longer-range LD may be due to selection forces restricting recombination events required to retain favorable alleles in the D subgenome for broader adaption of CIMMYT, ICARDA, and Australian genotypes when breeding for geographically diverse environments. Longer-range LD in the D subgenome compared with the A and B subgenome of hexaploid wheat have also been attributed to recent introgressions and population bottlenecks (Chao et al., 2010).
Significant genotype-by-environment interactions within and between years supported conclusions that QTL may harbor concomitant disease-related genes that respond differently to isolate infection across environments (Francki et al., 2018). GWAS provided further evidence of multiple host genes on the same and different chromosomes that responded independently to different environments, isolates or both. Moreover, low to moderate correlations for PLAD scores across locations in successive years in this study were consistent to those previously reported for field-based SNB response (Shankar et al., 2008;Francki et al., 2011;Lu and Lillemo, 2014;Ruud et al., 2019) and QTL for SNB resistance detected in one WA environment may not necessarily have a significant effect on SNB response in another (Shankar et al., 2008;Francki et al., 2011). New loci on chromosomes 5A and 6A in this study provided further evidence of alternative genomic regions harbored disease-related genes that responded to specific WA environments and isolates.
Genome wide association studies for APR to SNB in multienvironments in Norway based on the iSelect Infinium 90K genotyping array (Ruud et al., 2019) enabled a tentative comparison of marker-trait associations across global field-based studies. SNP markers on chromosomes 1B (206.1 cM) and 5B (214.95 and 252.96 cM) detected in WA environments using the genetic consensus map location were associated with disease response either in Nordic field infection in 2011 or inoculation with individual Norwegian isolates NOR4 or 201618 (Ruud et al., 2019). Further comparative studies will determine whether these regions share common genes or whether distinctive diseaserelated gene clusters respond differently to pathogen infection across diverse environments as proposed by Francki et al. (2018). SNP markers for the remaining QTL detected in this study, however, did not correspond to other QTL in Nordic environments or isolates (Ruud et al., 2019;Lin et al., 2020) indicating that in some instances, alternative loci probably respond to SNB in different environments and against isolates from diverse global origins.
The co-location of QTL with the position of known SnTox-Snn and Tsn loci provides a means to identify the role of NEhost interactions in response to SNB disease in the field. The marker fcp620 (co-segregated with IWB14942) on chromosome 5B indicated that Tsn1 was associated with SNB resistance but only at Northam and Manjimup in 2018 and not the remaining environments. Therefore, it appears that the SnToxA-Tsn1 interaction is inconsistent across environments in WA and corroborates a similar conclusion drawn from a previous QTL study using a bi-parental mapping population (Francki et al., 2011). Moreover, QTL did not co-locate with other known Snn loci indicating that historical and contemporary SNB isolates from WA may secrete alternative effector proteins or environmental influences may have an effect on known host-NE interactions. The inconsistent association with several known NE host susceptibility loci across environments, the preponderance of environment-specific QTL and minimal common loci across QTL from global field-based studies makes it reasonable to assume that alternative and distinct biological mechanisms are largely influenced either by the environment, isolates or their interactions. Biological mechanisms other than known SnTox-Snn interactions, including a myriad of biochemical and physiological process during infection, penetration and colonization (Bellincampi et al., 2014;Presti et al., 2015) may have a significant bearing on host resistance and susceptibility.
Additional studies are needed to identify other possible genes involved in SNB response but would require extensive field evaluation across wider locations and contemporary isolates to elucidate their role in specific environments. In doing so, it would provide a comprehensive analysis of the role of multiple genes on chromosomal blocks that respond to SNB infection and provide further insights on environmental factors and isolates affecting different biological pathways. In the meantime, the challenge for breeding SNB resistance is to transfer chromosomal blocks from donor resistant parents to increase the probability of an individual gene contributing suitable SNB resistance in any specific environment and against contemporary isolates. EGA Blanco, 6HRWSN125, 30ZJN09, and ZJN12 Qno 25 would be logical donor parents and the SNP markers identified on chromosomes on 1B, 5A, 5B, and 6A in this study will be important for marker-assisted selection and tracking chromosomal blocks harboring resistance genes.

CONCLUSION
The wheat response to SNB infections against historical and contemporary isolates under multi-environment field conditions showed considerable genotype-by-environment-by isolate interactions within and between years. Although phenotypic correlation between wheat genotypes sourced from Australian cultivars, inbred lines from CIMMYT and ICARDA and some landraces was low, four genotypes expressed stable phenotypes with consistently low SNB symptoms at two locations across three years and identified as suitable donor parents for breeding SNB resistant wheat for WA production environments. GWAS identified 20 QTL with the majority detected in only one environment, confirming the complex genetic control for SNB response and the role of different environments and diverse isolates on expression of minor host genes. It appears that a majority of SnTox-Snn interactions are not evident with SnToxA-Tsn1 being variable in different Western Australian environments and SNB response may involve other multiple complex biological mechanisms.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material. Disease response for 232 lines at 6 locations and genotype data for SNP markers associated with SNB resistance are available in Supplementary Table S6.

AUTHOR CONTRIBUTIONS
MF designed the experiments, contributed to data acquisition, analysis and interpretation, and a major contributor to writing the manuscript. EW contributed to data acquisition, analysis and interpretation, and a major contributor to writing the manuscript. CM and WM contributed to trial designs, planting, maintenance, and data acquisition. All authors read and approved the final manuscript.

FUNDING
This study was financially supported by the Grains Research Development Corporation through project DAW00248 and gratefully acknowledged.