Genome-Wide Association Studies in Diverse Spring Wheat Panel for Stripe, Stem, and Leaf Rust Resistance

Among several important wheat foliar diseases, Stripe rust (YR), Leaf rust (LR), and Stem rust (SR) have always been an issue of concern to the farmers and wheat breeders. Evolution of virulent pathotypes of these rusts has posed frequent threats to an epidemic. Pyramiding rust-resistant genes are the most economical and environment-friendly approach in postponing this inevitable threat. To achieve durable long term resistance against the three rusts, an attempt in this study was made searching for novel sources of resistant alleles in a panel of 483 spring wheat genotypes. This is a unique and comprehensive study where evaluation of a diverse panel comprising wheat germplasm from various categories and adapted to different wheat agro-climatic zones was challenged with 18 pathotypes of the three rusts with simultaneous screening in field conditions. The panel was genotyped using 35K SNP array and evaluated for each rust at two locations for two consecutive crop seasons. High heritability estimates of disease response were observed between environments for each rust type. A significant effect of population structure in the panel was visible in the disease response. Using a compressed mixed linear model approach, 25 genomic regions were found associated with resistance for at least two rusts. Out of these, seven were associated with all the three rusts on chromosome groups 1 and 6 along with 2B. For resistance against YR, LR, and SR, there were 16, 18, and 27 QTL (quantitative trait loci) identified respectively, associated at least in two out of four environments. Several of these regions got annotated with resistance associated genes viz. NB-LRR, E3-ubiquitin protein ligase, ABC transporter protein, etc. Alien introgressed (on 1B and 3D) and pleiotropic (on 7D) resistance genes were captured in seedling and adult plant disease responses, respectively. The present study demonstrates the use of genome-wide association for identification of a large number of favorable alleles for leaf, stripe, and stem rust resistance for broadening the genetic base. Quick conversion of these QTL into user-friendly markers will accelerate the deployment of these resistance loci in wheat breeding programs.

Among several important wheat foliar diseases, Stripe rust (YR), Leaf rust (LR), and Stem rust (SR) have always been an issue of concern to the farmers and wheat breeders. Evolution of virulent pathotypes of these rusts has posed frequent threats to an epidemic. Pyramiding rust-resistant genes are the most economical and environmentfriendly approach in postponing this inevitable threat. To achieve durable long term resistance against the three rusts, an attempt in this study was made searching for novel sources of resistant alleles in a panel of 483 spring wheat genotypes. This is a unique and comprehensive study where evaluation of a diverse panel comprising wheat germplasm from various categories and adapted to different wheat agro-climatic zones was challenged with 18 pathotypes of the three rusts with simultaneous screening in field conditions. The panel was genotyped using 35K SNP array and evaluated for each rust at two locations for two consecutive crop seasons. High heritability estimates of disease response were observed between environments for each rust type. A significant effect of population structure in the panel was visible in the disease response. Using a compressed mixed linear model approach, 25 genomic regions were found associated with resistance for at least two rusts. Out of these, seven were associated with all the three rusts on chromosome groups 1 and 6 along with 2B. For resistance against YR, LR, and SR, there were 16, 18, and 27 QTL (quantitative trait loci) identified respectively, associated at least in two out of four environments. Several of these regions got annotated with resistance associated genes viz. NB-LRR, E3-ubiquitin protein ligase, ABC transporter protein, etc. Alien introgressed (on 1B and 3D) and pleiotropic (on 7D) resistance genes were captured in seedling and adult plant disease responses, respectively. The present study demonstrates the use of

INTRODUCTION
Among the many foliar diseases of wheat, rusts are the economically most significant fungal diseases threatening the food security of the world's growing population. There are three types of rusts in wheat, stripe, or yellow rust (YR) caused by the fungus Puccinia stritiformis Westend. f.sp. tritici Eriks. (Pst), leaf, or brown rust (LR) caused by Puccinia triticina Eriks. (Pt), and stem or black rust (SR) caused by Puccinia graminis Pers. f. sp. tritici Eriks. & Henn. (Pgt). More than sixty wheat-producing countries distributed in all continents other than Antarctica have encountered the rusts proving their widespread presence. The damage caused by YR can be as high as 70% in case cultivars are susceptible and the climatic conditions are favorable for an early infection (Chen, 2005). LR has a widespread geographical presence causing considerable yield losses (Marasas et al., 2002). At an early onset, it causes much more damage to crop resulting in yield losses as compared to stem and stripe rusts (Bolton et al., 2008;Huerta-Espino et al., 2011). SR inflicts losses up to USD 1.12 billion worldwide which results particularly due to a reduction in yield and hampered end-use quality of the crop (Pardey et al., 2013). All three rusts pose a major threat to Indian farmers as well. North Indian conditions support the survival and spread of the YR fungal spores. LR is more prevalent in the whole of India as the disease is favored by intermediate temperatures. Also, LR is the most widely distributed amongst rusts and commonly visible in all wheat growing areas during the season. SR is contained mostly in the southern states of the country and survives throughout the year in the Nilgiri hills of southern India Nagarajan and Joshi, 1985).
In disease management practices, planting resistant varieties is the most economical, efficient, and ecologically acceptable tool to manage wheat rusts worldwide (McIntosh et al., 1995;Wiesner-Hanks and Nelson, 2016). Wheat cultivars with diverse resistance are deployed in different areas keeping in mind the pathotype distribution of three Puccinia species on wheat. To date, the reported number of officially designated resistance genes are more than seventy against both YR (78) and LR (77) in addition to about sixty designated resistance genes against SR (McIntosh et al., 2017). Two recently identified Yr genes Yr79 (Feng et al., 2018) and Yr80 (Nsabiyera et al., 2018) have been added to the rust-resistant gene library. Resistance to multiple rusts can be broadly categorized as all stage or seedling resistance (ASR) and adult-plant resistance (APR). Of these reported genes, many genes provide resistance at the seedling stage of plants (McIntosh et al., 2017). Since seedling genes are pathogen race-specific, the cell death phenomenon is activated due to plant hypersensitive response preventing pathogen spread (Ellis et al., 2014;Mondal et al., 2016). This also puts intense selection pressure on the pathogen for its survival. In such cases, pathogen evades and evolves itself which in result renders the deployed all stage resistance gene to be ineffective in a very short time (Line and Qayoum, 1992;Burdon et al., 2014;Li et al., 2014;Niks et al., 2015). The rust resistant wheat with ASR favors the selection of new pathotypes that multiply without any competition on the resistant host resulting in susceptibility of the cultivar/gene. Several all stage resistant Yr genes have become susceptible for instance Yr2, Yr6 to Yr9, Yr17, and Yr 27 . Other genes such as Lr1, Lr13, Lr24, Lr26, Lr37; Sr6, Sr8a, and Sr11 against LR and SR also succumbed (Kolmer et al., 2007;McIntosh et al., 2009McIntosh et al., , 2014Ellis et al., 2014). A well-known outbreak causing great loss of yield was witnessed in the year 1998 due to the new Pgt virulent race Ug99. This resulted in the failure of ASR Sr genes Sr24, Sr31, Sr36, and SrTmp in subsequent years with the emergence of its new variants (Pretorius et al., 2000(Pretorius et al., , 2010Jin et al., 2007Jin et al., , 2008Jin et al., , 2009Visser et al., 2011;Newcomb et al., 2016). On the other hand, APR resistance is usually governed by multiple genes and quantitatively gets less influenced by race-specific pathogens. The involved genes provide non-race-specific partial resistance to all the pathotypes of a given pathogen species, thus making it more durable (Lagudah, 2011;Burdon et al., 2014). Despite the fact that incorporating APR into new cultivars can be difficult when compared to ASR, it was found that many wheat cultivars possessing APR showed durable resistance (Mcintosh, 1992;Boyd, 2005;Navabi et al., 2005;Singh et al., 2005;Ren et al., 2012b;Chen, 2013;Randhawa et al., 2018). Some APR genes when used in combinations have been known to possess durable pleiotropic resistance against multiple wheat rusts and powdery mildew, i.e., Lr34/Yr18/Sr57/Pm38 (on chromosome 7DS), Lr46/Yr29/Sr58/Pm39 (on chromosome 1BL), and Lr67/Yr46/Sr55/Pm46 (on chromosome 4DL) (Lagudah, 2011;Risk et al., 2012;Ellis et al., 2014), of which Lr34 has been studied extensively in different crops including rice, barley, maize, and sorghum (Krattinger et al., 2013(Krattinger et al., , 2019. Deployment of seedling (all-stage resistance), adult plant (non-race specific and race-specific) and slow rusting resistance is to create diversity for resistance (Bhardwaj et al., 2019b). The wheat lines possessing rust resistance at the seedling stage will remain resistant during the whole life of a wheat variety. Therefore, both types of resistance can be amalgamated to strengthen more extended and durable resistance in the future.
As most of the ASR genes go ineffective in the long run and APR genes have only minor effects, more resistance associated genomic regions are needed to be identified and utilized in wheat genetic improvement for rust resistance. Not all of the QTL identified so far were explored for use in marker-assisted selection (MAS). Therefore, finding user-friendly markers becomes even more necessary (Zeng et al., 2019). With the advent of NGS technologies, it became possible to develop high-throughput Single Nucleotide Polymorphism (SNP) markers which are abundant, co-dominant, and present throughout the wheat genome (Allen et al., 2013;Rasheed et al., 2016;Wu et al., 2017Wu et al., , 2018. These markers have advantages over traditional PCR based markers (Chen et al., 1998;Röder et al., 1998;Ren et al., 2012b;Rosewarne et al., 2013;Zhou et al., 2014) for being high-throughput, low cost for genotyping, high efficiency, and allele specificity (Gupta et al., 2008;Colasuonno et al., 2014;Semagn et al., 2014;Long et al., 2017). Genotyping by sequencing (GBS) and array-based SNP genotyping platforms are utilized globally. Several SNP arrays are now available in the market viz. 9K (Cavanagh et al., 2013), 15K (Boeven et al., 2016), 35K (Allen et al., 2017), 50K (Bayer et al., 2017), 90K , 55K, 660K (Jia and Zhao, 2016), and 820K (Winfield et al., 2016). Genome-wide association studies (GWAS) using these chips have certain advantages over bi-parental QTL mapping. Genetic architecture of complex traits in diverse germplasm collections can be studied using GWAS, which detects the genomic regions present in linkage disequilibrium (LD) with genes associated with the trait under study (Hall et al., 2010;Zhao et al., 2011;Riedelsheimer et al., 2012). Due to the accumulation of historical chromosomal recombinations over several generations in a natural population, QTL can be identified using GWAS at a higher mapping resolution (Yu and Buckler, 2006;Semagn et al., 2010). GWAS has been successfully used to study various traits in wheat such as grain yield (Sukumaran et al., 2018), eyespot disease resistance (Zanke et al., 2017), pre-harvest sprouting resistance (Zhou et al., 2017), 36 agro-morphological traits (Sheoran et al., 2019) and so on.
The current study was undertaken with the objective of performing a large-scale association study for identifying genomic regions responsible for resistance to the three rusts from a very diverse set of germplasm comprising of improved genotypes, released varieties, exotic collection, genetic stocks, landraces, and some mutant lines. The study was focused on both seedling and adult plant stage rust response under controlled and natural field conditions, respectively. To identify seedling resistance, the most virulent and predominantly prevalent pathotypes were used to challenge the wheat material in this study. Care was taken to select pathotypes with varied avirulence and virulence structure which would knock down most of the genes in present-day wheat material (Bhardwaj, 2011). Keeping the aforesaid points into consideration seven pathotypes of Pgt (Prasad et al., 2018), five pathotypes of Pst (Gangwar et al., 2019), and six pathotypes of Pt (Bhardwaj et al., 2019a) were selected for screening wheat material at seedling stage under controlled conditions. The optimum conditions ensure foolproof evaluation and selection of rust resistant material. The SNP marker enabled the identification of genomic regions associated with disease resistance will be further firmed up so as to harness less exploited resistant genes for broadening the resistance base thereby avoiding any substantial losses due to rust diseases.

Plant Materials and Genotyping
The diverse panel of 483 genotypes mentioned in our previous report (Kumar et al., 2020) was used for this study. The germplasm comprised of the exotic collection (34), genetic stocks (44), improved genotypes (120), landraces (96), mutant lines (2), and varieties (187). Genotypes in the panel are adapted to different agro-climatic zones of India. Henceforth, this diverse panel will be recalled as rust association mapping panel (RAMP). Genotyping of the panel was done using 35K Axiom R Wheat Breeder's Array (Allen et al., 2017, Affymetrix UK Ltd., United Kingdom) as per manufacturer's guidelines with 35,143 SNP markers. The genotyping data was filtered by retaining markers with minor allele frequency (MAF) > 5%, genotypes with <10% missing SNP calls, and markers with <10% missingness. The obtained markers were further ordered according to the genetic map available at CerealDb Allen et al., 2017). This resulted in a selection of 14,650 polymorphic SNP markers as described previously (Kumar et al., 2020) for subsequent genetic analyses.

Phenotypic Evaluation for Stripe Rust, Stem Rust, and Leaf Rust
The RAMP was evaluated for stripe rust (YR), stem rust (SR), and leaf rust (LR) both at the seedling stage in controlled conditions and at the adult plant stage in the field under natural disease pressure.
Field evaluation of YR was done at two locations for two consecutive years viz. 2017-18 and 2018-19. By counting each location and year as one environment, a total of four environments (E1-E4) were considered. The plants were sown at seed farm, Uchani (29 • 42 48.7 N, 76 • 59 51.3 E), Haryana,  Table S2). The seeds were planted in a non-replicated augmented block design with single row of 1 m and the distance between two rows was 0.3 m. The planting was done in the first fortnight of November at Indore and Wellington and in the second fortnight at Karnal and Uchani, each year. A mixture of check lines susceptible to multiple rusts was planted as infector rows (at every 20th single row) and in spreader rows (perpendicular to the 1 m rows) surrounding the plot for establishing sufficient inoculum and uniform disease development. To ensure uniform disease distribution, spores were collected from the early infections that appeared naturally in the spreader rows and were used to inoculate the infector rows. The response to rust was recorded using disease severity (DS) and infection response (IR) as the two measures. DS was measured using the modified Cobb scale (Peterson et al., 1948) as an estimation of percentage coverage (0,5,10,20,40,60,80, and 100) of rust pustules (uredinia) over the flag leaf. IR was scored as a host reaction to rust pustules and converted to a 0-1 scale (Roelfs et al., 1992). The lines showing the mixed response of moderately resistant to moderately susceptible or vice-versa, were considered as the fifth category other than mentioned in Roelfs et al. (1992). Therefore, five scoring categories considered for the evaluation were: Resistant (R) = 0.2, Moderately Resistant (MR) = 0.4, Mixed response (M) = 0.6, Moderately Susceptible (MS) = 0.8, and Susceptible (S) = 1. Data were recorded at weekly intervals for three times when the flag leaves of the susceptible checks showed a disease score of 60S (DS: 60; IR: S). Out of these multiple scores of a test line, the one with the score tending toward susceptibility was kept for the study. These were further considered for GWAS in each environment by combining the two measures into a single value as coefficient of infection (COI). It is the product of DS and IR on a 0-100 linear scale (Loegering, 1959;Roelfs et al., 1992). Genotypes with COI scores of 0 to 20 were considered as resistant, with score ≥ 60 as susceptible and the remaining as moderate ( Table 1). COI was assumed to be suitable for GWAS and as a primary trait for the identification of significant marker-trait associations (MTAs) since it combines both the information from DS and IR for rust response (Yu et al., 2012;Gao et al., 2016Gao et al., , 2017Mihalyov et al., 2017).

Statistical Analysis
The phenotype data (IT, IR, DS, COI) for the three rusts were visualized and considered for studying Pearson's correlation in R statistical programming. Correlation plots for each rust describing the correlation of disease scores between different pathotypes and environments were created using corrplot R package. Mean comparison tests were performed for IT and COI among population structure based groups using Levene's test (Levene, 1960) dependent two-way t-test at a significance level of P < 0.05 in the R program. The normality of the original data was tested using the Shapiro-Wilk test in IBM SPSS Statistics v.22. A Restricted Maximum Likelihood (REML) approach (Corbeil and Searle, 1976) in R package "lme4" (Bates et al., 2015) was used for estimating variance components for IR and DS from field experiments for the three rusts. A linear mixed model was fitted by considering the overall mean as fixed effect and other factors as random effects. The random effects included genotype (g), location (e), genotype × location interaction (g × e) and year (as replication, r). Following model was applied to estimate variance components: y pqr = µ + g p + r q + e q + ge pq + e pq Where, y pq is the observation for the genotype p at location q in season r, µ is the overall mean, g p the effect of the genotype p, e q the effect of location q, r q the effect of year (season) r, ge pq the interaction between accession p within location q, and e pq the residual.
Broad sense heritability (H 2 ) was estimated by using the following equation: Where, σ 2 g is the genotypic variance, σ 2 e is the environmental variance, σ 2 g×e is the genotype by environment interaction variance, and σ 2 error is the residual error variance and t is the number of years for each location or location by year for the estimates of heritability across all environments. To estimate the stability of genotypes in response to YR, LR, and SR infection in their respective four environments, COI scores were used to perform GGE biplot (genotype × environment interaction) in R package GGEBiplotGUI (Frutos et al., 2014). The COI scores were reversed with reference to maximum score so that the resistant genotypes can have higher scores aiding the interpretation of stable resistant genotypes across environments.

Linkage Disequilibrium (LD), Population Structure, and Genetic Diversity
For the complete RAMP, LD analysis was performed across A, B, and D genomes separately. Intra-chromosomal pairwise marker LD as squared allele-frequency correlations (r 2 ) values were calculated in TASSEL v5.2 (Bradbury et al., 2007) using a sliding window approach with default parameters. As a function of genetic distance, the estimated r 2 − values for significant SNP marker pairs were plotted to understand the extent of LD. A second-degree "loess" function (Cleveland, 1981) in the R statistical program was fitted to estimate the rate of LD decay over genetic distance (cM). Critical r 2 − values were estimated as the 95 th percentile of square root transformed r 2 − values of the unlinked SNP marker pairs (Breseghello and Sorrells, 2006) showing a distance of more than 50 cM for each genome. It indicates the point beyond which the LD is caused by genetic linkage (Breseghello and Sorrells, 2006). The genetic distance at which LD fell below the critical r 2 − value was considered as the confidence interval (CI) of quantitative trait loci (QTL) in this study. In other words, the point of intersection between LD decay loess curve and critical r 2 was used to adjust the QTL-CI in terms of genetic distance.
The population structure and genetic diversity for the RAMP were previously described by Kumar et al. (2020). Briefly, the population structure was estimated with LD-based pruned markers using a Bayesian based model. The optimal number of subpopulations (i.e., K = 2) was determined using parallel runs in STRUCTURE 2.3.4 (Pritchard et al., 2000) with the Linux based python program "StrAuto" (Chhatre and Emerson, 2017). In addition, principal component analysis was also performed on the study panel using 14,650 SNP markers using 'prcomp' function in the R statistical programming language.

Marker-Trait Association Analysis
The seedling response under greenhouse condition and field evaluation for adult plant response were considered for finding out associations. In order to identify the loci associated with the response, 14,650 SNP markers and phenotypic trait values for seedling (IT) and adult plant (COI) responses were used to conduct genome-wide association analyses. For the estimation of marker-trait associations (MTAs), Genomic Association and Prediction Integrated Tool (GAPIT) (Lipka et al., 2012) was used by implementing compressed mixed linear model (CMLM) Zhang et al., 2010) in R environment. A markerbased VanRaden kinship (K) matrix (VanRaden, 2008) for the 483 accessions was also generated using GAPIT. The K and Q (Kumar et al., 2020) matrices were considered as random and fixed components, respectively, to avoid any spurious associations caused by population structure. Some additional association testing models were analyzed to compare the observed probability deviations from the expected distribution based on the Q-Q plots. A general linear model with kinship only and no correction for population structure (GLM_K), GLM (GLM_PC3.K), and mixed linear model (MLM_PC3.K) with kinship and the first three principal components, MLM with kinship and correction for population structure (MLM_Q.K), and MLM with kinship only (MLM_K) were considered. The comparison was done based on P-values and respective Q-Q plots for the MTAs obtained from GAPIT. All the comparative Q-Q plots and circular Manhattan plots were generated using the CMplot R package 1 .
GWAS was conducted on each rust pathotypes separately to identify race-specific seedling resistance loci. In the case of identifying resistance loci in the adult plant stage, COI was considered for each rust in four environments. Each environment (Supplementary Table S2) was considered as a unique set of phenotypic data to be considered for GWAS. The percentage of phenotypic variation explained by the MTA (R 2 ) was calculated as the difference of R 2 without SNP from R 2 with SNP of the GAPIT model (Abed and Belzile, 2019;Garcia et al., 2019). High stringency was observed in GAPIT for false discovery rate (FDR) adjusted p-values. SNP markers were declared significantly associated at p ≤ 0.001 {-log 10 (p) ≥ 3} in the selected model. Hence, a liberal approach similar to previous studies (Pasam et al., 2012;Zegeye et al., 2014;Gao et al., 2016;Visioni et al., 2018), was considered so as to reduce the chance of neglecting any significantly associated marker annotating for disease resistance. A putative QTL was designated for intra-chromosomal SNPs of the MTAs detected and fall within the range for QTL-CI defined by LD. A representative SNP for such putative QTL was selected with the lowest p-value as the most significantly associated SNP.
Significant markers observed were further subjected to in silico annotation. The flanking sequence of these markers was obtained from the EnsemblPlants database spanning 1 kb both upstream and downstream of the SNP position 2 . In order to obtain the reference physical map positions of these markers, the flanking sequence was used to make a query against IWGSC RefSeq v1.0 (Appels et al., 2018). The flanking region of the significant SNP was then used to explore candidate gene overlapping with this region by using JBrowse (Skinner et al., 2009) from URGI (Unité de Recherche Génomique Info/research unit in genomics and bioinformatics). Annotations of found overlapping candidate genes were obtained from IWGSC Annotation v1.1. For candidate genes with unavailable annotations, the orthologous genes from Triticum urartu and Aegilops tauschii with highly similar sequences were considered for prediction of gene function 1 https://github.com/YinLiLin/R-CMplot 2 http://plants.ensembl.org/Triticum_aestivum/Info/Index in wheat by implementing BLASTn 3 . Further, annotations for SNPs (e.g., intergenic variants) with no overlapping gene were searched using a similar approach. The annotations were then confirmed against the protein sequences for determining putative molecular functions in wheat using BLASTx with default parameters in the Blast2Go v5.2 tool (Conesa and Götz, 2008). This could provide aid in the identification of putative candidate genes for disease resistance.

Phenotypic Evaluation of Seedling and Adult Plant Stage
The phenotypic data recorded for disease response to pathotypes of three Puccinia species on wheat at seedling stage (IT) under controlled conditions and adult plant stage under field conditions (IR, DS, and COI) have been provided in Supplementary  Table S3. In the panel, rust gene postulation was successfully done in about 240 genotypes based on response to multipathotype testing at seedling stage. Known resistance genes for the three rusts were postulated mostly in varieties (139) followed by improved genotypes (93) and few genetic stocks (12) (Supplementary Table S3).

Seedling Response to P. stritiformis Pathotypes
On the linear scale of 0-9, IT scores ranged from 0 (as most resistant) to 8 (as susceptible) for the Pst pathotypes. Of the 5 Pst pathotypes, based on IT scores, 93, 63.7, 61.4, and 56.9% of the tested lines were found susceptible toYR_110S119, YR_238S119, YR_110S84, and YR_T, respectively (Table 1 and Figure 1A). This indicates YR_110S119 being the most virulent in the current study panel. On the other hand, pathotype YR_46S119 was found to be avirulent with 48.1% of the tested lines as resistant and 36.6% as susceptible. The mean IT scores for YR_110S119, YR_238S119, YR_110S84, YR_T, and YR_46S119 were 7.56, 5.64, 5.26, 5.11, and 3.51, respectively. Except for YR_46S119, the IT scores for other pathotypes were skewed toward susceptibility (Supplementary Table S4). The phenotypic correlation coefficients among pairs of the five pathotypes were observed to be significant (at P < 0.01, 0.001) and in the range of 0.09 (YR_110S119-YR_T) to 0.23 (YR_46S119-YR_110S84) (Supplementary Table S5 and Supplementary Figure S1).

Adult Plant Field Response to Stripe Rust
Infection response (IR) with an average score of 0.5 was recorded at environments YR_E1 and YR_E2, whereas at YR_E3 and YR_E4, it was 0.6 (Supplementary Table S3). The distribution of IR was found to be non-biased for any single response category (Figure 2A). An average disease severity (DS) score of 62.5 at YR_E3 was found maximum among the four environments. DS with scores of 60-80 was found in most of the tested lines in all environments ( Figure 2B and Supplementary  Table S3). Based on the coefficient of infection (COI) as the disease response at four environments (COI_YR_E1-E4), stripe rust-resistant lines were found to be falling in more than one environment. In the RAMP, 37.9, 45.0, and 45.1% of the genotypes were recorded as resistant with a COI score from 0 to 20, in environments COI_YR_E1, COI_YR_E2, and COI_YR_E4, respectively (Table 1 and Figure 3A). Highly significant (at P < 0.001) Pearson's correlation coefficients between IR, DS, and COI across four environments showed positive correlations (Supplementary Table S6 and Figure 3B). Strong correlation coefficients (R = 0.8-0.9) were observed between COI with IR and DS within the environment, whereas, it ranged from 0.6-0.8 between IR and DS. For COI among environments, strong values of R were detected between YR_E1 with YR_E2 and YR_E3 with YR_E4.

Adult Plant Field Response to Leaf Rust
An average high score of IR (0.8) and DS (21.7) was recorded in environment LR_E2 (Supplementary Table S3). The IR with score 1.0 (S) was seen in most of the genotypes irrespective of the environment (Figure 2C). DS with pustules coverage score 0-20 was present in more than 300 genotypes in each environment ( Figure 2D). Likewise, more than seventy percent of the tested lines were observed to be resistant in the panel as per the COI score in all four environments (COI_LR_E1-E4) ( Table 1 and Figure 3C). All correlations across four environments for IR, DS, and COI were found positive and highly significant (at P < 0.001) (Supplementary Table S6 and Figure 3D). Interestingly, unlike most genotypes featuring susceptible IR score, a very strong correlation (R = 0.99) was found between DS and COI per environment. This may suggest that a large proportion of the tested lines were able to resist the spreading of rust pustules. IR and DS showed a correlation coefficient value of R = 0.5-0.6 and COI with a correlation coefficient of 0.7, within and among the four environments.

Adult Plant Field Response to Stem Rust
Like leaf rust, stem rust also had a similar IR frequency distribution in the panel for the score 1.0 (S) with an overall average of 0.7 in each environment (Supplementary Table S3 and Figure 2E). Pustules coverage was found in a limited number of genotypes beyond the DS score of 40, where only 4 genotypes (HTW-6, IC212182, IC28617, and LGM69) attained rust coverage of 100 (Supplementary Table S3 and Figure 2F). The average DS was higher in SR_E3 followed by SR_E1. In the study panel, ∼62% of the tested lines were observed resistant to stem rust in environment COI_SR_E1 and E4 (Table 1 and Figure 3E). Results of correlations between DS and COI were found similar to those for leaf rust field response, whereas IR and DS showed R ranging from 0.4 to 0.6 within the environment. It varied from 0.5 to 0.7 for COI across the environments (Supplementary Table S6 and Figure 3F). All correlation coefficients were observed to be highly significant at P < 0.001. For each of the three rusts in four environments, it was observed that the susceptible check lines present at regular intervals of the test genotypes showed maximum disease response with a score of 100S. This indicates that any possible variations observed in the reaction of the genotypes in RAMP would be genetic in nature and not due to the environment or inoculum load. The disease response at seedling and adult plant stages were observed to have less than 10% missing observations in the study panel. The distribution of phenotypic disease response (IT and COI) were found deviating from the normal distribution. No improvement in the normality was observed after using logarithmic and square root transformation. Therefore, the original data were used for subsequent genetic analyses.

Variance Components Estimation and Broad-Sense Heritability (H 2 )
Estimation and analysis of variance components, using a linear mixed model approach, revealed highly significant (P < 0.001) differences for IR, DS, and COI among the genotypes (g) across all environments in each rust type. Significant differences for location (e) and genotype-location interactions (g × e) were not observed in all cases (Table 2), therefore, by and large, indicating ample inoculum load and congenial atmospheric conditions for the appearance of the disease. High H 2 estimates were observed across the four environments for each rust type, ranging from 0.78 to 0.92 for IR, 0.78 to 0.90, and 0.83 to 0.90 for COI ( Table 2). The first two principal components explained together about 95.05, 88.68, and 86.76% of the total variation in GGE Biplot analysis for YR, LR, and SR field response, respectively. In reference to ranking resistant genotypes across environments, 22.15% (for YR), 61.49% (for LR), and 43.89% (for SR) were found in proximity to the ideal genotype (Supplementary Figure  S2). In this study, a genotype would be considered as an ideal genotype, which would show a uniformly low and stable disease response in the tested four environments for all the three rusts.

SNP Markers, Population Structure, and Linkage Disequilibrium
Of the 35K SNP markers, 14,650 polymorphic mapped markers were used for LD analysis and subsequently for GWAS. Genome D was found evident with lower marker density as compared to A and B genomes. Population structure analysis on RAMP was previously described by Kumar et al. (2020). The two subpopulations observed comprised of 106 (in SP1) and 377 (in SP2) genotypes (Supplementary Table S3). Genotypes that have been involved mostly in the breeding programs were observed as a major subpopulation. Principal component analysis (PCA) concurred with the result of STRUCTURE. Two separate clusters were observed in the PCA where PC1 (principal component 1) and PC2 contributed with 13.04 and 4.11%, of the total variation, respectively (Supplementary Figure S3). The effect of population structure was evident in the seedling response to multiple rust races. All pathotypes had a significant mean difference (at P < 0.05) for subpopulation based seedling disease response except for the pathotypes YR_110S119, YR_238S119, and YR_46S119 (Supplementary Table S7). The mean of IT scores for genotypes in subpopulation one (SP1) was higher than that of those in subpopulation two (SP2). This shows that in the current study genotypes were more resistant in SP2 than in SP1 to most of the pathotypes. A similar trend was observed in adult plant disease response where SP2 was significantly more resistant than SP1 based on COI (Supplementary Table S7). Therefore, the influence of population structure on rust infection was considered as covariates for subsequent association analyses.
genome. The significant marker pairs at P < 0.01 were considered for the study. Genome A had 69.56% of significant marker pairs with an average r 2 − value of 0.250. Chromosome 2A had the highest (21.2%) and chromosome 4A had the lowest (10.1%) percentage of significant marker pairs of LD r 2 . Genome B had 73.67% of significant markers with an average r 2 − value of 0.212. All chromosomes of B genome had a similar percentage (∼16-18%) of significant marker pairs except for chromosome 4B (5.2%) and 7B (9.6%). Chromosome 1D (31.9%) and 2D (29.4%) had the highest number of significant marker pairs from the D genome with chromosome 4D (2.7%) bearing the lowest count ( Figure 4A).
The critical r 2 -values were 0.126, 0.219, and 0.215 for genome A, B, and D, respectively. The difference in LD decay analysis was observed between A, B, and D genome. Whereas, predicted genome-wide LD decay was below critical r 2 = 0.201. As LD beyond the critical value is considered to be caused by genetic linkage, chromosome-wise significant marker pairs with r 2 > 0.201 were also estimated (Figures 4B-D). At a genomewide scale, 36.35% of significant SNP pairs showed LD beyond critical r 2 (> 0.201). Among the marker pairs in LD due to linkage, chromosome 5B contained the highest percentage (11.96%) of these markers, while chromosome 4D contained the lowest percentage (<< 1%). The percentage of these marker pairs in the A, B, and D genomes were 39.35, 47.21, and 13.43%, respectively. The map distance at which the fitted decay curve intersected with the critical r 2 provided an estimate of QTL-CI. For genome A, B, and D, the estimated QTL-CI was 5, 3, and 8 cM, respectively (Figures 4B-D).

Marker Trait Associations for the Three Rusts Field Response
Under natural disease pressure in respective environments, 78, 151, and 252 significant (at −log 10 p ≥ 3) MTAs were observed using the CMLM for YR, LR, and SR field response in adult plant stage, respectively (Figure 6, Supplementary Table S10, and Supplementary Figure S5). In silico annotation of these MTAs associated with adult plant disease response are also reported. Several resistance associated genes were found in adult plant responses such as E3 ubiquitin-protein ligase, NB-LRR protein, Lr10 disease resistance locus receptorlike protein kinase, serine/threonine kinase protein, DEADbox ATP-dependent RNA helicase, disease resistance protein RGA3, etc. (Supplementary Table S10). LD based putative QTL identification revealed 29 distinct loci on all chromosomes except for 3D, 4D, 5D, and 7B to be associated with YR response, where 16 loci were observed in more than one environment. In case of LR response, 45 distinct loci were observed on all chromosomes except for 2D, 4D, 5D, and 7A. For SR response, none of the 44 distinct loci were observed on chromosomes 2D, 5D, and 7D. Eighteen and twenty-seven loci were observed in more than one environment for LR and SR, respectively ( Table 4 and Supplementary Table S11). Among the identified representative SNPs, haplotype analysis using Haploview v4.2 (Barrett et al., 2005) markers with LD value r 2 > 0.8 were observed in the case of LR and SR but not in YR. For LR, only one set of marker pairs comprising two markers was in strong LD, whereas for SR, three sets of marker pairs comprising of 2, 7, and 11 SNPs were observed in LD (Supplementary Figure S6 and Supplementary Table S11).
For YR field response, QYr.ramp-2D.1 was the only locus observed in all the four environments (YR_E1-E4), whereas the remaining 15 loci represented at least two environments. Most   and (C) SR environments SR_E1, SR_E2, SR_E3, and SR_E4 were plotted from inside to outside, respectively. A multi-track Q-Q plot for each case is presented at the right upper corner of the circular Manhattan plot. The threshold value at -log 10 (p) ≥ 3 is indicated as red-colored circle for each environment. A rectangular version of the plots is shown in Supplementary  Figure S5. significant association was seen for QYr.ramp-6B.4 in YR_E2 followed by QYr.ramp-1B.4 in YR_E4. For LR field response, three loci were seen associated in all four environments (LR_E1-E4) namely QLr.ramp-4A.4, QLr.ramp-6B.3, and QLr.ramp-7B.1. Most significant association was found for QLr.ramp-1B.3 in LR_E3 which also holds the most numbers of MTAs observed. In the case of SR, the most number of loci (nine) affiliated with all four environments (SR_E1-E4) were observed. Amongst these and other loci, QSr.ramp-3D.4 holds the most number of MTAs followed by QSr.ramp-1B.4, where the most significant association was observed for QSr.ramp-3B.7 in SR_E3 ( Table 4).
Among YR field response associated loci, QYr.ramp-6B.4 and QYr.ramp-1B.4 exhibited similar and higher contribution for COI with R 2 ranging from 2.42 to 3.55% in at least two environments. Thirteen loci were associated with YR field response in single environment, including 6 loci in YR_E4, 3 loci in YR_E3, and 2 loci each in YR_E2 and YR_E1 (Supplementary Table S11). For LR field response, of 45 loci, QLr.ramp-1B.3 accounted for the higher level of R 2 of 5.03% while contributing to COI. To LR field response in single environment, 27 loci were seen associated where 10 loci in LR_E1, 7 loci in LR_E2, 6 loci inLR_E3, and 4 loci in LR_E4 were observed (Supplementary Table S11). In case of SR field response, QSr.ramp-3B.7 explained the higher level of contribution for COI with R 2 ranging from 1.96 to 4.18%. Seventeen loci were associated with a single environment, where 8 loci in SR_E1, 4 loci each in SR_E2 and SR_E4, and single loci in SR_E3 were observed (Supplementary Table S11).
Like seedling response, some representative SNPs were found common between two rust type field responses. Among those SNPs, AX-94638655 and AX-95176310 representing loci for YR and LR; AX-94874313 and AX-94877284 representing loci for LR and SR showed a bidirectional allelic effect. Whereas, AX-94805944 representing loci for LR and SR showed unidirectional allelic effect (Supplementary Tables S10, S11). Despite having different representative SNP markers, certain QTL observed for resistance in the adult plant stage were seen to be sharing QTL-CI based similar to identical genomic region with those observed for seedling stage resistance (Supplementary Table S12). In such case, QTL observed for both seedling and adult plant disease responses were considered as ASR loci and those observed only for the adult plant disease response were considered as APR loci. There were few ASR loci which had identical representative SNP marker for the QTL observed in seedling and adult plant stages. Of 12 YR ASR loci, QYr.ramp-3B.4 was represented by SNP AX-94680284. Similarly, of 15 LR ASR loci, QLr.ramp-3B.3 and QLr.ramp-3D.4 were represented by SNPs AX-94741529 and AX-94874313, respectively. In the case of SR, with the most number of ASR loci (31), there were fourteen such representative SNP markers (Supplementary Table S12). The most number of APR loci were observed against LR (30) disease response followed by that of YR (17) and SR (13). Several genomic regions in terms of QTL-CI were found common in more than one rust field response. These robust genomic regions bearing QTL-CI based co-localized ASR and APR loci can be harnessed for further exploration and usage ( Table 5). Among these genomic regions, seven of them been found on

Pyramiding Effects of Favorable Alleles on Field Response
The cumulative effect of favorable MTAs for field resistance was studied based on the estimated number of favorable alleles in each specific genotype. For this, SNPs representing the distinct loci for each rust were considered. Alleles associated with a reduction in disease response were considered as favorable alleles at each locus of the representative SNP. With 29 SNPs for YR, the number of favorable alleles ranged from 6 to 26. Similarly, for LR and SR, it ranged from 10 to 40 and 4 to 39 with 45 and 44 SNPs, respectively (Supplementary Table  S13). The accumulating resistant alleles were expected to be strongly correlated with the disease response. A significant negative correlation (at p < 0.0001) was observed in the case of YR (-0.4825), LR (-0.6256), and SR (-0.5775) between the number of favorable alleles in individual genotypes and averaged disease response under their respective environments. This observation was further supported by the linear regression coefficients (R 2 ) found significant at p < 0.0001 in all three cases (Figure 7). The total SNPs (118) estimated as field response based QTL representative markers were examined using a hierarchical clustering approach by selecting 72 genotypes having average COI score ≤ 20 in all the three rusts. Clearly, these genotypes were found to have the most number of favorable alleles compared to moderately resistant to susceptible genotypes except for some resistant genotypes viz. E3257, H957, HUW679, PBW138, IC321918, and LBRL4. The genotypes were grouped in same clusters based on the SNP array pattern having favorable alleles in those genotypes (Supplementary Figure S7). Pedigree based closeness was also visible in genotypes showing allelic similarity in terms of favorable alleles, e.g., HPW360-HPW361, DBW50-DBW88, and VL802-VL804. The present finding provides valuable information to the wheat breeders in deciding the selection of parents for the crossing program, particularly those to be utilized as source of rust resistance. Incorporation of diverse alleles can be ensured if parents are chosen from the different clusters (Supplementary Figure S7).

DISCUSSION
Crop scientists face dual challenges of increasing productivity and sustaining it by managing stresses occurring due to biotic and abiotic factors. Among biotic factors, the three rusts of wheat cause continuous threat through the rapid evolution of new races. The development and deployment of resistant cultivars is the most economical and environmentally safe method to control diseases. Primary gene pool including indigenous collections comprising of landraces, old cultivars and breeding lines are considered as a valuable genetic resource for providing new and durable resistance that can be exploited for the development of current day high-yielding varieties (Mujeeb-Kazi et al., 2013). Utilization of primary gene pool is advantageous over secondary (comprising of wild relatives) as they carry homologous genomic regions having a better opportunity as combiner with hexaploid wheat (Wulff and Moscou, 2014).
To widen the spectrum of rust resistance in modern wheat varieties, finding lines having new sources of resistance along with newer alleles for known resistance genes, is the way forward. High throughput genotyping technologies coupled with bioinformatics enable us to mine valuable information from genome database paving the way for a better understanding of vast germplasm collections. With this objective, a diverse panel of 483 genotypes found suitable for GWAS (Kumar et al., 2020) was utilized for identifying genomic regions associated with rust resistance in the form of diverse resistant alleles and their combinations for durable resistance. This panel comprised diverse genotypes both in temporal and spatial dimensions, i.e., genotypes adapted to various agro-climatic conditions spanning India and varieties released during pre and post green revolution era. This panel encompassing of hitherto less engaged germplasm such as landraces offers a rich and new accessible source for disease resistance. Studies have supported the utilization of landraces for exploring novel resistant alleles (Muleta et al., 2017b;Pasam et al., 2017;Riaz et al., 2018). The same diverse panel was evaluated against all the three rusts without aiming any biasness toward specific rust. The phenotypic distribution of disease response did not follow a normal distribution, and original phenotypic data was utilized as such, similar to the previous studies Edae et al., 2018). In this study, more than 50% of genotypes in the panel showed resistance to five out of seven SR pathotypes. Similar proportion of genotypes were found resistant to SR in four environments at the adult plant stage under field disease situations. This observation strongly suggested the presence of ASR loci in the panel. Whereas, in the case of LR, the presence of APR loci was supported by the fact that more than 70% of the lines were resistant in each of the four LR environments in spite of the panel showing susceptibility to five out of six LR pathotypes under seedling stage disease response. This was in concurrence with most number of ASR and APR loci observed against SR and LR disease response, respectively. Significant variations were observed among the genotypes in the panel as indicated by analysis of variance (ANOVA). High heritability estimates observed in this study for field disease responses were supported by significant positive correlations among traits and occurrence of consistent data recorded across environments. This indicated reproducibility of the data reliable for GWAS (Godoy et al., 2018).
Wheat breeder's 35K array derived from wheat 820K SNP array provides breeder oriented informative markers (Allen et al., 2017). It has an added advantage of capturing information on germplasm derived from secondary and tertiary gene pools as the array is developed using exome capture of wheat and its wild relatives (Allen et al., 2017;Rasheed and Xia, 2019). In our study, LD decay was found strongly different in the three genomes with B genome to have the quickest decay followed by A genome, which is supported by previous reports Riaz et al., 2018). This supports the practical use of D genome markers suitable for association tests even though less in number . D genome was found to have the highest LD than other genomes which has also been mentioned in several previous studies (Nielsen et al., 2014;Wang et al., 2014;Zegeye et al., 2014;Voss-Fels et al., 2015;Riaz et al., 2018). The number of significant marker pairs was found lowest in the chromosome 4D similar to a previous study (Muleta et al., 2017a).
Population structure in the panel showed two subpopulations (Kumar et al., 2020) which was further supported by PCA. A significant effect of population stratification was observed in differences based on disease response at seedling and adult plant stage. The resistant subpopulation of the two subpopulations, comprised of majorly improved genotypes and released varieties. It appeared that landraces in the minor subpopulation were found mostly susceptible. The landraces under study have originated from central and northwestern states of India namely Gujarat, Rajasthan, and Madhya Pradesh (Kumar et al., 2020). Landraces adapted to these regions are mostly tolerant to abiotic stresses viz. drought, heat, and salinity, e.g., genotype KHARCHIA LOCAL (origin: Rajasthan, India) is adapted to long-term salinity stress (Mahajan et al., 2017(Mahajan et al., , 2020. We hypothesize that landraces adapted to harsh environments, such as abiotic stresses and low-input conditions give an eluded appearance of diseases. But under high input conditions, we observed that landraces appeared to be more susceptible when compared to improved genotypes. Since not all landraces were found susceptible (Supplementary Table S3), they may still be beneficial in the identification of novel resistant alleles. Of 96 landraces in the panel, twenty-one were found resistant in field response against LR and eleven each against YR and SR. Landrace IC321998 (origin: Uttarakhand, India) was found resistant for all three rusts. This suggests an indirect selection for disease resistance in spite of the fact that landraces are traditionally selected by the farmers preferring better agronomic traits (Zeven, 2002).
Care was taken to avoid spurious association due to population stratification by treating it as a covariate along with kinship accounting relationship among genotypes for association studies in this panel (Stich et al., 2008;Lipka et al., 2012;Yang et al., 2014). The selection of CMLM applied in this study was found suitable from the comparative study of different association models. In most of the cases, based on the Q-Q plots, CMLM was found to have minimum deviations in observed probabilities with respect to expected probabilities (Supplementary File S1). Several studies have been successfully conducted using the CMLM model (Arruda et al., 2016;Muleta et al., 2017a,b;Sheoran et al., 2019). It also provides efficiently high statistical power over and is an improvement to the MLM model for a large sample size (Zhang et al., 2010).
Three ASR loci QYr.ramp-1A.2, QLr.ramp-1A.2, and QSr.ramp-1A.2 were found co-localized on chromosome 1A against each rust located at 72.25-75.74 cM. Except for QYr.ramp-1A.2 exhibiting minor effect, other loci were found associated with disease response to multiple pathotypes and environments. This was in agreement with the QTL observed to have a minor effect but also showed seedling resistance on chromosome 1AL (Rosewarne et al., 2012;Jighly et al., 2015). QLr.umn-1AL was reported to be associated with both field and seedling disease response mapped at 149.8 cM on chromosome 1AL . In this study, QLr.ramp-1A.2 was observed proximal to this previously reported QTL. APR for LR is not known to be present on chromosome 1AL  and the chance of having an association with Lr59 is highly unlikely as source of such alien introgression (Marais et al., 2008) is not present in the current study panel. A MTA exhibiting seedling resistance against Pgt race TTKSK represented by DArT marker wPt-2014 was reported on chromosome 1AL mapped at 70.5 cM (Laidò et al., 2015). This could be present in close proximity to the loci QSr.ramp-1A.2 identified in this study.
LR APR loci QLr.ramp-1B.2 was mapped on the short arm of chromosome 1B located at 8.24-9.93 cM. It was identified in close proximity to SSR marker Xswm271 linked to known APR gene Lr75 (Singla et al., 2017) and proximal to another APR gene Lr71 . Genomic region 23.90-26.22 cM (1B) was found to have a large number of MTAs in all the three rusts for seedling disease response ( Figure 5 and Table 3). This region located on chromosome 1BL harbors several known Yr genes such as Yr9, Yr10, Yr15, Yr 64, Yr65, etc. In concurrence ASR loci QYr.ramp-1B.3 was found in a close proximity to Xbarc119 (27.4 cM) and Xgwm413 (29.9 cM) which are linked to Yr15 and YrH52 (Maccaferri et al., 2015a). This region also showed a possible association with wheat-alien translocation (1RS.1BL) derived from Secale cereale (Schlegel and Korzun, 1997). Resistance genes for all the three rusts are known to be associated with this translocation viz. Yr9 (Ren et al., 2009), Lr26, andSr31 (Mago et al., 2005). Race-specific gene postulations demonstrated the presence of these genes in the panel, possibly being reflected as the co-localized ASR QTL QYr.ramp-1B.3, QLr.ramp-1B.3, and QSr.ramp-1B.4. The most virulent Pgt race Ug99 and its lineage have not been observed in India and its neighboring countries inferring the effectiveness of Sr31 in Indian germplasm (Prasad et al., 2018;Bhardwaj et al., 2019b).
So far, five QTL have been reported for resistance against YR on chromosome 1D, with three present at the distal end (Zwart et al., 2010;Ren et al., 2012a;Vazquez et al., 2012) and two present at the proximal end of its short arm (Maccaferri et al., 2015b;Godoy et al., 2018). Not many reports are present for ASR genes against YR at locations close to the centromere on chromosome 1DS where Godoy et al. (2018) reported QTL present at 89.58 cM. Two potential ASR loci QYr.ramp-1D.1 (AX-94508418: 1.03 Mb) and QYr.ramp-1D.3 (AX-94444583: 34.6 Mb) were observed both at distal and proximal regions of 1DS, respectively, in this study. The chromosome 1DS is also known to harbor Lr42 gene conferring resistance against LR at both seedling and adult plant stages and recently it has been fine mapped on 1DS flanked by markers TC387992 and Xwmc432 to a 3.7 cM genetic interval (Gill et al., 2019). QLr.ramp-1D (3.38 cM) was observed very close to this location at adult plant stage only. QSr.ramp-1D.1 located at 1.69 cM on 1DS was observed in close proximity of gene Sr33 which is known to be flanked by Xbarc152 and Xcfd15 at a distance of 1.8 cM on either end on 1DS (Sambasivam et al., 2008). The gene was introgressed in bread wheat from its wild relative Aegilopes taushii providing resistance against Ug99 group races (Periyannan et al., 2013).
The genomic region 102.12-104.59 cM was observed to have ASR loci for each rust. This region was present on the chromosome 2BL supported by the physical positions of QTL associated SNPs spanning ∼612-722 Mb region. Corresponding to this physical region, Aoun et al. (2019) identified QSr.ace-2B conferring resistance to Pgt races TTKSK and JRCQC (virulent to Sr9e) linked to Sr9 gene flanked by SNP markers IWA6399 and IWA7955 (92.7-105.3 cM) on chromosome 2BL. SR ASR loci QSr.ramp-2B.6 was identified within this genomic region supporting possible affiliation to Sr9. In co-localization to which YR ASR loci QYr.ramp-2B.4 was identified in a single field environment possibly due to its minor effect. It was located 5 cM proximal to another minor effect QTL QYrAvS.wgp-2BL . It was further located 3 cM proximal to SNP marker IWA638 (724 Mb: 107.58 cM) (Blake et al., 2016) which has been identified in a close proximity of YrSp gene (Feng et al., 2015). Another ASR loci QLr.ramp-2B.7 against LR was identified at 103.81-104.59 cM on chromosome 2BL. Gao et al. (2016) reported a QTL 2B_3 at 102.28-108.35 cM interval on chromosome 2B in response to field experiments only which reduced any possible similarity. Two formally known LR resistance genes Lr50 and Lr58 are present on the terminal region of chromosome 2BL (Brown-Guedira et al., 2003;Kuraparthy et al., 2007). SSR markers Xgdm87 (Lr50) and Xcfd50 (Lr58) linked to these genes are present distant to the identified QTL QLr.ramp-2B.7 based on the integrated consensus map (Maccaferri et al., 2015a). In a previous report, multiple resistance QTL against YR and LR QYr.ifa-2BL/QLr.ifa-2BL was identified on chromosome 2BL. It was flanked by DArT markers wPt-3378 and wPt-7360 (Crossa et al., 2007;Buerstmayr et al., 2014). However these markers were also present distant (Maccaferri et al., 2015a) to QLr.ramp-2B.7, leaving it for being an unexplored ASR loci.
Two APR and one ASR loci were identified co-localized on chromosome 6AL located at 215.37-221.75 cM (∼605-615 Mb). QYr.ramp-6A.5 was the ASR loci which provided resistance against YR in this region. It has been found distal to a previously reported ASR QTL tagging SNP IWA2129 on chromosome 6AL mapped at 212.2 cM (596.89 Mb) (Blake et al., 2016;Muleta et al., 2017a). Chromosome 6AL is also known to harbor the YR ASR gene Yr38 which is linked with another ASR gene Lr56 for LR (Park, 2016). None of these genes were postulated in the current study panel. In a previous study, Lr64 gene was mapped on chromosome 6AL tightly linked to the KASP marker which was developed for the SNP marker IWB59855 . The LR APR QTL QLr.ramp-6A.4 identified was represented by SNP marker AX-95147877 (615.46 Mb) in this study, found to be present close to SNP marker IWB59855 (614.17 Mb). Therefore, QLr.ramp-6A.4 could be strongly associated with the gene Lr64. Interestingly, Lr64 has not been characterized as ASR or APR gene so far (Park, 2016), whereas QLr.ramp-6A.4 was identified exhibiting APR response. This could provide supporting evidence for the identification of the true resistance nature of Lr64 in the future. An ASR QTL QSr.fcu-6A was mapped earlier on chromosome 6AL to the genomic region harboring Sr13 gene conferring resistance against Pgt races TMLKC, TRTTK, and TTKSK in emmer wheat . An SNP-based semi-thermal asymmetric reverse PCR (STARP) marker rwgsnp7 was developed from the SNP IWB34398 (6AL: 615.45 Mb) to map QSr.fcu-6A. Validation of marker rwgsnp7 in hexaploid wheat provided evidence that chromosome 6D harbors a homoeologous allele similar to that present on chromosome 6A for this marker . In this study, QSr.ramp-6A.5 was identified as APR loci mapped at 215.37-221.75 cM (∼605-615 Mb) interval on chromosome 6AL where the SNP IWB34398 overlaps with this region. Therefore, QSr.ramp-6A.5 could be present in the genomic region harboring Sr13. Since none of the Pgt races mentioned above were used in this study, seedling resistance could not be confirmed governed by this region in this study.
Chromosome 6B was observed to have co-localized APR loci for the three rusts. Based on MTAs associated with these QTL, an estimated physical location of ∼712-714 Mb was obtained (IWGSC RefSeq v.1.0) representing the long arm of the chromosome. Among several other reported YR APR QTL on 6BL (Mu et al., 2019), three QTL viz. QYr.wsu-6B.1 (Bulli et al., 2016), QYr.ucw-6B (Maccaferri et al., 2015b), and QYr.wgp-6B.1 (Santra et al., 2008) spanning ∼631-720 Mb genomic region were found overlapping to some extent with QYr.ramp-6B.4 identified in this study. LR APR QTL QLr.ramp-6B.6 was identified distally close to a recently reported APR QTL QLr.cim-6BL mapped in a durum mapping population on chromosome 6BL . It was located in the interval of SNP markers AX-95155193 and AX-94562707 spanning 691.07-698.26 Mb genomic region. Lr3 (Kolmer, 2015) and Lr9 (Schachermayr et al., 1994) are well known genes present on chromosome 6BL. QLr.ramp-6B.6 was unlikely to be linked with either of these genes as Lr3 is an ASR gene and Lr9 has not been deployed in the current study panel. Similarly, Sr11 gene is also present on chromosome 6BL but exhibits ASR (origin: Triticum turgidum ssp. durum) (Park, 2016). It was mapped within the marker interval of Tdurum_contig55744_822 (IWB72471: 709.53 Mb) and BS00074288_51 (IWB10724: 715.97 Mb) (Blake et al., 2016;Nirmala et al., 2016). This region was found to be overlapping with those of QTL QSr.ramp-6B.3. Several genotypes were postulated to have this gene in the study panel. Since Sr11 provides a race-specific ASR against Pgt race TKTTF which has not been used in this study, QSr.ramp-6B.3 could have been mistakenly identified as an APR loci. Although, without further investigation a conclusion cannot be made.
The region (17.65-18.61 cM) harboring one ASR loci against YR and SR each, and one APR loci against LR were identified on the long arm of chromosome 6D (∼459-460 Mb). Gebrewahid et al. (2020) reported an APR QTL QYr.hebau-6DL located at 464.6-472.0 Mb interval flanked by SNP markers AX-108848475 and AX-109273869. QTL QYr.ramp-6D.2 identified in this study was present proximal to QYr.hebau-6DL, showing ASR property. There was a clear difference in the nature and position of these two QTL. LR APR QTL QLr.ramp-6D.1 was observed in a single environment depicting a possible minor effect. Although, chromosome 6DL harbors a known ASR gene Lr38 inheritably linked with SSR marker Xwmc773 (Mebrate et al., 2008) could be present relatively close to QLr.ramp-6D.1. Of the four known Sr genes on chromosome 6D, only the ASR gene Sr29 is present on the long arm (Dyck and Kerber, 1977). In a report, an APR QTL QSr.cim-6DL was mapped on chromosome 6DL. In our study, QSr.ramp-6D was identified as an ASR loci which could be in close proximity of Sr29.
The telomeric region of chromosomes 2AL and 7DS bears significant importance in YR resistance along with other resistant systems for barley yellow dwarf virus, powdery mildew and fusarium head blight (Boukhatem et al., 2002). Apart from the three known YR seedling resistance genes viz. Yr1 (Bansal et al., 2009), Yr32 (Eriksen et al., 2004), and YrJ22 , a recently identified and temporarily designated YR seedling resistance gene YrH9017 (Dong-fang et al., 2019) has been found on chromosome 2AL. None of the genotypes in the panel harbored these genes based on their pedigree and racespecific gene postulations. Instead, an unexplored YR APR loci QYr.ramp-2A.4 (179.61 cM) was observed at the telomeric region represented by SNP AX-94463225 (764.1 Mb). In co-localization to which, a LR APR loci QLr.ramp-2A.5 (AX-94995671: 779.6 Mb) was observed 33.07 cM and 73.07 cM distal to APR QTL QLr.ubo-2A  and QLr.hebau-2AL (Zhang et al., 2017), respectively, in this study. APR loci QYr.ramp-7D.2 and QLr.ramp-7D.2 on chromosome 7D (6.96-18.37 cM) were found associated possibly with pleiotropic resistance gene Lr34/Yr18/Sr57, where SNPs representing these QTL viz. AX-95205886 (5.38 Mb) and AX-94674448 (72.94 Mb), respectively being in high LD (D' = 0.87) are present on the short arm of chromosome 7D based on the IWGSC RefSeq v1.0 (Appels et al., 2018) reference physical map. The association to this pleiotropic gene can be further supported by the fact that forty genotypes in the current study panel were characterized for Lr34 in a previous report (Pawar et al., 2013). When compared to 90K SNP array, 35K SNP array used in this study was able to capture association with this pleiotropic resistance gene, since none of the 90K SNP markers are known to be linked to this gene (Muleta et al., 2017b). Since this array is capable of capturing conserved regions associated with yield-related traits present on the 1RS.1BL translocation region (Luján Basile et al., 2019) and was designed using wild relatives (Allen et al., 2017), provided an advantage over the most utilized 90K SNP array.
Evidence of true QTL can be derived from the results of having significant associations seen in at least two pathotypes or environments for a rust type (Edae et al., 2018). Validation of markers associated with these QTL can be done by comparing them with the positions of the markers of previously mapped QTL (Pasam et al., 2012;Hwang et al., 2014;Kertho et al., 2015;Maccaferri et al., 2015a). In this study, MTAs found associated with multiple pathotypes may be broadly useful in the enhancement of future germplasm and breeding programs (Liu et al., 2017). Accumulation and detection of favorable resistance alleles often pose challenges in the breeding program upon consideration (Riaz et al., 2018). In our study, a large number of favorable resistant alleles were detected in some landraces and these showed moderate to resistant behavior against the three rusts. As a separate observation, SNP AX-94394356, AX-94674448, and AX-95176310 being the QTL representing markers in field response against LR, had a poor percentage of resistance alleles accumulated in improved genotypes and released varieties when compared to landraces, in the panel. Similarly, AX-94608698, AX-95226309, AX-94485764, and AX-94805944 observed for SR field response, had more accumulated favorable alleles in landraces than that in other genotype categories. Such observations were not distinctively visible in case of YR response. This observation is possible when these favorable alleles have accumulated in landraces over time with varying frequencies (Riaz et al., 2018). The concept of 'Speed breeding' introduced (Watson et al., 2018) with the aim of utilizing latest advancements in phenotyping and high throughput data generation, can be helpful in developing rust-resistant wheat cultivars at an accelerated rate (Hickey et al., 2012;Riaz et al., 2016). To determine the identified QTL to be potentially novel or are associated with previously mapped QTL, allelism tests will be required.
Breeders engaged in wheat genetic improvement have a keen interest in looking for desirable parents or donor lines to enrich their crossing block. While grain yield is prime focus area along with end product making quality parameters, yet stabilization of yield potential through tolerance against abiotic and biotic stresses will always be a major concern. Therefore, any information on diverse or novel sources of resistance, in this case against the three rusts, is of considerable importance to the breeders. Since wheat cultivation in India corresponds to many of the agro-climatic regions growing wheat globally, the information becomes universally important in combating rust diseases. Our finding of hitherto less utilized sources of stripe, leaf and stem rust viz. IC321918, IC212153AMB, IC321998, H957, LBRL4, and some exotic lines shall attract the attention of the wheat breeders. When present alone, the APR genes/QTL do not confer adequate resistance especially under high disease pressure; however, combinations of 4 or 5 minor genes usually result in "near-immunity" or a high level of resistance . Increasing the frequency of QTL for rust resistance through divergent crossing in the population improvement program of wheat can provide durable rust resistance by broadening the genetic base for resistance. The pleiotropic QTL can confer slow rusting resistance against the three rusts of wheat. Bi-parental and multi-parent populations developed using resistant landraces and modern varieties shall enable toward this endeavor.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
RT, AR, and DiK conceived the theme of study. DeK, AK, SS, SJ, MI, and UA analyzed the data. DeK, VC, and RT drafted the manuscript. DeK conducted the field condition phenotyping at Karnal, India. OG and SB conducted the phenotyping of the seedling resistance test. MS conducted the field condition phenotyping at Wellington, India. SP and TP the conducted field condition phenotyping at Indore, India. HK, RS, and PS provided insights into reported genes for rust resistance. VC, RT, GS, AR, GPS, and DiK provided overall guidance and edited the manuscript. All authors read and approved the final version of the manuscript.

ACKNOWLEDGMENTS
We are thankful to the Indian Council of Agricultural Research, Ministry of Agriculture and Farmers' Welfare, Govt. of India for Advanced Super Computing Hub for Omics Knowledge in Agriculture (ASHOKA) facility at ICAR-IASRI, New Delhi, India created under National Agricultural Innovation Project, funded by World Bank at ICAR-IASRI, New Delhi. Support provided by the Department of Biotechnology, Ministry of Science and Technology, Government of India in the form of DBT-JRF fellowship program to DeK was acknowledged. Financial support from CABin Grant of ICAR-IASRI, New Delhi in the form of fellowship to AK was also thankfully acknowledged. Drs. Pramod Prasad and Subodh Kumar from the regional station, ICAR-IIWBR, Shimla for help in phenotyping against rust pathotypes. Support of Dr. Arun Gupta from the GRU, ICAR-IIWBR is highly acknowledged.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2020.00748/ full#supplementary-material FIGURE S1 | Heatmap of correlation coefficients showing the correlation between the 18 pathotypes representing the three rusts used in the study. The magnitude of the correlation is represented by the color gradient legend on the right side of the heatmap. Significant correlations are marked as * (p < 0.05), * * (p < 0.01), and * * * (p < 0.001).
FIGURE S2 | GGE Biplot representing variations explained in the first two principal components for field disease response against YR, LR, and SR. Concentric circles represent the ranking of genotypes with respect to being an ideal genotype stable in corresponding environments.       (5), stem (7), and Leaf (6) rust pathotypes.        FILE S1 | Comparative Q-Q plots of six association models for multiple rust pathotypes and four environments each for YR, LR, and SR. The CMLM was observed as the best fit model.