A meticulous process was employed to select representative genotypes from 33 diverse countries, encompassing various subgroups within the 3K rice genome dataset to develop a core sub-set from 3K-RGP. This process involved curating a set of 147 accessions (28 aus, 28 ind1A, 26 ind1B, 35 ind2, 15 ind3, 55 indx, 2 trop, 1 aro and 1 admix) that collectively represents all the subpopulation of 3K rice genome panel. These accessions were chosen to facilitate the evaluation of their resistance against blast and bacterial blight diseases. To measure the genetic diversity, an analysis was conducted utilizing 0.5 million SNPs. This involved employing the Unweighted Pair Group Method with Arithmetic Mean (UPGMA) algorithm, which computed genetic distances based on molecular markers. The outcomes of the analysis were utilized to construct a UPGMA tree, helping in the careful selection of genotypes for subsequent screening. The genetic diversity assessment confirmed the variability within the genotypes, reinforcing the suitability of the chosen core sub-set for further analysis. A core sub-set of 147 accessions of 3K RG Panel from 33 different countries which include 28 aus, 28 ind1A, 26 ind1B, 35 ind2, 15 ind3, 55 indx, 2 trop, 1 aro and 1 admix were used in this study for phenotyping blast and bacterial blight resistance.

Genetic diversity and population structure was done using genome-wide 0.5 million SNP data points of all the accessions. The tree was constructed using neighbor-joining method and the calculated dissimilarity index. Structure admixture analysis was analyzed through admixture software with k=3 as cross-validation error shown low. A Neighbor-joining tree was constructed using TASSEL Software version 5.0. The Phylogenetic tree is visualized through an unrooted tree using treeio R package. The Principal Component Analysis is done using SNPRelate R Package and Genetic Diversity and Analysis of Molecular Variance (AMOVA) is done using poppr R package.

Genome wide association study was performed using more than 0.5 million SNPs and the average SES score of three replications for blast and the average lesion length of five BLB clipped leaves. The analysis was conducted using five different multilocus GWAS models, including mrMLM (Wang et al., 2016), FASTmrMLM, FASTmrEMMA (Tamba and Zhang, 2018 ), ISIS EM-BLASSO, and pLARmEB (Zhang et al., 2017b), using an R-based methodology (https://CRAN.R-project.org/). The criterion for significance to find peak associations (MTAs) with the target trait is a LOD score of ≥5. The phenotypic allele effect (ai) was determined by the formula given by Zhang et al., 2013 and the favorable alleles of each trait.

To identify the candidate genes underlying the MTAs, all annotated genes contained within 50kbp LD from the peak MTA of the Oryza sativa reference sequence (Os-Nipponbare-Reference-IRGSP-1.0; RAP database: http://rapdb.dna.affrc.go.jp/download/irgsp1.html) were taken into consideration. Additionally, known blast R genes or literature-based QTLs were compared to the sites of the significant MTAs (Sharma et al., 2012; Ashkani et al., 2016) and recent literature (Kang et al., 2016; Ma et al., 2015, 2015; Zhao et al., 2011; Xu et al., 2014; Su et al., 2015; Mgonja et al., 2016; Zhu et al., 2016; Jiang et al., 2020; Yang et al., 2022).

The total candidate genes were taken into the haplo-pheno analysis for the identification of superior haplotypes. Haplotype analysis was carried using in-house scripts in R programming. Then the analysis included One-way ANOVA with haplotype as a fixed factor. Subsequently, Duncan’s multiple range test (DMRT) and ANOVA were used to test the statistical significance among the mean of haplotype groups using the Agricolae package in R (De Mendiburu, 2018). All reported P values underwent adjustment using either the Benjamini-Hochberg FDR or the Bonferroni method when multiple comparisons were performed, and a statistically significant threshold of adjusted p < 0.05 was considered.

For the whole investigation, we used the ‘3K filtered’ SNP set from the SNP seek database. The following filtering criteria were used to generate the filtered from the Base SNP set: (1) alternative allele frequency at least 0.01, (2) fraction of missed calls per SNP at most 0.292, and this SNP set was previously present in the SNP seek database, which was directly used in this work (http://snp-seek.irri.org/download.zul). Only nonsynonymous SNPs and indels in the exon region that cause an amino acid change were taken into account in the haplotype analysis for the candidate locus. The SNP-seek database (http://snp-seek.irri.org/_download.zul) was used to gather data on haplotypes and their diversity in order to identify the best haplotype by classifying haplotypes according to the genotype trait means for the related genes. The R-package (https://CRAN.R-project.org/) was used to create boxplots with a significance threshold of p <= 0.05.

The SNP-based diversity analysis provided valuable insights into the diversity exhibited by the panel ( Supplementary Figure 1 ). The admixture analysis showed a notable distinction at k=2, distinctly separating the indica and aus groups. The optimal population structure emerged at k=3, revealing population divergence within aus, ind1A, and ind1B ( Supplementary Figure 1A ). The structure analysis revealed that indx has a blending of all subpopulations, signifying a historical interplay of domestication and selection between indx and other clusters within the population. Remarkably, while the admixture analysis highlighted the optimal k value as 3, it couldn’t distinguish the subpopulations of ind1B and ind3. However, the cladogram analysis provided confirmed relationships and divergences within these populations ( Supplementary Figure 1B ). The cladogram, illustrating the genetic relationships between samples, suggested the findings of identity-by-descent and neighbor-joining tree construction methods. The cladogram unequivocally reinforced the admixture results, underlining how indx’s presence is a consequence of controlled selection and domestication of various varieties. Notably, indx accessions demonstrated a heightened relatedness between ind3 and ind1B, elucidating the dispersion of their unique diversity and their positioning within the broader subgroups of ind2 and ind1B.

Employing principal component analysis (PCA) on our population, we extracted 15% of the total variation, effectively representing two distinct clusters. These clusters shed light on the contributions of aus and indica lines to the population’s genetic makeup. Despite the relatively modest variance captured by the PCA analysis, the second cluster, attributed to indica lines, elucidates the divergence in allele composition in a blended manner ( Supplementary Figures 1C, D ). To further delve into genetic diversity, the poppr r package was utilized, as displayed in Supplementary Table 1 . Among the subpopulations, ind1A displayed the highest Nei’s genetic distance at 0.34, while indx exhibited the lowest at 0.23. As for the Shannon diversity index, indx emerged with the highest value of 3.9, in contrast to ind1A, which recorded the lowest at 2.6. An examination of molecular variance unveiled that 87.74% of the variance exists within samples, while the remaining 12.25% of the variance is accounted for between the distinct population groups.

Each of the accessions within the 3K-RGP sub-set underwent comprehensive screening for resistance against blast and BLB in multiple replications. Disease reaction showed minimal variation across replications, so the average score was used for the subsequent analysis ( Supplementary Figure 2 ). In terms of blast severity, the spectrum of severity scores spanned from 0, indicating high resistance (e.g., LALSAITA::IRGC 43915-1), to 5, indicative of high susceptibility (e.g., RTS 16::IRGC 8235-1). Among all accessions, 12.92% showed resistance with a severity score (SES) of 0-2, another 12.92% demonstrated moderate resistance with SES scores above 2 to 3, while the majority, comprising 72.78% of genotypes, exhibited susceptibility with SES scores above 3 to 5. The reference susceptible check HR12 recorded a score of 5, while the resistant check Tetep scored 0. Evidently, the phenotypic variation in blast severity was considerable among the genotypes within the panel, revealing a pronounced divergence ( Supplementary Figures 2A–D ).

In the case of BLB screening, all accessions were inoculated during the maximum tillering stage. Lesion length averages across the panel ranged from 1.16cm (e.g., SUWEON 311: IRGC 61890-1) to 27.2cm (e.g., BPI 76 NON SENSITIVE (GREEN)::IRGC 9790-1). The susceptible check TN1 exhibited a lesion length of 21.5 cm, while the resistant check IRBB60 displayed a length of 0.3 cm. Among the accessions, 17.07% displayed resistance with mean lesion lengths varying from 1.13 to 3.00cm. Additionally, 42% were moderately resistant, characterized by mean lesion lengths spanning 3.08 to 6.00cm. Around 14.58% were moderately susceptible, with mean lesion lengths ranging from 6.02 to 8.92cm, and 25% of genotypes showed susceptibility, their mean lesion lengths varying from 9.02 to 27.2 cm ( Supplementary Figures 2B, E, F ).

We used GWAS analysis using 0.5 million SNPs from the selected accessions and their average scores against blast and BLB to unveil meaningful MTAs. Employing various multilocus GWAS methods (mrMLM, FASTmrMLM, and ISIS EM-BLASSO), we sought to identify significant associations by considering a LOD value ≥5 along with the phenotypic variance (PVE) as a threshold for significance. A total of 23 significant MTAs identified from the analysis, comprising nine for blast resistance and 14 for BLB resistance ( Table 1 ). The soundness of the model is evident from the well-fitted Manhattan and Q-Q plots ( Figure 1 ). The MTAs linked to blast resistance were dispersed across chromosomes 1, 4, 6, 7, 8, 11, and 12, with PVE ranging from 2.60% to 11.96%. Correspondingly, the BLB resistance revealed 14 MTAs, distributed across chromosomes 1, 2, 5, 6, 8, 11, and 12, encompassing PVE values ranging from 3.4% to 19.13%. These associations were identified through two distinct multi-locus GWAS models.

A comprehensive search yielded 107 candidate genes from the nine significant MTAs associated with blast resistance. Similarly, for BLB, 210 candidate genes emerged from the 14 significant MTAs. To enhance our understanding, all identified candidate genes for both blast and BLB were annotated using the RAP-DB database ( Supplementary Table 2 )

We compared the SNP positions of the prime MTAs pertaining to the target trait with those of established blast and BLB resistance genes or QTLs. By assessing candidate genes positioned within a range of −/+ 50kbp sequences from the peak SNP marker, we sought connections with known R-genes/QTLs. It’s noteworthy that all the significant MTAs for both blast and BLB were distinct, as there were no prior reports of R-genes or QTLs from these regions, except for the presence of the blast resistance gene Pi40 (16274830–17531111) in close proximity to MTA-17556914 on chromosome 6. Nonetheless, several other genes and transcription factors associated with defense response were successfully identified. For instance, among the 107 candidate genes related to blast resistance, 18 corresponded with already documented genes associated with defense response mechanisms. In contrast, for BLB, a substantial number of 54 candidate genes aligned with previously reported defense response genes were identified. However, when we expanded the search to include sequences within −/+ 100 kb from the significant MTAs, we found four overlapping resistant genes on chromosome 11 associated with bacterial blight: Xa23 (LOC_Os11g37620), Xa39 (LOC_Os11g37759), Xa10 (LOC_Os11g37570), and Xa46 (LOC_Os11g37540). These findings contribute to a deeper understanding of the genetic basis of blast and bacterial blight resistance, offering valuable insights for the development of improved disease-resistant varieties through targeted breeding programs.

All MTAs associated with blast and BLB resistance underwent a detailed exploration to uncover favourable SNP alleles. This investigation categorized SNPs demonstrating a negative SNP effect (ai) leading to increased resistance against blast and BLB diseases as “favourable alleles.” Conversely, alleles resulting in positive SNP effect (ai) to blast and BLB diseases were classified as “unfavorable alleles. Allele with more than 5% PVE with high phenotypic allele effect (ai) were selected as favourable allele. For blast and BLB traits, MTAs showcasing a reduction in phenotypic allele effect were recognized as harboring favorable alleles. For instance, for blast resistance, (MTA-193194162 and MTA-13878941 exhibiting high PVE with 11.5% and 11.96%, carrying ‘T’ and ‘T’ alleles denoting a decrease in phenotypic allele effect by -0.88 and -0.56 respectively. Similarly, for BLB, (MTA-335140024 and MTA-183761269), possessing the ‘A’ and ‘A’ alleles respectively, demonstrated a reduction in phenotype by -0.56 and -0.35 with high PVE of 13.38% and 19.13%. However, in case of blast and BLB, a total of four and nine MTAs were identified respectively as favourable alleles with high PVE. This study effectively illustrated that genotype possessing favourable SNPs for the target trait exhibited a substantial decrease in the phenotypic effect of the trait ( Table 1 ).

To determine the haplotype diversity, haplotype analysis was conducted on the 3K RG panel sub-set, correlating phenotypic outcomes against SNP data for all 317 identified candidate genes (107 for blast and 210 for BLB). This comprehensive analysis aimed to delineate the diversity of haplotypes within the 3K-RG panel. Significant haplotype diversity across 107 genes for blast and 210 genes for BLB were observed ( Supplementary Figure 5 ).

For a better comprehension of the phenotypic performance associated with specific haplotypes, the trait mean was calculated for each individual haplotype. The objective was to identify superior haplotypes, defined as those leading to significantly lower blast scores or lesion lengths for BLB compared to other haplotypes. Statistically validated by Duncan’s test, all candidate loci possessing more than two haplotypes were scrutinized to pinpoint superior haplotypes ( Table 2 ). A number of superior haplotypes linked to resistance traits were successfully detected. In the case of blast, eight superior haplotypes were defined by blast scores ≤2. Among the identified superior haplotypes, eight candidate genes exhibited blast scores ranging from 0.00 to 1.33. LOC_OS12G39700 and LOC_Os06g30440 stood out with superior haplotypes H4 and H33, showing the lowest average blast scores of 0.00-0.67 ( Figure 2 ; Supplementary Figure 3 ). Furthermore, the superior haplotype H29, associated with candidate gene LOC_Os01g18910, displayed an average blast score of 0.33-1.00. Similarly, H18 emerged as the superior haplotype for LOC_Os01g24540, LOC_Os01g24640, and LOC_Os01g24750, manifesting an average blast score within the range of 0.33-1.33. Additionally, candidate gene LOC_Os04g23960 featured the superior haplotype H17, while LOC_Os06g21040 harbored H177 and H60 as superior haplotypes with scores of 0.67-1.33.

For BLB, superior haplotypes were characterized by average lesion lengths ≤3. Out of the 210 candidate genes analyzed, four genes located on chromosomes 2, 5, and 6 yielded superior haplotypes ( Table 2 ). Among these genes, a total of five superior haplotypes with BLB scores spanning from 1.52cm to 4.86cm were identified. Notably, LOC_Os02g12660 contained the superior haplotype H39, showcasing the lowest average lesion length of 1.88 – 2.06cm ( Figure 3 ; Supplementary Figure 4 ). Identifying these superior haplotypes in the context of the studied traits holds promising implications for developing next-generation disease-resistant cultivars.

A selection of 10 accessions was successfully identified carrying superior haplotypes corresponding to eight candidate genes associated with blast resistance. These accessions fall into the aus and indx subpopulations across five distinct countries: Bangladesh, Brazil, India, Myanmar, and Pakistan. For instance, accession AUS 439:IRGC 29221-1 (LOC_Os01g24540-H18, LOC_Os01g24640-H18, and LOC_Os01g24750-H18) was recognized as bearing the superior haplotype for three candidate genes:. Similarly, two accessions, namely AUS 329:IRGC 29116-1 (LOC_Os01g24640-H18, LOC_Os01g24540-H18) and TAK::IRGC 73124-1 (LOC_Os06g30440-H33, and LOC_OS12G39700-H4) carried the superior haplotype of two candidate genes. The remaining accessions were carriers of superior haplotypes from a single gene ( Figure 4 , Table 3 ).

Eight accessions coming from six diverse countries (Bangladesh, China, India, Myanmar, Pakistan, and Sri Lanka) and spanning four subpopulations (aus, ind1A, ind2, and ind3) were identified as carrying superior haplotypes for BLB resistance. The average lesion length of these accessions ranged from 1.54 to 3.00 cm. Notably, accession MIN ZAO 6::IRGC 63772-1, possessing the superior haplotype of the candidate gene LOC_Os05g35490-H49, exhibited the lowest average lesion length of 1.54 cm. Three of the eight accessions—MIN ZAO 6::IRGC 63772-1 (1.54 cm), KUTTA::IRGC 52184-1 (1.70 cm), and KURULU WEE WHITE::IRGC 66518-1 (1.88 cm)—displayed the lowest average lesion lengths for BLB ( Figure 4 , Table 3 ).

The superior haplotypes identified through this study hold promising potential for the development of next-generation disease-resistant varieties via haplotype-based breeding approaches.