Abstract Background Rice blast caused by Magnaporthe oryzae is the most severe and devastating disease in rice results in serious losses worldwide. Based on this, the interaction between rice and M. oryzae has been studied extensively for decades, but the pathogen always has a negative effect on the new and emerging rice varieties. Results The present study employed comparative transcriptome strand-specific RNA sequencing and genome approaches of Diantun rice susceptible (D502) and resistance (D506) lines (leaves) in the presence of blast fungus, M. oryzae. Overall differential expression genes (DEGs) displayed 5838 and 3719 DEGs in D502 and D506, respectively 24hpi, however, the expression of DEGs in the former line was 5113, and in later line it was 4794 after 48hpi. Interestingly, only 2493 and 2418 DEGs were similar at both time hour points in both lines, respectively. Among DEGs, mostly exhibited down-regulated expression only in D502 major pathways, including plant hormones signal transduction and starch and sucrose metabolism at both time hours, suggesting susceptibility D502 on upon pathogen infection. Additionally, protein-protein interaction network analysis based on DEGs was performed between both varieties to find possible connections and increase interaction network complexity at 24h to 48h in D506, that might result in resistance to M. oryzae. We found many up and down-regulated DEGs only in D506 after pathogen infection, which might have a significant role in PTI and ETI immunity response. Next, through genomic analysis, different non-synonymous single nucleotide polymorphisms (nsSNPs) were identified between both D502 and D506 rice varieties. Here, four up-regulated genes, including WAK1, WAK4, WAK5, and OsDja9 harboring nsSNPs were found only in resistant D506 variety. Following alignment of open reading frame (ORF) region sequences revealed that the exonic SNPs lead the amino acid variation. Conclusion Our study proved that SNPs in these four genes were related to providing resistance in D506 line upon pathogen infection. In summary, we conclude that above-targeted rice defense and resistance genes identified through gene transcripts and modern genomic approaches could help us provide robust rice breeding and agricultural practices in future. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-025-06357-5. Keywords: Rice, Pathogen, Resistance, Defense, Transcriptomics, Genomics Background Rice (Oryza sativa) is a staple crop in Asian countries, serving as a significant source of calories for more than 50% of the population globally [[40]1]. Rice blast, caused by Magnaporthe oryzae, a hemibiotrophic fungus, is a highly destructive and damaging disease in rice at any stage and lead to more than 30% yield losses yearly worldwide [[41]2]. Due to its significant role in affecting rice production, the interaction with M. oryzae has been studied extensively for decades and has become a model crop system for plant-pathogen interactions [[42]3]. However, exploring novel and effective strategies for rice blast disease management is necessary. Over the past twenty years, comprehensive genetic and molecular research on plant-to-microbe interaction has reported that plants have developed two immunity stages. The first defense stage is regulated through identification of pathogen-associated molecular pattern (PAMP)-triggered immunity (PTI) by the extracellular domain of pattern recognition receptors (PRRs). In addition, host-derived elicitors known as damage associated molecular pattern (DAMPs) are commonly found in apoplast following pathogen attack [[43]4]. Furthermore, the second layer, known as effector-triggered immunity (ETI), that recognize and defends against different pathogens, leading to disease resistance [[44]4, [45]5]. This resistance is regulated by intracellular receptor molecules that possess nucleotide-binding (NB) and leucine-rich repeat (LRR) domains, which specifically target the effector proteins produced by M. oryzae [[46]6]. PTI and ETI initiate several similar immune responses such as; activation of mitogen activated protein kinase (MAPK), initiation of reactive oxygen species (ROS), temporary calcium influx, callose deposition, transcriptional changes, and regulation of plant hormones [[47]4, [48]7]. These responses, which lead to suppressed disease by inhibition of pathogen growth, that are successfully introduce virulence proteins referred to as effectors into the plant cell. Among them, mostly effectors seem to work by disrupting the events related to host immunity activation by recognizing PAMPs [[49]7]. In addition, wall associated kinase (WAK) are receptors like protein kinase, play an important role in plant immune response by acting as DAMPs and are also involved in many physiological processes, including plant growth, increased host immunity, and defense response against M. oryzae through regulation of the cell elongation [[50]8, [51]9]. It has been documented that WAK enhances resistance through positive or negative regulation toward different phytopathogens in many plants, including rice [[52]10, [53]11]. For instance, a study reported that three OsWAKs positive and one OsWAK negative regulator enhance rice blast disease resistance through eliciting ROS production and regulating the expression of defense genes [[54]11]. Therefore, WAK genes display several function diversities and have various mechanisms to enhance immunity in host plants by suppressing diseases caused by different pathogens. But, due to the genetic variability and pathogenic nature of M. oryzae, the resistance regularly decreases after a few years. In this regard, many approaches such as; RNA sequencing and microarrays between susceptible and resistant varieties, have been employed to examine the transcriptome profile of rice in response to M. oryzae [[55]12, [56]13]. Also, genomic analysis on the base of single nucleotide polymorphisms (SNPs) has been an emerging approach in recent decades to find out genetic variation and resistance mechanisms between susceptible and resistant varieties [[57]14]. Previously, transcriptome analysis allowed us to identify specific gene expressions between rice to M. oryzae interaction, which enhanced our understanding. For instance, a study found potential differentially expressed genes (DEGs) in susceptible and resistant varieties at different time hours exhibited up-regulated expression in defense, metabolism, plant-pathogen interaction, signaling and transport pathways of resistance variety after pathogen infection [[58]15]. Meanwhile, the DEGs between susceptible (parent line) LTH and their resistance monogenic lines were involved in signaling pathways [[59]16]. Moreover, to effectively identify genes that are responsible for phenotypic diversity, researchers commonly combine comparative genomic analysis with DEGs, function of homologous genes, genetic differences, and transgenic method studies for functional classification [[60]17]. The combined use of genomic and transcriptomic analysis could be an effective approach to identify resistance and defense-related genes in rice [[61]18]. However, to better understand rice immune responses against M. oryzae, it is important to perform an integrated analysis of transcriptome profile and SNPs variation of genes associated with pathogen signaling recognition and examine defense-related genes. However, Zhefang Diantun 502 (D502) rice, is one of the most important variety in Yunnan Province, China, and plays an essential role in rice breeding [[62]19]. Zhefang tribute rice D502 is a specialty of Zhefang town, Mang city, Yunnan Province, known for its bright color and good taste as well as considered as a long grain, soft, and not sticky rice. In 1623, King Duo Tian of Zhefang introduced its first time in Beijing. Later, it was considered as local rice, which was cultivated in China [[63]20]. In the last decades, the cultivation of this variety has been reduced due to the increasing population, but at the end of the economic era, it was cultivated again. Afterward, the Luxi Grain Bureau applied for a trademark of this variety, which was approved in 2004, and later, the cultivation and production of this variety significantly increased [[64]21]. In 2017, the Ministry of Agriculture observed D502 rice to register it as a Geographical Indication for Agricultural Products, leading to its recognition as a protected product under the National Geographical Indication system. It has been widely used for more than 30 years and is known for its high quality and good production of rice [[65]22]. However, the production is drastically reduced due to rising known and unknown fungal diseases. Developing new variety through breeding, which possess genetic resistance, is most effective and useful in mitigating diseases [[66]23]. Due to this, a new rice variety Diantun 506 (D506), was developed from Diantun 502, through rejuvenation research employing molecular methods, which later showed strong resistance to rice blast and has the potential to become a high-quality product in Yunnan Province, China [[67]24]. However, till now the mechanisms of resistance is unknown. In the present study, we reveal rice defense genes and pathways through analyzing gene transcripts and mining genome to explore the significant resistance genes based on nsSNP between susceptible Diantun 502 (D502) and its new development resistance line Diantun 506 (D506) upon M. oryzae. By analyzing different nsSNPs in the genome of both varieties, we successfully found four genes that showed high defense and resistant expression in the D506. Finally, the present work gives us enough information that the integration of advanced multi-omics approaches could be used to identify defense or resistance genes in crop plants. Materials and methods Plant material and seed growing condition The new rice variety Diantun 506 (D506), developed through rejuvenation research employing molecular methods, exhibits strong resistance to rice blast and has the potential to become a high-quality product in Yunnan Province, China. Therefore, the present study used Diantun rice (D502 parent line) susceptible variety of Yunnan province, and (D506), a resistance line for rice blast resistance evaluation. The compatible rice line D502 and its new developed line D506 were planted in a similar conditions according to previous study [[68]16]. Rice blast fungal strains and inoculation A total of 13 strains of Magnaporthe oryzae, such as I2, I5, I6, I8, YN125, YN132, YN199, YN358, YN396, YN541, YN626, Guy11, and 617 were used for inoculation into D502 and D506, respectively. The conidial suspensions 1×10^5 conidia mL^-1 with 0.02% Tween-20 of tested strains were inoculated by spraying method on 20 days old-grown D502 and D506 rice seedlings and distilled sterilized water (ddH[2]O) were sprayed as the negative control (mock) plants of D502. A total of seven rice plants were used in each replicate, and three replicates were used in each treatment. The mock and fungus-inoculated rice seedlings were placed in dark chambers with 100% humidity at 25℃ and then transferred to the greenhouse after 24 hours. The disease score was recorded at 7dpi based on the severity of lesions [[69]25]. Gene transcripts analysis of tested rice varieties upon pathogen infection The fully expanded second and third leaves of 7 rice seedlings from each group were harvested and immediately frozen in liquid nitrogen for RNA-seq at 24hpi and 48hpi, with three biological replicates in each treatment. Total RNA was extracted from rice leaves using a plant RNA kit (OMEGA, USA) and degradation and contamination was observed on 1% agarose gels, and then evaluated for purity using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using the Qubit® RNA Assay Kit in a Qubit® 2.0 Fluorimeter (Life Technologies, CA, USA), and the integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). For library preparation, 1 µg of total RNA per sample was processed as input material for the RNA sample preparations, and libraries were created using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer's protocol. mRNA was enriched using magnetic beads conjugated with poly-T oligos. Fragmentation and priming of mRNA were performed using NEBNext First Strand Synthesis Reaction Buffer (5X), and first-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). Second-strand cDNA synthesis was carried out with DNA Polymerase I and RNase H, followed by adenylation of 3’ ends and ligation of NEBNext Adaptors. AMPure XP beads (Beckman Coulter, Beverly, USA) were used to purify cDNA fragments ranging from 250-300 bp. After USER enzyme (NEB, USA) treatment at 37°C for 15 minutes, followed by heat inactivation at 95 °C for 5 minutes. PCR amplification was performed using Phusion High-Fidelity DNA polymerase, Universal primers and Index primers. Final libraries were purified with AMPure XP beads and validated using the Bioanalyzer 2100 system before sequencing on the Illumina Nova-seq 6000 platform, generating 150 bp paired-end reads. Raw sequencing data were processed to remove low-quality reads, and Q20, Q30, and GC content metrics were calculated. Clean reads were mapped to the reference genome using HISAT2 v2.0.5, and gene counts were quantified with feature Counts v1.5.0-p3. The sequencing data has been submitted to the NCBI under BioProject number PRJNA1205520 available at ([70]https://www.ncbi.nlm.nih.gov/search/all/?term=PRJNA1205520). Analysis of differentially expression genes The DESeq2 package (1.42.0) in R software was used to identify differentially expressed genes (DEGs) in digital gene expression data through models according to binomial distribution in both varieties. The DEGs were successfully assigned with a log2FC value > 1 and a p-value ≤ 0.05 as differentially expressed between mock- and fungus-treatment in D502 and D506, respectively. EdgeR was used with single-scale normalization factors to perform differential expression analysis in cases without biological replicates. A TB tool II was used to show the expression of DEGs at different time hours. Gene ontology and KEGG pathway analysis The Gene ontology and KEGG pathway of the DEGs were performed using clusterProfiler (4.8.1) R package, which included a correction for gene length bias. Moreover, a significant p-value ≤ 0.05 GO terms, including Biological process (BP), Cellular component (CC), and Molecular functions (MF) were identified in both varieties. KEGG (Kyoto Encyclopedia of Genes and Genomes) assists a complete database, which is designed to elucidate the advanced functions of biological systems. However, the KEGG database ([71]http://www.genome.jp/kegg/) was used to determine the related pathways and functions of DEGs. Furthermore, clusterProfiler R package was used to test the statistical enrichment of DEGs in KEGG pathways. Gene expression profile The maSigPro package in R software was used to explore the expression profile of genes resulting from normalized time course RNA-seq data (log2 FPKM +1) between the susceptible D502 and resistant D506 rice species. The genes with significantly different expressions were visualized by plotting their expression profiles and characterizing them into distinct clusters. After that, the clusters three genes, that showed high expression in D506 were used for GO enrichment and KEGG pathways by using DAVID 6.8 [72]https://david.ncifcrf.gov/ online software. Protein-protein interaction network analysis To explore interaction between DEGs, we carried out protein-protein interaction network analysis through STRING database. For species listed in the database, we build the network by retrieving the target gene list from the available database. The resulting network was made up of nodes (protein) and edges (interaction), illustrating the relationships among proteins. The PPI networks were visualized by using Gephi software (V.0.10), and the hub genes were selected based on topological analysis including degree, closeness and betweenness centrality. In addition, the DEGs from each module were utilized for GO enrichment and KEGG pathways analysis by using DAVID 6.8 tool ([73]https://david.ncifcrf.gov/) and significantly identified with p-value < 0.05. Genome sequencing Rice genomic DNA samples were extracted, and their quality was assessed. A total of 0.2μg genomic DNA per sample (3 biological replicates) was performed using CTAB method for DNA library preparations. Furthermore, sequencing libraries were prepared using NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, USA). The genomic DNA samples were fragmented by sonication to a size of 350 bp. After that, the fragment underwent end-polished, A-tailed, and was ligated with a full-length adapter for Illumina sequencing with further PCR amplification. After purifying the PCR products through (AMPure XP system Beverly, USA), the libraries qualities were evaluated for size distribution by (Agilent 5400 system, Agilent, USA) and quantified by QPCR (1.5 nM). The qualified libraries were combined and sequenced on Illumina platform with PE150 strategy, based on the optimal library concentration and the necessary data volume. Furthermore, the original sequencing data was obtained from Illumina platform and was transformed in (raw read or raw data) by base calling and recorded in FASTQ format, and sequence adapter and low-quality data were cleaned. Burrows wheeler aligner (BWA) software was used for mapped the sequence read to the reference genome and original results were obtained in BAM format file. Then the duplication was removed by SAMtools and Picard ([74]http://broadinstitute.github.io/picard/) software. SNP calling GATK2 (v3.7) software was used to perform SNP calling. Raw vcf files were filtered with GATK standard filter method and other parameters (cluster:3; WindowSize:35; – gcpHMM 10 -stand_emit_conf 10 -stand_call_conf 30, and QD < 2.0, FS > 60.0, DP < 10, MQ < 30.0, HaplotypeScore > 13.0). The SnpEff software was used to annotation for the Variable site. CNVs were identified using CNVnator to determine possible deletions and duplications (parameter: -call 100). Structural variants (SVs) were identified with BreakDancer and the functional annotation of variants was identified by ANNOVAR. The known genes from UCSC were utilized for annotating genes and regions. Genomic analysis Genomic analysis was performed to find out the genes carrying different non-synonymous single nucleotide polymorphisms (nsSNPs) in both D502 and D506 varieties. A Venny 2.2 online tool ([75]https://bioinfogp.cnb.csic.es/tools/venny/) was used for identify the different SNPs between the susceptible and resistance varieties. The genes harboring different nsSNPs were subjected to the function enrichment analysis by using DAVID 6.8. The significantly GO and KEGG pathways were presented in a graph. Sequence alignment and 3D Protein structure analysis The NCBI open reading frame (ORF) finder ([76]https://www.ncbi.nlm.nih.gov/orffinder/) was used for the ORF region of four candidate genes. The sequences of four genes were aligned by MEGA X and Bioedit software. The protein sequence analyses were conducted by PSIPRED v4.0 ([77]http://bioinf.cs.ucl.ac.uk/psipred/). The Alphafold3 online available at ([78]https://alphafoldserver.com/), was used for three-dimensional 3D prediction of protein structures of selected genes, and the pictures of the structures was generated by PyMOL. Statistical analysis The experiment was carried out in the greenhouse conditions with six treatments and each treatment have three replicates. Statistical parameters such as mean, standard deviation, analysis of variance, LSD multiple comparison tests were calculated by using IBM SPSS statistics 27. The GraphPad prism 8 version was used to create the graph, TB tool II utilized for generating the heatmap of DEGs, and the graph and heatmaps were edited with the help of Adobe illustrator CC 2019. Results Response of Diantun 502 and Diantun 506 seedlings to Magnaporthe oryzae Diantun 502 (D502) is highly susceptible to almost all rice blast fungal strains, which have been used as a recurrent parent to develop new resistance variety such as Diantun 506 (D506) against rice blast disease. The genetic relationship between D502 and D506 was very close. Therefore, in this study, the susceptible (parent) D502 and their resistance line D506 were used to verify the rice blast resistance. The disease score was recorded after 7 days of inoculation with 13 fungal strains, and it was found that the D502 was very susceptible to all 13 strains of blast fungus compared to D506. The disease lesions on D502 ranged from 3-5, whereas only one strain (YN541) showed disease lesions ranging from 3-4 on D506 (Fig. [79]1A). The disease index (DI) of susceptible and resistance variety was recorded after inoculation with 13 fungal strains. It was found that the DI of D502 was higher, ranging from 25.05 -78.88, while the DI of D506 was found lower between 4.33-30, whereas only one strain (YN541) showed 61.04 disease index on D506, which was also lower than D502 (Fig. [80]1B). These results revealed that D506 showed high resistance to almost all strains of blast fungus as compared to D502. Fig. 1. [81]Fig. 1 [82]Open in a new tab Response of Diantun rice seedlings toward 13 different strains of Magnaporthe oryzae. A Disease symptoms of the pathogenic strains on leaves of D502 and D506 after inoculation in greenhouse. B Disease index representing the severity of the disease on the rice leaves. The data are presented in mean standard deviation (M ± SD, n> 21 leaves). Error bars represents ± SD and different letters showed differ significant (p-value ≤ 0.05) High throughput RNA-seq analysis for differentially expressed genes (DEGs) To obtain the transcripts from D502 and D506 inoculated leaves with 13 strains of Magnaporthe oryzae, the high-throughput RNA-Seq was carried out. More than 803 million raw reads were generated from 18 libraries, and 97.42% were clean. The 96.82 and 92.14 averages were calculated of Q20 and Q30 of these data, respectively. All clean reads were linked with the rice genome, about 92.99% could be mapped to the rice reference genome of the Oryza sativa Indica group (ASM465v1), and 89.48% of reads were uniquely mapped (Table [83]1). The results indicated that these data were highly dependable, and expression levels of the obtained transcript could be utilized for further analysis. Table 1. Statistical data of the RNA-Seq reads for samples Samples Raw reads Clean reads Raw bases(G) Clean bases(G) Q20% Q30% Mapped Reads% Concordant Pairs% D502_0_1 46102338 44813312 6.92 6.72 96.8 92.18 92.95 89.41 D502_0_2 43352018 42533016 6.5 6.38 95.48 89.07 93.2 89.62 D502_0_3 43566520 42625160 6.53 6.39 96.92 92.31 92.92 89.17 D502_24_1 40673120 39585930 6.1 5.94 96.86 92.18 93.1 89.68 D502_24_2 47770862 46542680 7.17 6.98 96.99 92.54 93.07 89.67 D502_24_3 45423032 43968984 6.81 6.6 97 92.51 92.95 89.39 D502_48_1 42082688 40732666 6.31 6.11 97.02 92.59 93.25 89.67 D502_48_2 44774154 43709146 6.72 6.56 97.01 92.6 92.8 89.68 D502_48_3 42944632 42495160 6.44 6.37 96.85 92.15 93.02 89.61 D506_0_1 45954468 45188236 6.89 6.78 96.77 91.99 93.82 90.44 D506_0_2 46344380 45048506 6.95 6.76 96.79 92.17 92.9 89.17 D506_0_3 50470710 49288016 7.57 7.39 96.64 91.76 93.48 89.75 D506_24_1 44624476 43393764 6.69 6.51 96.81 91.98 91.71 88.57 D506_24_2 41516852 39920676 6.23 5.99 97.2 92.93 92.6 89.02 D506_24_3 44577642 43384212 6.69 6.51 96.94 92.4 93.08 89.58 D506_48_1 44491444 43271422 6.67 6.49 96.96 92.49 92.84 89.18 D506_48_2 46779556 45550552 7.02 6.83 96.94 92.39 93 89.4 D506_48_3 42281344 41019046 6.34 6.15 96.95 92.42 93.27 89.75 AVERAGE % 44651679.78 43503915.78 6.69 6.52 96.82 92.14 92.99 89.48 TOTAL 803730236 783070484 120.55 117.46 [84]Open in a new tab To investigate the comparison and grouping among all samples, we performed principal component analysis (PCA) based on the fragments per kilobase of transcript per million mapped reads (FPKM) value of all obtained genes. The samples were organized according to PCA1, which accounts for (32.19%), and PCA2, which describes (17.65%) of the variances. In addition, the mock samples of both varieties were clustered nearly in one plot of PCA, indicating the samples are more similar. In contrast, the treatment samples of D502 showed high similarity and clustered together in plot 1, whereas the samples of D506 were clustered closely in plot 3 in terms of PCA analysis and then separated based on time points after M. oryzae infection, respectively (Fig. [85]2A). Hierarchical clustering analysis was carried out according to the FPKM value of all identified genes, which revealed that the mock and treated samples of both varieties were clustered into seven groups. Among them, the mock samples of both varieties were clustered into two different groups, while the biological replicate of D502 at 24h and 48h were clustered into three groups (Fig. [86]2B), suggesting that M. oryzae continuously induced highly transcriptional changes in the susceptible D502. Interestingly, the biological replicate samples of D506 at 24h were clustered into two groups, but in the term at 48h, the samples were clustered into one group (Fig. [87]2B), indicating that minimum transcriptional changes after M. oryzae infection in D506. Fig. 2. [88]Fig. 2 [89]Open in a new tab Detail of biological replicate samples analysis and transcripts of differentially expressed genes (DEGs) between Diantun rice D502 and D506 in response to Magnaporthe oryzae at different time hours interval. A Principal component analysis (PCA) of treated and mock samples data at 24h and 48h in D502 and D506. The treated samples were inoculated with M. oryzae and mock ddH[2]O water. B Hierarchical clustering analysis of 18 samples based on FPKM value of all identified DEGs. C Expressed DEGs in between D502 and D506 varieties of rice. D Venn diagram showing core and unique expressed DEGs at both time interval. E The up-regulated and down-regulated expressed DEGs in D502 and D506 at different time. Function annotation analysis of DEGs between susceptible and resistance variety To better understand the susceptible and resistant mechanisms in D502 and D506, transcriptome sequencing was carried out from the leaves after infection of Magnaporthe oryzae at 24hpi and 48hpi. The differential expression genes (DEGs) were identified on the base of log2FC >1 and -value ≤ 0.05. Therefore, a total of 5838 DEGs (2975 up-regulated and 2863 down-regulated) and 3719 (1612 up-regulated and 2107 down-regulated) were identified at 24hpi, while 5113 (2605 up-regulated and 2508 down-regulated) and 4794 (1753 up-regulated and 3041 down-regulated) were found at 48hpi in D502 and D506, respectively (Fig. [90]2C). Among identified DEGs, only 2493 (35.3%) and 2418 (32.3%) were similar at both time hour points in D502 and D506 (Fig. [91]2D). Out of these, only 587 (12.8%) up-regulated and 830 (16.1%) down-regulated of the transcripts were commonly regulated at all time points in both variety D502 and D506, representing the genetic similarity based on their expression profiles (Fig. [92]2E). The obtained DEGs in D502 and D506 were used for further GO and KEGG enrichment analysis. GO and KEGG pathway analysis The GO enrichment and KEGG pathways analysis of up-regulated and down-regulated DEGs were conducted between D502 and D506 at different time hours. Biological process (BP) of up-regulated DEGs, including ribosome biogenesis (GO:0042254) and ribonucleoprotein complex biogenesis (GO:0022613), were significantly enriched in resistant variety at both time hours, but in susceptible only at 24h. In addition, the nucleoside metabolic process (GO:0009116) and glycosyl compound metabolic process (GO:1901657) were significantly enriched in the susceptible at only 24hpi. In contrast, while rRNA processing (GO:0006364), rRNA metabolic process (GO:0016072), and ncRNA processing (GO:0034470) were highly associated with resistance at 24h and 48h. Additionally, preribosome (GO:0030684) and mitochondrial part (GO:0044429) in the cellular process were significantly enriched in D506 at both time hours, but in D502 at only 24h. The ribonucleoprotein complex (GO:1990904) and ribosome (GO:0005840) were present in susceptible and resistant at 24hpi. At the same time, the up-regulated DEGs were significantly involved in the term of mitochondrion (GO:0005739) at 24h only in resistant variety. Moreover, the up-regulated DEGs were significantly regulated in only one iron ion binding molecular function at all time hours in both cultivars, while RNA binding (GO:0003723), structural molecule activity (GO:0005198), pseudouridine synthase activity (GO:0009982) and serine-type peptidase activity (GO:0008236) were present only in resistant at different time hours. Furthermore, the KEGG result revealed that only one similar up-regulated pathway, including ribosome biogenesis in eukaryotes (osa03008) was highly enriched in D506 at both time hours, then in D502 at 24h, while the DEGs were significantly regulated in term of ribosome only at 24h in both varieties (Fig. [93]3A). Fig. 3. [94]Fig. 3 [95]Open in a new tab GO and KEGG pathways enrichment of expressed DEGs between both varieties. A Up-regulation of DEGs in GO and KEGG pathways; (B) down-regulation of GO and KEGG pathways in susceptible and resistant varieties of rice. (C and D) Heatmap displaying the expression of up and down-regulated DEGs in five different GO and KEGG pathways. The function annotation of DEGs was analyzed and revealed that the DEGs were down-regulated in three similar biological processes such as; drug metabolic process, pyridine nucleotide metabolic process and nicotinamide nucleotide metabolic process at 48h in both varieties, while no DEGs were significantly regulated in cellular component process. Furthermore, the DEGs were down-regulated at all time hours in two similar molecular functions including, DNA-binding transcription factor activity and transcription regulator activity in both D502 and D506 after blast pathogen infection. However, the KEGG pathways results revealed that most DEGs were significantly down-regulated in D502 at both time hours, while minor in D506 at only 48h. The DEGs were enriched in the term of MAPK signalling pathway-plant, and plant-pathogen interaction were down-regulated in D502 at both time hours than in D506 at only 48h, while in the plant hormones signal transductions and starch and sucrose metabolism pathways were down-regulated only in susceptible at both time hours. Minimum DEGs were down-regulated in five KEGG pathways, including Glycolysis / Gluconeogenesis, pentose phosphate pathways, carbon fixation in photosynthetic organisms, carbon metabolism, and biosynthesis of amino acids at 48h in D506 (Fig. [96]3B). Taken together, these results suggest that when the rice attack to blast fungus, many DEGs were showed down-regulation in response to M. oryzae infection in D502 at both time hours points, but in minor and incompatible in D506 at 48h. Expression of Up and down-regulated DEGs in KEGG pathways To explore the expression of DEGs in both up and down-regulated KEGG pathways, initially analyzed the ribosome, ribosome biogenesis in eukaryotes, plant hormone signal transduction, MAPK signaling pathway-plant, starch and sucrose metabolism and plant-pathogen interaction involving genes. A total of 236 genes were found in these six pathways at both time hours, of which more than 34 genes in ribosome and ribosome biogenesis in eukaryotes were found up-regulated in D506 (Fig. [97]3C). Among down-regulated DEGs, many genes were showed highly downstream expression in four down-regulated pathways only in D502 at different time hours points, these including OsEIL4, OsGH3-13/OsTLD1, OsSMG, OsRR26, OsPYL11, OsbHLH105, OsPP2C49, OSERS2, OsETR2, OsIAA5, OsHP3, OsSAURs, OsNPR3, OsMKKs, OsRTH2, OsCDPs, OsRLCK123, OsKCS3, OsSAPKs, OsFLS2, OsCERK1 and OsCEBiP (Additional file 2; Table [98]S1). Out of 236 DEGs, 23 were found down-regulated in plant hormone signal transduction pathways only in susceptible at 24h and 48h, which belong to the defense-related hormones, disease resistance proteins, and TF families such as; indole-3-acetic acid (IAA), mitogen-activated protein kinase (MAPK), Cytokinin, abscisic acid (ABA), bHLH, Ethylene receptor (ETR), NPR1 Disease resistance (PR), Auxin response and Serine/threonine protein kinase. Consequently, the starch and sucrose metabolism pathways displayed more than 20 down-regulated DEGs, related to the starch biosynthesis and starch metabolism (Additional file 2; Table [99]S1). Furthermore, a total of 39 DEGs were involve in MAPK signaling pathway-plant and half of them showed highly down-regulated expression in D502 including, OsMKK4, OsPYL11, OsRTH2, OsPP76, OS-ERS2, OsEIL4, OsMKK5, OsbHLH105, OsSAPK2, OsFLS2 and OsSAPK7 at various time hours. In addition, 12 DEGs were found in plant-pathogen interaction pathways, which were highly expressed down-regulated upon M. oryzae infection in susceptible. These included OsCDPK21, OsSMG1/OsMEK6, Pti1, OsCDPK12, OsMKK5, OsKCS3, OsCERK1, and OsFLS2, whereas OsCEBiP (chitin elicitor binding protein precursor) which is a key regulator of OsMYB22 defense response and its relating protein were found up-regulated in resistance at 24h and down-regulated in susceptible at 48h, respectively (Fig. [100]3D). These results suggest that, the down-regulation genes including OsCEBiP in these three pathways play a significant role in the susceptible and could involve in the weekend host defense response against M. oryzae. Moreover, we also analyzed phenylpropanoid biosynthesis, MAPK signaling pathway-plant, and plant-pathogen interaction pathways in resistance variety and didn’t find any important defense-related genes. Expression profile of DEGs in D502 and D506 at different time course To check the different expression of genes among both rice varieties, we conducted mRNA-seq at 24- and 48-time hours after inoculation by M. oryzae. The maSigpro package was used to identify four clusters of gene with similar expression. A total of 10,002 genes were identified among four cluster groups. However, the cluster 3 genes showed high expression in D506 compared to D502 (Fig. [101]4A). Meanwhile, the other three cluster genes expression are reversed and noncompatible (Additional file 1; Fig. [102]S1 and Additional file 2; Table [103]S2). Therefore, cluster 3 genes were used for GO enrichment and KEGG pathways by using DAVID bioinformatics online tool. The results revealed that sixteen genes were significantly involved in term of biological process plant defense response such as; defense response to other organism (GO:0098542), response to other organism (GO:0051707) and defense response (GO:0006952). In addition, nine genes participated in cellular components, and six enriched in molecular functions, respectively (Fig. [104]4B). In the KEGG pathways analysis, 11 genes were significantly enriched in terms of metabolic pathway, 3 associated with cysteine and methionine metabolism and 2 involved in biosynthesis of unsaturated fatty acids, respectively. After that, to find out more information and understand the role of expressed genes after M. oryzae infection, the expression pattern of cluster 3 gene was analyzed, which were regulated in GO and KEGG pathways, respectively (Additional file 2; Table S3). Fig. 4. [105]Fig. 4 [106]Open in a new tab Differentially expressed genes in the cluster profile of Diantun susceptible and resistance varieties of rice. A Expression of profile cluster three were determined from time-course RNA-Seq data of leaf samples at 24h and 48h post-inoculation with Magnaporthe oryzae. B Go and KEGG pathways enrichment analysis of cluster three DEGs by using DAVID bioinformatics online tool. Protein-protein interaction network analysis of DEGs To examine the interaction between DEGs, we performed a PPI network analysis of susceptible and resistant varieties of DEGs. A Gephi software (V. 0.10) was used to visualize the PPI interaction network of DEGs. A total of 4 modules were selected, including orange, purple, dark green, and blue, which showed similar DEGs. In D502, the network interaction complexity was higher at 24h but decreased at 48h in all modules. Meanwhile, at 24h in D506, the interaction network complexity was also higher and stable or increased at 48h in all modules after M. oryzae infection. Therefore, the maximum DEGs were observed in orange, purple, dark green and blue module at both time hours in D506. While in D502, the orange, purple, dark green, and blue modules contained high DEGs at 24h, but the number of DEGs were continuously decreased in similar modules at 48h (Fig. [107]5A, Additional file 2; Table S4). The highest 1833 nodes and 11966 edges were identified at 24h, and 1404 nodes and 4388 edges were found at 48h in D502. Subsequently, the 1046 nodes and 5740 edges were observed at 24h and 1558 nodes and 6568 edges at 48h in D506 (Fig. [108]5B). The maximum average degree and minimum modularity were found at 24h in both varieties, but the modularity was significantly increased at 48h in D502 and D506 (Fig. [109]5C). These results suggested that the loose connection and low interaction network complexity of DEGs between 24h to 48h could lead to susceptibility after M. oryzae infection in D502. Furthermore, the functional enrichment analysis was performed using DEGs from each module in both varieties at different time hours. Fig. 5. [110]Fig. 5 [111]Open in a new tab Protein-protein interaction network and functional enrichment analysis of DEGs between susceptible and resistance verities. A Interaction network of DEGs in different modulaes at different time hours in both varieties. B Nodes and edges; (C) average degree and modularity; (D and E) GO enrichment and KEGG pathways analysis of four module DEGs in D502 and D506 at 24h and 48h. (A, D and E) Different colours indicate the module number between both varieties. A total of 19, 18, 16, and 16 GO and KEGG pathways were significantly regulated at 24h, and 20, 11, 8, and 10 pathways were found at 48h in different modules of D502. Whereas 20, 14, 19, and 9 pathways were identified at 24h, and 20, 16, 20, and 14 GO, and KEGG pathways were significantly found at 48h in various modules of D506 (Additional file 1; Fig. [112]S2). Among them, mainly DEGs were involved in defense-related and enzymatic GO and KEGG pathways in different modules of D506, such as cell wall macromolecule catabolic process, chitin catabolic process, abscisic acid biosynthetic process, chitinase activity, chitin binding, ethylene-activated signaling pathway, regulation of ethylene-activated signaling pathway, MAPK signaling pathway–plant, and plant hormones signal transductions at 24h, and tryptophan biosynthetic process, L-phenylalanine catabolic process, kinase activity, and phenylalanine, tyrosine and tryptophan biosynthesis at 48h (Additional file 1; Fig. [113]S2). It was found that many genes were significantly regulated in more than 20 similar GO and KEGG pathways in D502 at both time hours, but the gene numbers were highly decreased in some pathways at 48 hours, respectively (Fig. [114]5D and Additional file 2; Table S5). While, the DEGs in D506 were significantly involved in more than 30 similar GO and KEGG pathways at both time hours, and the gene number highly increased in many pathways at 24h to 48h upon pathogen infection (Fig. [115]5E and Additional file 2; Table S5). Meanwhile, it is suggested that the increasing number of DEGs after M. oryzae infection, which was regulated in defence and immunity-related pathways at different time hours, might show resistance in D506. DEGs involve in PTI and ETI signalling pathways To combat pathogen attacks, plants have developed two innate immunity systems, referred to as PTI and ETI, which are activated by cell surface PRRs or NLR intracellular receptors [[116]26]. Therefore, to better understand the functions of up and down-regulated DEGs and their role in the activation of PTI and ETI against M. oryzae were checked. Some DEGs were involved in enzyme activity encoding P450 cytochrome, and transcriptional factors (TFs) such as AP2/ERF, bZIP, MYB, NAM, and WRKY in D506, but not in D502 (Fig. [117]6A). TFs and P450 genes play an important role in plant immunity, and also expressed various defense genes against M. oryzae infection in rice [[118]16]. The functions of some DEGs in defense response were unknown, which were designated as others. Disease resistance (R) proteins play a significant role in activation of plants immunity. Mostly, R protein encoded a nucleotide-binding domain (NB-ARC), which is involved in the activation of ETI upon pathogen infection [[119]27]. There were thirteen up-regulated, and six down-regulated NB-ARC, and NBS-LRR DEGs were found at different time hours only in the resistant variety. These DEGs are suggested to activate ETI in response to M. oryzae infection in D506. The cell surface PRRs are commonly related to receptors like kinase or protein (RLKs and RLPs), where they understand P/M/DAMP to trigger the PTI signaling pathways [[120]3]. However, three LRR-RLPs, two LRR-RLKs, Six RLCK-RLK, three WAK-RLK including WAK1, WAK4, and WAK5 and two SDR-RLK were found up-regulated. In comparison, two LRR-RLPs, one LRR-RLKs, four RLCK-RLK and one WAK-RLK were expressed down-regulated upon M. oryzae infection at 24h and 48h only in D506 (Fig. [121]6A). Moreover, one respiratory burst oxidase homolog (RBOH) and one mitogen-activated protein kinase (MAPK) regulate upstream, and at the same time, four DEGs in MAPK were identified as down-regulated upon M. oryzae infection in D506 at different time hours. In addition, sixteen up-regulated and eleven down-regulated defense-related TFs were found at different time hours in D506 encoding AP2/ERF, bZIP, MYB, NAM and WRKY, which are associated with different defense-related hormones, including salicylic acid (SA), jasmonic acid (JA), abscisic acid (ABA), gibberellic acid (GA) and ethylene (ET), respectively (Fig. [122]6A). Calcium plays a crucial role against various stress of plants. Two up-regulated and one down-regulated ca^+ DEGs were found, which were associated with two pore calcium channel proteins (TCP), calmodulin (CAM) and cyclic nucleotide-gated ion channel (CNG), which suggested that the DEGs were also involved in calcium influx. Interestingly, one disease defence-related TF OsbHLH65, and two genes OsRNG3, OsSDI1 were up-regulated in D506, but down-regulated in D502 after infection of pathogen (Fig. [123]6B). Taken together, the up and down-regulation of defense genes and TFs, which associated with defense-related plant hormones could also affect the PTI and ETI activation in resistance variety as compared to susceptible after M. oryzae infection. Consequently, a total of 50 (60.97%) up-regulated defense-related genes were present only in D506. Among them, 13 were up-regulated at 24h and 31 at 48h, while 6 showed similar expression at both time hours (Fig. [124]6C), which suggested that continued expression could also be essential in response to pathogen infection. Fig. 6. [125]Fig. 6 [126]Open in a new tab Expression of up and down-regulated DEGs invloved in increase plant defense and immunity signalling pathways in resistant variety. A Up and down-regulated expression log2FC value of DEGs in response of M. oryzae infection. B Three DEGs exhibit opposite expression between both varities at different time hours. C Venn diagram of only up-regulated DEGs in (A) at both time hours in D506 Genome analysis between susceptible and resistance varieties To explore the genetic relationships between D502 and D506, we sequenced the genome of both varieties. The raw reads were clean, and 276.07 and 260.1 million reads were mapped from total reads with 98.08 and 98.03 mapping rate % to the rice reference genome of the Oryza sativa Indica group (ASM465v1). The mean genome coverage was higher in the D506 and lowest in D502, while the average depth and coverage percentage were found mostly similar in both varieties (Additional file 2; Table S6). Additionally, more than 67% of SNPs were present in intergenic region of both varieties, and the rest are associated with genic region. Among the genic region, many SNPs were associated with intronic region (Additional file 2; Table S7). In addition, the SNP mutation type distribution, and frequency were higher of all six combinations in D502, compared to D506 (Additional file 1; Fig. S3A Additional file 2; Table S8), which were found in CDS region, and divided into two categories such as; synonymous and non-synonymous. In comparison to synonymous, the non-synonymous SNPs (nsSNPs) were higher in both varieties. According to their functional annotation, all SNPs were classified into four groups such as; high, low, moderate and modifier impact. Among these groups, the modifier impact contains a maximum of 56.98% SNPs in D506, then 45.49% in D502. This was followed by moderate, low, and high impact in both varieties, respectively (Additional file 1; Fig. S3B and Additional file 2; Table S9). Therefore, the minimum number of high impact SNPs were detected, that were disturbing the transcript sites (splice site acceptor and donor), and stop codon (stop site gain and loss). Identification of genes harboring nsSNPs and their functional analysis Based on nsSNPs analysis, we identified 2668 and 1247 genes carrying different nsSNPs in susceptible D502 and resistance D506. The annotation analysis from both variety genomes revealed that a large number of SNPs were present in the intergenic, intronic, and exonic regions of these genes (Fig. [127]7A). These genes were used for functional enrichment analysis using the DAVID bioinformatics online tool. The enrichment analysis revealed that the genes carrying nsSNPs, mostly down-regulated in D502, were significantly associated with defense and immunity-related GO enrichment and KEGG pathways compared to D506 variety. Remarkably, the DEGs were mainly involved in phosphorylation, protein phosphorylation, response to other organisms, and diterpenoid metabolic process in D502, harboring nsSNPs in different regions. Moreover, mostly genes were significantly regulated in nucleus cellular component groups, ATP binding, protein binding, and protein serine/threonine kinase activity molecular functions of GO terms. The KEGG pathway results revealed that mostly DEGs were significantly involved in the biosynthesis of secondary metabolites followed by DNA replication and tropane, piperidine and pyridine alkaloid biosynthesis, carried nsSNPs, and exhibited down-regulated expression in D502, respectively (Fig. [128]7B). In comparison, the minimum genes were regulated in the BP, CC, MF, and KEGG pathway group of D506 (Fig. [129]7C). These results suggested that the maximum number of genes carrying nsSNPs present in GO enrichment, and KEGG pathways in D502 are likely associated with susceptibility to M. oryzae infection compared to D506. Fig. 7. [130]Fig. 7 [131]Open in a new tab Identification of genes harboring non-synonymous single nucleotide polymorphisms in D502 and D506 and their functional enrichment. A Overall number of genes carrying nsSNPS in different regions of both varieties. B and C GO and KEGG pathways enrichment analysis of genes, unique in both D502 and D506, classified into three categorizes such as; Bilogical Process (BP), Cellular Component (CC), and Molecular Function (MF), and (D) selected defense genes expression in the resistant variety D506 Combined transcriptome and genome analysis revealed resistance in Diantun 506 Combined transcriptome and genome analysis were conducted to understand the effects of nsSNPs on susceptible and resistance and to investigate the expression of genes harboring different nsSNPs in both varieties. A 37 up-regulated and 17 down-regulated genes were identified exclusively in D506, carrying nsSNPs at different positions in D502 and D506. These including 8 NB-ARC, 5 LRR, 10 RLCK, 4 WAK, 10 P450 cytochrome, 1 squamosa promoter-binding-like protein (SPL16), 5 Putative pentatricopeptide repeat-containing protein (PPR), 1 OsDnaJ, 5 transcriptional factors and 3 E3 ubiquitin protein ligase genes (Additional file 1; Fig. S4). From them, four candidate genes, viz., WAK1, WAK4, WAK5 and OsDja9, which are associated with disease resistance and immune responses, were found to be up-regulated only in D506. The WAK1 and WAK4 genes exhibited more than 5-6-fold change up-regulated expression in D506 at 24h upon pathogen infection. Whereas WAK5 showed a more than 1.5-fold change increase in expression, while OsDja9 expressed more than 1-fold change at 48h only in D506 (Fig. [132]7D). Furthermore, the open reading frame (ORF) region sequences of four genes were aligned. It was found that the exonic nsSNP (G/C) of WAK1 gene leading Cystine-Serine amino acid variation in the ORF region sequence (Table [133]2). Meanwhile, a WAK4 harboring large number of different intronic and exonic nsSNPs. Of which, the (G/A, A/G, G/A, T/C, and T/C) were located in the exon region, and exhibited the amino acid variation, such as Arginine-Glutamine, Lysine-Glutamate, Arginine-Histidine, Tyrosine-Arginine, and Cystine-Arginine. Furthermore, the exonic nsSNPs (A/T) of WAK5 was change the Tyrosine-Aspartate amino acid in the ORF region sequence, and the (C/G) exonic nsSNPs of OsDja9 was also responsible for Glycine-Arginine amino acid variation in the ORF region sequence. Afterward, the secondary and three-dimensional structures of selected genes were constructed. It was observed that the 3D-structure of all genes was changed due to variation of amino acid in the ORF region caused by exonic nsSNPs (Additional file 1; Fig. S5 A-D), which suggested that the SNPs in this region might be responsible for inactivation of these genes in D502. These results revealed that SNPs in WAK1, WAK4, WAK5, and OsDja9 could related to cope with M. oryzae infection in D506. Table 2. Analysis of four candidate genes harboring different single nucleotide polymorphisms between susceptible D502 and resistance D506 rice genomes. Locus_tag Description Variety SNP Position gene-OsI_08173 Hypothetical protein 117 390 D502 T *C (Ser) D506 C *G (Cys) Reference C *G (Cys) gene-OsI_24264 Hypothetical protein 179 445 452 513 520 601 735 897 982 1033 D502 *G (Arg) *A (Lys) *G (Arg) A *T (Try) *T (Cys) C T C C D506 *A (Gln) *G (Glu) *A (His) G *C (Arg) *C (Arg) T A A T Reference *G (Arg) *A (Arg) *G (Arg) A *T (Try *T (Cys) C A A T 1057 1136 1405 1585 1699 2043 2078 2105 2153 2232 D502 T C G G A A A T C T D506 C T A A G C T A T G Reference C C G G A C T T C T 2361 2379 2458 2462 2477 D502 G C G C G D506 A A A T A Reference G C G C G gene-OsI_21839 Hypothetical protein 301 912 2541 D502 *T (Tyr) C C D506 *A (Asp) A T Reference *A (Asp) A T gene-OsI_21393 Hypothetical protein 225 243 250 382 596 955 1535 1559 2092 2820 D502 T A *G (Gly) C G T G T G G D506 C T *C (Arg) G A C A A C T Reference C T *C (Arg) C G C A A G T 2834 2893 2945 3475 3952 3961 4028 4052 4078 D502 C G G T C G T C A D506 T T A C A A G A G Reference T T A C A A G A G [134]Open in a new tab Position are numbered based on NCBI reference sequence. The red astrike (*) indicates exonic SNPs lead amino acid variation in the susceptible variety, whereas the green astrike (*) represent exonic SNPs lead amino acid variation in the resistance variety. Discussion Rice blast caused by Magnaporthe oryzae is a highly destructive and damaging disease, resulting in severe losses to rice production. During favorable temperatures, the pathogen causes 70-80% yield losses in field [[135]2]. To mitigate the losses caused by this disease, there is a need to develop new disease resistance variety through breeding, which could harbor resistance gene with more effectiveness and have a sustainable approach to control not only rice blast pathogen but also other emerging pathogens [[136]23, [137]28]. Therefore, the present study employed two rice varieties, Diantun 502 (D502) and its bred Diantun 506 (D506), due to similar genetic background of both varieties, we further focus on understanding the mechanisms of broad-spectrum resistance to rice blast pathogen. Gene expression profile in rice D502 and D506 varieties Plants displayed diverse defense responses to invading pathogens in the form of signal molecules, defense genes, or multiple pathways in order to reduce their impact and survive in this unavoidable situation. Based on this, the present study was carried out to study the rice host defense transcripts between susceptible D502 and its new development resistance line D506 in the presence of M. oryzae during different time periods. Initially, a large number of DEGs 4911 (more than 67%) were found similar at both time hour points in D502 and D506, which were higher than from other researchers [[138]16, [139]29]. The GO enrichment analysis demonstrates that, after pathogen infection, mainly DEGs were significantly up-regulated in many biological processes and cellular component group at both time hours in D506, but not in D502 at 48h. The KEGG pathways analysis of DEGs showed that the genes were regulated in ribosome biogenesis in eukaryotes at both time hours in resistance variety, and in ribosome at only 24h in both varieties. Therefore, in response to M. oryzae, no DEGs were regulated in similar pathways at both time hours in D502, suggesting that the activation of DEGs were suppressed after pathogen infection. Moreover, the enrichment analysis of down-regulated DEGs were showed that genes were mainly involved in biological processes and molecular functions at 48h in both varieties. Therefore, an interesting result of KEGG pathways revealed that the DEGs of D502 were significantly down-regulated in many major pathways at both time hours, including plant hormone signal transduction, starch and sucrose metabolism, MAPK signaling pathways, and plant-pathogen interaction in response to M. oryzae. This result suggests that when the rice attacked by blast fungus, many DEGs were showed down-regulation upon M. oryzae infection, which might result in loss of resistance in D502. Previously, similar findings have also been reported that the down-regulated DEGs were regulated in major KEGG pathways of susceptible (parent line) variety and exhibited susceptibility to pathogen, as compared to (bred lines) resistance [[140]30]. This could be due to the initial stage of pathogen infection typically involves a complex exchange, perception, and signal transductions [[141]31]. For instance, a study reported that many KEGG pathways were down-regulated in disease resistance in susceptible maize parent lines as compared to their resistance lines after pathogen infection, revealing susceptibility in maize [[142]32]. Starch and sucrose metabolism are associated with starch synthesis, which is essential for growth and development in rice, while phenylpropanoids are secondary metabolites, produced from phenylalanine, serving different roles in plants such as structure and signaling functions, leading production of lignin, phytoalexin, and phenolic compounds [[143]33]. However, it is possible that the down-regulation of these pathways might be involved in susceptibility in D502. Furthermore, the gene expression profile in different clusters, obtained using an RNA-seq approach, suggested higher gene expression in cluster 3 after 24h and 48h of pathogen infection in resistance variety. The function analysis of DEGs revealed that they were mainly associated with defense pathways, including defense response to other organism, response to other organism, defense response, ADP binding, and metabolic pathways. Similarly, another study demonstrated that 231 genes in 1, 2, 4, and 12 clusters exhibited higher expression in resistance variety [[144]34]. Similar finding were reported in other recently studied cases, where only 4 clusters of DEGs exhibited high expression, which was regulated in many GO and KEGG pathways [[145]35]. Interaction network analysis on gene module detection and expression patterns between susceptible and resistant varieties Network analysis is a novel approach which employ to detect potential gene modules with the highest connectivity among genes [[146]36]. This analysis involves constructing a coexpression network that reveals correlations among genes across samples [[147]37]. It has been widely used in many studies to discover genes and explore the connection between find out genes and to explore the connection among gene expression patterns and related phenotypes [[148]36, [149]38, [150]39]. Therefore, we carried out a network analysis of DEGs of susceptible and resistant varieties at both time hours to classify different clusters based on their expression in both varieties. The nodes and edges were increased at 48h in all modules of D506, and decreased in D502 suggests that nodes and edges within modules were strongly correlated. Our data proved that the network interaction complexity was significantly increased in resistance variety, but decreased in the susceptible variety. The maximum average degree and lowest modularity were observed at 24h in both varieties, but the modularity significantly increased at 48h in D502 and D506. Furthermore, the DEGs were increased at 48h in all modules of D506 including (orange, purple, dark green and blue). Of this, the orange module contained maximum DEGs in D506, and out of which mostly were showed up-regulated expression. The DEGs that were found in purple module showed similar expression patterns at different hours. Many DEGs that were associated with dark green module exhibited down-regulated expression in D506. The DEGs that were related to the blue module displayed different expressions in D506 at both time hours. The functional enrichment analysis revealed that many DEGs were involved in similar GO and KEGG pathways, which were significantly decreased in susceptible variety at 48h, but increased in resistance D506 at the same time hours. While, many DEGs were regulated in defense and enzyme-related GO and KEGG pathways, including cell wall, chitin, abscisic acid biosynthetic process, ethylene-activated signaling pathway, MAPK signaling pathway – plant, plant hormones signals transductions, phenylalanine, tyrosine and tryptophan biosynthesis and phenylalanine metabolism at different time hours. The plant cell wall is the first barrier protecting rice from infection by the blast fungus through pathogenesis-related (PR) protein [[151]40]. However, there are only two PR proteins have been studied, including β-1,3-glucanases and chitinases, which belong to the GH-18 and GH-19 families [[152]41, [153]42]. Previously, the overexpression of chitinases 1 PR gene exhibit resistance in rice against rice blast fungus [[154]43]. Chenault et al. found that the chitinase and glucanase genes were highly expressed in resistance line, compared to their parent susceptible line, suggesting they showed resistance after pathogen infection [[155]44]. Similarly, the genes that were associated with ethylene pathway, MAPK signaling pathway – plant, plant hormones signal transductions, phenylalanine, tyrosine and tryptophan biosynthesis and phenylalanine metabolism showed similar expression patterns in both varieties. However, only five genes were expressed at 48h in D506, including phospho-2-dehydro-3-deoxyheptonate aldolase 2 (AROG2), macrophage migration inhibitory factor (MIFH), tryptophan synthase alpha chain (OsTSA), tryptophan synthase beta chain-like (TRPB2), tryptophan biosynthetic process (OASA2). The phospho-2-dehydro-3-deoxyheptonate aldolase enzyme is associated with phenylalanine and tryptophan biosynthesis, and both are essential for the production of chorismate biosynthesis [[156]45]. The chorismate biosynthesis is a vital intermediate in the biosynthesis of some important metabolites and is also associated with lignin biosynthesis [[157]45, [158]46]. Macrophage migration inhibitory factor (MIF) is a proinflammatory cytokine encoded by a genetically diverse locus. Recently, MIF has been identified as a versatile cytokine secreted by a variety of cell types, playing a crucial role in immune responses and various physiological processes [[159]47]. The tryptophan biosynthetic pathway supplies precursors for plant defense-related secondary metabolites and is recognized for its induction by pathogens and elicitors [[160]48]. It has been reported that the TSB1 gene is involved in plant defense against pathogens and is also a key regulator of plant growth [[161]49]. Our results are also supported to these statement that after infection of M. oryzae, mainly DEGs were regulated in defense and immunity-related pathways at different time hours, which might show resistance in D506. Transcriptome analysis revealed DEGs might involve in PTI and ETI signaling Plants have developed two immunity systems for combating the pathogen attack. The first system is PRRs, which identify the conserved PAMPs, and the second layer is ETI, which primarily acts within the cell and responds quickly to the host resistance by interacting with the pathogen effector and AVR gene, which shows defense reactions through direct or indirect [[162]1]. This ETI is related to a hypersensitive response (HR), which characterizes a system of programmed cell death and triggers an established of innate immunity signaling pathways that result in ion fluxes, creation of reactive oxygen species (ROS), and release of nitric oxide (NO) [[163]50]. However, the present study identified many DEGs with different expressions in D506, including P450 cytochrome, transcriptional factors (TF), NB-ARC-LRR, WAK, RLCKs, MAPK, RBOH, and CAM-CNG, which might be involved in the activation of both immunity system against M. oryzae attack in rice. Similarly, a study found many expressed DEGs only in near isogenic lines, encoded in NB-ARC, MAPK, RBOH, TF, and CAM [[164]51]. Whereas previous evidence suggests that the CYP71P1 member of P450 cytochrome family showed tryptamine 5-hydroxylase, catalyzing the transformation of tryptamine to serotonin, which increased the resistance against M. oryzae [[165]52]. Transcriptional factors (TFs) play a significant role in regulating the expression of many defense-related genes. For example, a study demonstrated that TFs were involved in the plant hormone signaling and regulated the expression of PR and disease resistance (R) genes [[166]53, [167]54]. In rice, many disease resistance genes have been reported and contain NB-ARC and LRR domains. These encode NLR protein and are associated with defense response and enhanced immunity against various phytopathogens [[168]27]. In the present study, we also found NB-ARC-NBS-LRR genes at 24h and 48h only in D506, suggested that these genes might be associated with enhanced immunity after M. oryzae infection. A recent study investigates the role of the NB-ARC protein RLS1 with cysteine-rich receptor-like protein RMC in rice, demonstrated that the activation of RLS1 triggers cell death, and increases disease resistance and immunity [[169]55]. Therefore, RLPs, RLKs, WAK, and SDR were also found only in D506. Cell surface PRRs, including RLK, RLP, and WAK, recognized P/M/DAMP to initiate the PTI signaling pathways in rice [[170]3]. The plasma membrane-localized WAK genes, which are strongly associated with plant cell walls, and have the ability to identify and release oligogalacturonides from the cell wall, recognize DAMPs, and activate immune response after M. oryzae infection in host [[171]34]. Additionally, RBOH and MAPK genes were also identified in this study, which showed up and down-regulated expression in response to M. oryzae. ROS and MAPK play an essential role in the first and second defense stages of plants, known as PTI and ETI. Upon recognizing PAMPs, PRRs regulate many immune responses, including MAPK and robust generation of ROS by OsRbohB and other components that enhance defense and immunity in rice [[172]56, [173]57]. These statements suggest that the identified DEGs in present study could be involved in PTI and ETI signaling, which enhance the disease resistance and immunity against M. oryzae in D506. Identification of non-synonymous SNPs between susceptible and resistant varieties and their effect on gene expression reveal resistance in D506 Integrating genome-wide polymorphisms data and transcriptomic reprogramming is becoming more common for classifying targeted genes associated with specific traits [[174]58]. Among the identified DEGs carrying nsSNPs between both varieties, only 41 genes were found with a different expression, which belongs to the defense and immunity-related families. Of these, only four genes, including WAK1, WAK4, WAK5, and OsDja9, which have been reported previously as disease resistance, exhibited up-regulated expression in D506. Wall associated receptor like kinase (WAKs) are known to have positive and negative regulators, that play a significant role in resistance and plant immune response by acting as DAMPs. They are also involved in many physiological processes, including defense response against M. oryzae and activated the cell elongation [[175]11]. Previously, it has been reported that WAK genes are essential contributors to enhancing immunity against many phytopathogens through various defense mechanisms [[176]8]. In Arabidopsis, the AtWAK1 is identified as cell-wall–derived oligogalacturonides (OGs) and regulates OG-mediated defense responses against different pathogens [[177]59]. Furthermore, it was documented that the overexpression of two WAKs, including WAK1 and WAK25 in transgenic rice plants enhances disease resistance to rice blast fungus [[178]60, [179]61]. In addition, DnaJ proteins have been reported in various organisms and play a significant role in the immunity of plants. A study reported that the overexpression of OsDja9 or knockout of OsDRP1E enhances resistance to the M. oryzae in rice transgenic plants [[180]62]. Therefore, we aligned the sequence of the ORF region of these genes and found that all genes harboring nsSNPs in this region and showed amino acid variations. Furthermore, the secondary and 3D protein structures were constructed, and it was found that the SNP altered 3D structures of selected genes were different. These results suggest that the SNP in the ORF region could be involved in up-regulated expression of all genes, that showed resistance against M. oryzae. Consequently, nsSNPs are single base changes responsible for amino acid variation in encoded protein sequences. They are mainly related to disease, which has been well-studied previously to check the impact of nsSNPs on individual proteins. For instance, a previous study investigated the effect of nsSNPs on protein-protein interaction and 3D structure variations, giving greater insight into the mechanisms by which nsSNPs can lead to disease [[181]63, [182]64], whereas the SNP are present on the active site of enzyme and binding site of receptors, might leading the change or loss of gene functions. It was documented that, the nsSNPs change amino acid in the protein sequence of CYP enzyme, and responsible for 3D structure variation and potentially disturb their metabolic activity [[183]65]. For instance, non-synonymous SNPs of the Pita allele was identified and compared them with reference and 47 Sri Lankan rice sequences, revealed resistance against M. oryzae [[184]66]. Another study demonstrated that Pib allele carrying nsSNPs, which leads to the variation of amino acids in the specific region, which might be involved in resistance [[185]67]. Molecular insights related to susceptible and resistant varieties In addition, we compare the important and significant molecular interactions involved in the form of four defense and immunity-related genes that harbor nsSNPs. These genes have a role in the activation of immune signaling and defense response upon pathogen infection. For instance, a gene transcript and genome analysis study between susceptible and resistant varieties demonstrated that SNP in the coding sequence of two up-regulated genes (CYSTM and NB-LRR) plays an important role in disease resistance through regulating salicylic acid (SA) pathways against Phytopathora infestans in tomato [[186]68]. The recent evidence demonstrated that the genes harboring nsSNPs displayed up-regulated expression in plant hormone signal transduction and phenylpropanoid pathways at an early stage of Fusarium ear rot pathogen infection in maize resistance line as compared to susceptible (parent) line [[187]69]. While, integrated analysis in maize revealed two resistance genes encoding in calmodulin-like protein and LRR-RLK harboring nsSNPs in resistance cultivar, which enhance the amount of carotenoids and SA at the early stage of pathogen infection and showed defense response against grey leaf spot pathogen in resistance variety as compared to susceptible [[188]70]. Furthermore, the previous study applied GWAS and identified Pb4 NLR gene, which belongs to the WAK family, and serves as a key gene against M. oryzae by regulating PTI immune responses [[189]71]. It has been documented that WAK genes are involved in defense response through positive or negative regulation, which may be triggered by chitin derived from fungal cell walls [[190]11]. In recent study, the zmCPK39 resistance module was found in maize against multiple foliar diseases [[191]72]. Similarly, ZmWAKLY module enhances resistance in maize through regulating immune signaling and defense response after pathogen infection in resistance variety [[192]73]. However, this method of association mapping based on SNPs is effective for annotating gene functions and identifying genetic regulatory networks, making it valuable for investigating the roles and interactions of genes within the plant-pathogen interaction pathway. For instance, the SNP in the promoter and coding region sequence of up-regulated genes, including Pi68, RGA2, RGA4, NBS-LRR, and LRR-RLK located on chromosomes 1, 2, 3, and 6, play an important role in disease response through regulating various defense and immunity-related mechanisms and also enhance resistance against pathogen infection in resistance line as compared to susceptible parent line [[193]74]. This evidence from different crops suggests that the genes carrying different SNP between susceptible and resistance varieties have a possible role in activating defense mechanisms in rice Diantun 506 against M. oryzae. Conclusions Based on overall study, we concluded that both the susceptible (D502) and resistance (D506) varieties in response to Magnaporthe oryzae infection responds differently when confirmed through the integration of multiple omics-based approaches. Firstly, our study investigated four defense and disease resistance genes, which showed high up-regulated expression patterns after pathogen infection in resistance variety and carrying different nsSNPs. Secondly, on the basis of SNP and high expression pattern of genes, it is possible that the identified genes might be involved in resistance and enhance the immunity of resistant rice variety D506 against rice blast pathogen, M. oryzae. Finally, these findings could offer important genetic resources for molecular breeding in agricultural practices for disease management against rice pathogens. Supplementary Information [194]Supplementary Material 1. ^(2.6MB, docx) [195]Supplementary Material 2. ^(2.6MB, xlsx) Acknowledgments