Weighted Gene Co-Expression Network Coupled with a Critical-Time-Point Analysis during Pathogenesis for Predicting the Molecular Mechanism Underlying Blast Resistance in Rice

Rice blast, caused by the ascomycete fungus M. oryzae, is one of the most important diseases of rice. Although many blast resistance (R) genes have been identified and deployed in rice varieties, the molecular mechanisms responsible for the R gene-mediated defense responses are yet not fully understood. In this study, we used comparative transcriptomic analysis to explore the molecular mechanism involved in Piz-t-mediated resistance in a transgenic line containing Piz-t (NPB-Piz-t) compared to Nipponbare (NPB). Clustering and principal component analysis (PCA) revealed that the time-point at 24-h post inoculation (hpi) was the most important factor distinguishing the four time-points, which consisted of four genes of mitogen-activated protein kinases (MAPKs) signaling pathway, one gene related to WRKY DNA-binding domain containing protein, five pathogenesis-related protein (OsPR1s) genes, and three genes of R proteins involving in the most significant protein-protein interaction (PPI) pathway. Using weighted gene co-expression network analysis (WGCNA) to investigate RNA-seq data across 0, 24, 48, and 72 hpi, nine modules with similar patterns expression pattern (SEP) and three modules with differential expression pattern (DEP) between NPB-Piz-t and NPB across 0, 24, 48, and 72 hpi with KJ201 (referred to as Piz-t-KJ201 and NPB-KJ201) were identified. Among these the most representative SEP green-yellow module is associated with photosynthesis, and DEP pink module comprised of two specific expressed nucleotide-binding domain and leucine-rich repeat (NLR) genes of LOC_Os06g17900 and LOC_Os06g17920 of Pi2/9 homologous, three NLR genes of LOC_Os11g11810, LOC_Os11g11770, and LOC_Os11g11920 which are putatively associated with important agronomic traits, and a B3 DNA binding domain containing protein related genes (LOC_Os10g39190). Knockout of LOC_Os10g39190 via CRISPR-Cas9 resulted in plant death at the seedling stage. The research suggested that Piz-t and multiple NLR network might play important roles in the regulation of the resistance response in the Piz-t-KJ201 interaction system. The identified genes provide an NLR repository to study the rice-M. oryzae interaction system and facilitate the breeding of blast-resistant cultivars in the future.


Background
Rice blast is a major disease threatening global rice production, causing a 10-30% global crop yield losses. Utilization of rice-blast resistant cultivar is the most economical and environment-friendly way of minimizing crop losses. Currently, over 100 rice-blast R genes have been mapped and at least 32 have been characterized and cloned (Xiao et al. 2020;Zhao et al. 2018). However, only few broad-spectrum and durable resistant cultivars have been bred due to the highly complex and dynamic process associated with blast resistance response (Ashkani et al. 2016).
The completion of the whole genome sequences for rice and M. oryzae has made the pathosystem become a premier model for studying plant-fungal interactions, which increases our understanding of the mechanisms underlying fungal infections. The major members involved in rice-M. oryzae interactions include pattern-recognition receptors (PRRs) and NLR genes from rice, which trigger a diverse array of immune responses including energy metabolism, pathogen recognition, defense related proteins, hormone signaling, ROS, and redox homeostasis, especially on the changes of genes expression and transcriptional reprogramming. Of these, R genes are well studied as the famous gene for gene resistance; meanwhile, a large NLR proteins immune signaling network responses to invading pathogens and confers more durable resistance than single race-specific R gene . Understanding of these genes mediated defense response to blast can be useful in helping resistance breeding.
Previously, proteome and transcriptome analyses have revealed many genes involving in the defense responses of rice to M. oryzae (Wei et al. 2013;Wang et al. 2014;Zhang et al. 2016;Jain et al. 2017;Tian et al. 2018), however, little is known about the potential connection among these genes. Recently, several approaches have been developed to look at genes with similar expression patterns which may participate in specific biological functions. One approach is the method named WGCNA, which was implemented as a R package to identify significant associated genes modules with similar expressions. This reduces the dimension of complex data and simplifies it into several modules, thereby providing an overview of gene-gene inter-relationships at the system level. In addition, clustering is a useful exploratory technique for analysis of gene expression data, and PCA can help to reduce the dimensionality of the data set by transforming to a new set of variables to summarize the features of the data. Thus, integrating these two analysis techniques to analyze the four time-point data can contribute to determine the critical one (Yeung and Ruzzo 2001).
The rice blast R gene Piz-t is a member of the Pi2/9 multi-allelic gene family on chromosome 6 , which follows a gene-for-gene fashion to the M. oryzae AVR gene AvrPiz-t (Li et al. 2009a, b). This gene was once widely used in breeding programs for increasing resistance to M. oryzae Tian et al. 2020). Numerous studies also have examined the interaction profiles of Piz-t-AvrPiz-t in rice (Park et al. 2012(Park et al. , 2016Wang et al. 2016;Tang et al. 2017;Tian et al. 2018;Bai et al. 2019;Zhang et al. 2020). However, transcriptome profiling underlying the Piz-t-mediated rice resistance to blast fungus remains largely unknown.
Generally, the whole infection process begins at the first time infected at 24 hpi, and then invades the neighboring cells at 48 hpi (Cao et al. 2016). Subsequently, rice blast disease lesions will emerge about 72 hpi (Wilson and Talbot 2009). The objective of this study then is to identify a set of candidate defense genes mediated by Piz-t with M. oryzae using the clustering, PCA, and WGCNA approach across the four time-points. Notably, a larger number of R genes and regulators were identified during pathogenesis. The identified genes can be further studied to reveal the molecular mechanisms underlying rice-blast interactions in the future.

RNA-Seq and Validation of RNA-Seq Data by qRT-PCR
Leaf tissues were collected from each rice line at 0, 24, 48, and 72 hpi for RNA-seq. Some plants were retained for investigation of disease progression. There were some spindle-shaped lesions on NPB leaves and no obvious disease phenotype on Piz-t-KJ201 leaves after 7 days of incubation according to a disease rating of 5 (on a 0-5 scale) (Mackill and Bonman 1992) (Fig. 1a). RNA-Seq analysis of Piz-t-KJ201 and NPB-KJ201 were conducted as shown in Fig. 1b. A total of 1,242,379,052 high-quality reads were generated, and over 90% of the reads in each sample successfully aligned to RGAP7, with unique mapped reads at > 85%, and multiple mapped reads or fragments were < 6%. To confirm the gene expression patterns identified from the RNA-seq data, the transcript levels of 18 important genes were examined and validated by qRT-PCR (Additional Table 1). All the gene expression patterns obtained were found in consistency with the RNA-seq data (Fig. 1c).
Results from the venn diagram showed that the upregulated differentially expressed genes (DEGs) were 406, 175, 98, and 104, and the down-regulated ones with 50, 71, 33, and 280 at 0, 24, 48, and 72 hpi, respectively, and 49 common DEGs at the four time-points (Fig. 2a,  b). Clustering analysis also revealed that those DEGs at the four time-points have consistent expression levels among their three replicates, respectively. Although part of those genes expression levels at Piz-t-0hpi.3 were different with that at Piz-t-0hpi.1 and Piz-t-0hpi.2, they had higher level of consistency with Piz-t-0 hpi compared with NPB-0 hpi (Additional Fig. 1.A, B, C, D). Further comparisons of the DEGs across the four time-points showed that most of the DEGs were up-regulated at 0 hpi and down-regulated at 24 hpi. In contrast, almost opposite expression patterns were seen between 0 and 24 hpi. Also, the DEGs at 24 hpi were significantly different from the other three time-points as shown by Fig. 2c. PCA on 24 group data related to Piz-t-KJ201 and NPB-KJ201 with three replicates revealed that the first principal component (PC), PC1, provided the most information on time-points that determine expression difference (Fig. 3), suggesting that the time-point at 24 hpi is the most significant factor that triggers global defense response in gene expression.
Previous studies suggested that the pathogens finish the first time infected at 24 hpi, accordingly, the rice activates defense responses against M. oryzae before 24 h (Bagnaresi et al. 2012;Jain et al. 2017). Thus, this timepoint at 24 hpi is the most intriguing for further extensive exploration. The workflow for leaf transcriptome of Piz-t-KJ201 and NPB-KJ201 using RNA-Seq. c The qRT-PCR validation of the differentially expressed gene between NPB-Piz-t and NPB inoculated with KJ201 in the PPI pathway at 24 hpi. The data were log2 transformed for FC Differentially Transcriptional Profiles between Piz-t-KJ201 and NPB-KJ201 at 24 hpi At 24 hpi, 175 and 71 genes were identified to be upregulated and down-regulated in Piz-t-KJ201 compared to NPB-Piz-t, respectively (Supplemental Table 2). To provide valuable insight into the functional associations of the 246 DEGs at the critical time-point, GO analysis of these DEGs showed that the most significant categories of biological processes were associated with auxin catabolic processes, oxazole or thiazole biosynthetic and metabolic processes. The most significant cellular component categories correlated with extracellular region, intrinsic and integral component of plasma membrane, and plasma membrane, and the most molecular function categories related to ADP binding, monooxygenase activity, oxidoreductase activity, and ion binding (Additional Fig. 2A). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that the DEGs functions were mainly involved in PPI, lysine degradation, MAPK signaling pathway, cAMP  signaling pathway (Additional Fig. 2B). Several transcription factors (TF) such as FAR1, WRKY, NAC, and bHLH were also significantly enriched in the critical time-point (Table 1). Next, we selected the genes of the significant PPI pathway for further analysis. Notably, four genes of MAPK signaling pathway, one WRKY-related gene, five PR1 genes, and genes of eight R proteins were involved in this pathway. All of these genes had significantly expressed at higher levels in Piz-t-KJ201 compared to NPB-KJ201 except LOC_Os02g41710, which had higher expression level in NPB-KJ201 (Table 1, Fig. 1c). Interestingly enough, the expression ratios of LOC_Os06g17920 and LOC_ Os06g17900 between the Piz-t-KJ201 and NPB-KJ201 were very large, and both of these two genes had small amount of expression levels in the NPB-KJ201 across the four time-points (Table 2, Fig. 1c). In a similar vein, the other four genes also had higher or specific expression in Piz-t-KJ201 (Table 1, Fig. 1c).
As the MAPKs are conserved signaling molecules that transduce pathogen stimuli into intra-cellular responses, those up-regulated MAPKs in Piz-t-KJ201 contributed to the activation of PR1s and R genes, leading to the synthesis of antimicrobial compounds such as secondary metabolites that prevent pathogen development. Additionally, several active TFs from FAR1, WRKY, NAC, and bHLH families have been demonstrated to be involved in responses to biotic and abiotic stresses (Ramamoorthy et al. 2008), which were implicated in the regulation of transcriptional reprogramming associated with early response to M. oryzae. Taken these genes together showed that expanded transcriptional activation may be controlled by Piz-t in the PPI pathways during the early response of rice to M. oryzae.

Gene Co-Expression Networks
The WGCNA R package was used to construct a coexpression network consisting of the normalized read counts for the 24 samples based on 33,451 genes. The results showed 12 co-expression modules in gene membership, ranged from as low as three (grey) to as much as 5734 (turquoise) genes, as shown in Additional Tables 3,  4 , 5, 6, 7, 8, 9, 10, 11, 12, 13 and 14. These involved genes reflect that the rice-M. oryzae interaction is quantitative, and they represent the complexity of a cellular transcription network, from which nine SEP and three DEP modules were enriched, respectively, between NPB-KJ201 and Piz-t-KJ201 (Fig. 4). The nine SEP modules consisted of genes with similar expression profiles between NPB-KJ201 and Piz-t-KJ201 at the four time-points, suggesting the basal similarities between their transcriptional profiles; and the three DEP modules contained those genes with unique expression profiles in NPB-KJ201 or Piz-t-KJ201 at the four time-points, meaning that they may directly associate with resistance response.
The nine SEP was involved in four different types of expression pattern: (i) the turquoise and green modules: up-regulated at 0 and 24 hpi, and then down-regulated at 48 hpi, but up-regulated at 72 hpi again (Additional Tables 3 and 4); (ii) the blue, brown, black, and magenta modules: down-regulated at 0 and 24 hpi, and then upregulated at 48 hpi, but down-regulated at 72 hpi again (Additional Tables 5, 6, 7, 8); (iii) the red and yellow modules: up-regulated at 0, 24, and 48 hpi, but downregulated at 72 hpi (Additional Tables 9 and 10); (iv) the green-yellow: down-regulated at 0 and 24 hpi, and upregulated at 48 and 72 hpi, with an obvious higher expression level in Piz-t-KJ201 across the four time-points compared to NPB-KJ201 (Additional Table 11). On the other hand, the DEP modules were composed of the pink, purple and grey ones. All of these three DEP modules obviously had different expression patterns between Piz-t-KJ201 and NPB-KJ201 across the four time-points (Fig. 4,Supplementary Tables 12,13,14). Although the purple model included some obvious DEGs between Piz-t-KJ201 and NPB-KJ201, they didn't occur at the key time-point of 24 hpi (Supplementary Tables 13). Additionally, the grey modules contained as little as 3 complete reverse regulated DEGs between Piz-t-KJ201 and NPB-KJ201, and they also had lower levels expression in Piz-t-KJ201 compared to NPB-KJ201 at 24 hpi (Supplementary Tables 14). Intriguingly, the pink module contains a set of specific higher DEGs in Piz-t-KJ201 compared to NPB-KJ201, especially at 24 hpi, which suggests that this module involved in DEGs was unique to the resistant response (Supplementary Tables 12).

Differentially Transcriptional Profiles in the Pink and the Green-Yellow Module
We further analyzed the transcriptional profiles of the representative pink and green-yellow modules, and the results showed that 106 and 34 genes that were differentially regulated between Piz-t-KJ201 and NPB-KJ201 (Additional Tables 11 and 12). Gene ontology (GO) analysis of DEGs in the pink module showed that the most significant molecular function categories identified were ADP binding, adenyl ribonucleotide binding, adenyl nucleotide binding, ion binding, and nucleotide binding. Meanwhile, the most significant biological processes were associated with DNA integration, cell communication, and signaling transduction. In the case of DEGs in the green-yellow module, the most significant molecular function categories were chlorophyll binding, electron transporter and transfer activity, and tetrapyrrole binding. The most significant biological process categories were photosynthetic electron transport in photosystem II, protein-chromophore linkage, and photosynthesis process, while cellular component categories identified include those associated with photosystem, thylakoid, and membrane protein complex. KEGG analysis also revealed that DEGs specific to the pink module were mainly involved in pathways related to PPI, while DEGs specific to the green-yellow module were mainly associated with the photosynthesis (Table 3).
Five resistance genes of LOC_Os11g11810, LOC_ Os06g17920, LOC_Os11g11770, LOC_Os06g17900, and LOC_Os11g11920 involved in PPI pathway and one TF of LOC_Os10g39190 are the most significant element in the pink module (Table 2, Fig. 1c, Additional Fig. 3). Interestingly, LOC_Os06g17900 and LOC_Os06g17920, two putative pseudogenes of the Pi2/9 homologous showed the first and second highest values, respectively, and three activated NLR genes of LOC_Os11g11810, LOC_Os11g11770, and LOC_Os11g11920 in Piz-t-KJ201 have been reported to associate with important agronomic traits. Another TF of LOC_Os10g39190 perhaps involved in ABA signaling regulation was significantly up-regulated in Piz-t-KJ201 plants across four time-points, especially at 24 hpi (Table 2, Fig. 1c, Additional Fig. 3).
Next, we focused on the functional analysis of LOC_ Os10g39190. One gRNA locating in the 345 bp of Os10g053710 was designed to silence this gene (Additional Fig. 4). Although we obtained six CRISPR mutants, but all of these transgenic plants didn't survive at the seedling stage.

Discussion
In this study, we showed an approach in evaluating large transcription data using co-expression networks in combination with the critical time-point at 24 hpi analysis. Specifically, we applied this strategy in analyzing gene expression dataset from rice in response to M. oryzae infection and identifying genes possibly involved in rice against M. oryzae, In this study, we identified several significant DEGs related with TF, MAPK signaling pathway, OsPR1, and R proteins of genes to be possibly involved in the resistance response in the Piz-t-KJ201 interaction, facilitating understanding the interaction of rice-M. oryzae and breeding for new resistant cultivars.

Resistance Response Involved in the Regulation of Photosynthesis
Previously, it has been reported that both pathogenassociated molecular pattern (PAMP) -triggered immunity (PTI) and effector-triggered immunity (ETI) contribute to the resistance response to rice blast (Boller and Felix 2009), and NLR-mediated immunity is known to be modulated by many environmental factors such as light (Gao et al. 2020). In the present study, nine SEP between Piz-t-KJ201 and NPB-KJ201 were found, especially green-yellow module, in which the photosystem plays a vital role. Previous studies  confirmed that photosystem was easily targeted by different pathogenic species such as bacteria, viruses, fungi and oomycetes to produce cytoplasmic calcium bursts, jasmonic acid, salicylic acid, and reactive oxygen species, resulting in transcriptional reprogramming to mount a full defense response (Fondong et al. 2007;Jelenska et al. 2007;Rodriguez-Herva et al. 2012;Nomura et al. 2012;Li et al. 2014a, b;Petre et al. 2016;Sowden et al. 2018;Liu et al. 2018;Rosas-Diaz et al. 2018). Subsequently, the photosystem needs to be restored in time to prevent or reduce the damage and energy costs, which can be obviously shown in the green-yellow module with higher expression in Piz-t-KJ201.  (Monné et al. 2015). Although the other three MAPK signaling pathway genes were up-regulated in Piz-t-KJ201, complete understanding of their roles in rice disease resistance still remains unclear. Interestingly, five OsPR1 encoding SCPlike extracellular protein were identified to be upregulated in Piz-t-KJ201, suggesting that they may play particular roles in resistance response. Similarly, there were also evidences that the expression of PR1a and PR1b were related to disease resistance (Ponciano et al. 2006). Several TFs of WRKY, bHLH, AP2/ERF, Dof, MYB, and NAC etc. were significantly differentially expressed between Piz-t-KJ201 and NPB-KJ201. These members of the TF families have been confirmed to be involved in regulating rice defense responses to M. oryzae (Licausi et al. 2013;Yokotani et al. 2014;Cheng et al. 2015). Thus, they may play important roles in the regulation of defenseresponsive genes. The WRKY104 was a positive regulator enhancing the resistance to M. oryzae by up-regulating pathogenesis-related (PR) genes as also shown in recent reports (Hou et al. 2019). Additionally, the higher expression of R genes LOC_Os11g11990, LOC_Os11g43700, and LOC_Os12g17430 in the Piz-t-KJ201 interaction system were also observed, which have also been implicated in stress responses (Zheng et al. 2014).

Several Putative NLR Pairs with Piz-t Helper to Mediate Immune Signaling
Plant R genes detect effector proteins secreted by pathogens by indirectly or directly binding them via effector-targeted host proteins, with a complex scenario of sequential action of "sensor" NLRs to detect pathogen effectors and "helper" NLRs to initiate immune signaling . The pink module is the most representative DEP, which includes 6 significant DEGs between Piz-t-KJ201 and NPB-KJ201 across the four time-points. The alleles of LOC_ Os06g17920 and LOC_Os06g17900 of Piz-t locus were regarded as R pseudogenes in the rice Nipponbare genome (Wu et al. 2012). However, when the prediction was based on our present study, these two pseudogenes were actually not pseudogenes but always highly expressed in Piz-t-KJ201 compared to NPB-KJ201 (Table 2). Thus, they are maybe functionally redundant with Piz-t as helper NLR activated by sensors NLR. The specific expressions and characteristics of the alleles of LOC_Os06g17920 and LOC_ Os06g17900 in Piz-t-KJ201 deserve further explorations. In addition, LOC_Os11g11810 and LOC_Os11g11770 were identified to be known Pia and RPM1, respectively (Okuyama et al. 2011;Yoo et al. 2017;Onaga et al. 2017). Previous studies identified that Pia encodes a NBS-LRR disease resistance protein that are known to contribute to weak resistance to rice blast (Okuyama et al. 2011;Césari et al. 2013;Jung et al. 2013); while RPM1 encodes NB-ARC domain containing protein and can be activated by the defense protein RIN4 to sense effector, and then transmit signals downstream (Li et al. 2014a, b), suggesting that the two genes confer to resistance as the sensor NLRs. Moreover, both of LOC_ Os11g11810 and LOC_Os11g11770 were identified to be nonsynonymous single nucleotide polymorphism NLR genes, and another NLR gene of LOC_Os11g11920 was also down-regulated at developmental period in black rice (functional cultivar) (Seol et al. 2016), implying that they likely play an important role in determining agronomic traits differences. This is consistent with the previous report that a complex relationship exist between disease resistance and yield-related components .
The B3-domain transcription factor abscisic acidinsensitive 3 (ABI3) was a central regulator in ABA signaling (Zhang et al. 2005). However, to date, there are still a few examples of functional validation of the roles being played by B3 in plant defense. Unfortunately, in the functional identification of LOC_Os10g39190 via CRISPR-Cas9, we can't get the transgenic seed as all of the T 0 generation of transgenic plants died at the seedling stage. In general, the effect of knocking out 1 member of the paired NLRs was more severe in some sensor mutants compared to some helper mutants, which usually led to the death of the mutants , from which we can infer that LOC_Os10g39190 maybe as a potential important sensor member.

Validation of the Present Study by Comparing Previous Reports
In our earlier studies, four proteins of a receptor-like protein kinase (gi|59,800,021), a NADP dependent malic enzyme (gi|54,606,800), a putative transaldolase (gi|57, 900,129), and a putative bowman birk trypsin inhibitor (gi|53,792,234), three C3HC4-type Ring finger E3 ubiquitin ligase (APIP2, 6, 10), and a bowman-Birk type bran trypsin inhibitor (APIP4) were confirmed to be involved in Piz-t-KJ201 interaction (Park et al. 2012(Park et al. , 2016Tian et al. 2018;Zhang et al. 2020). In the present study, we observed that a receptor-like protein kinase precursor (LOC_Os10g33040), a NADP-dependent malic enzyme (LOC_Os05g09440), a transaldolase (LOC_ Os08g05830), a bowman-Birk type bran trypsin inhibitor precursor (LOC_Os01g03360), and two zinc finger C3HC4 type domain containing protein (LOC_ Os01g20910 and LOC_Os12g04590) were also differentially regulated in Piz-t-KJ201 (Additional Table 2). Therefore, parts of the present results were validated by the previous studies. Though we tried to identify the function of LOC_Os10g39190, it's regretful that the mutant can't survive at the seedling stage. This research reveals that many genes are involved in the defense network of Piz-t-KJ201, consistent with previous observations (Park et al. 2012(Park et al. , 2016Wang et al. 2019;Bai et al. 2019;Zhang et al. 2020). These genes and TFs identified in the present study were of high biological relevance. It was their pairs and interactive NLR network that confer blast resistance in the Piz-t-KJ201 interaction as described in the former researches Baggs et al. 2017;Wang et al. 2019).

Conclusion
In this study, we identified that the 24 hpi is the critical time-point during early stage pathogenesis, from which four genes of MAPKs signaling pathway, one gene related to WRKY DNA-binding domain containing protein, five OsPR1s genes, and three genes of R proteins were confirmed involved in the significant PPI pathway. In addition, nine SEP modules and three DEP modules between NPB-Piz-t and NPB were also identified, among which the most representative SEP green-yellow module is associated with photosynthesis, and DEP pink module comprised of two specific expressed NLR genes of LOC_Os06g17900 and LOC_Os06g17920 of the Pi2/9 homologs, three NLR genes of LOC_Os11g11810, LOC_Os11g11770, and LOC_ Os11g11920 which are putatively associated with important agronomic traits, and a B3 DNA binding domain containing protein related genes (LOC_Os10g39190). Knockout of LOC_Os10g39190 via CRISPR-Cas9 was done, however, transgenic plants failed to survive at the seedling stage. It suggested that Piz-t with these NLRs and regulators network might play important roles in the regulation of the resistance response in the Piz-t-KJ201 interaction system. The identified genes provide a valuable resource to study the rice-M. oryzae interaction system and facilitate breeding of new blast-resistant cultivars in the future.

Plant Materials and Blast Isolates
Rice (Oryza sativa L.) lines, including NPB and its transgenic line with Piz-t (NPB-Piz-t), and the M. oryzae isolates KJ201 was used in this study (Tian et al. 2018). Rice plants were grown in the greenhouse and kept under natural conditions about 2~3 weeks (3-4 leaves). Spores concentration in the suspension was adjusted to 5 × 10 5 conidia/mL. After spray-inoculated, the seedlings were maintained in the dark for 24 h at 28°C and then kept under high humidity over 95% for about a week for symptom evaluation.
Leaf tissues were collected from each rice line at 0, 24, 48, and 72 hpi and frozen in liquid nitrogen. Some leaves of NPB-KJ201 showed obvious disease lesion, while Piz-t-KJ201 had no disease symptoms at 7 days post inoculation (Fig. 1a).

Sample Preparation
The total RNA from the infected leaves was extracted using the RNAeasy kit (Qiagen, Germany) and treated with an RNase-Free DNase Set (Qiagen, Germany), according to the manufacturer's instructions. The RNA quality, library construction and size were assayed using a 2100 Bioanalyzer system (Agilent, USA). The libraries were synthesized using the TruSeq RNA Sample Preparation v2 kit (Illumina, USA). Total RNA from each treatment was pooled and then libraries were constructed as showed in Fig. 1b and used for sequencing. The samples were run in the NovaSeq system and raw sequences of paired 150-bp were obtained.

Sequences Processing and Analysis
FastQC was used to estimate the raw reads and assess their quality. Trimming of reads was carried in Trimmomatic and reads containing contaminant primer/adapters and long stretches of poor-quality bases were removed. Each sample was mapped back onto O. sativa reference genome from the MSU Rice Genome Annotation Project (MSU Rice Genome Annotation Project Release 7, RGAP7) (Kawahara et al. 2013). Redundancies were removed using Picard suite of tools (http://broadinstitute. github.io/picard/). Alignment was performed using the metrics module of SAMtools (Li et al. 2009a, b). And read counts were estimated using featuresCounts (Anders et al. 2015). Normalization and differential expression analysis was performed using EdgeR (Anders and Huber 2010).
For a given genotype, significant DEGs at each timepoint were identified based on fold change (FC) (ratio of number of transcripts in inoculated sample to the number of transcripts in the mock control ≥2 or ≤ 0.5 and a p-value < 0.05).

Bioinformatics Analysis
Clustering and PCA were done using the R package. Biological processes, cellular components, and molecular functions of DEGs were determined by GO database annotation (http://www.geneontology.org/). Protein signaling pathways were elucidated using the KEGG database (http://www.genome.jp/kegg/ pathway. html). Pathways enriched at P-value < 0.05 were considered significant.

Co-Expression Network Analysis
Using RSEP version 1.1.1 to quantify expression levels of genes in terms of FPKM (fragments per kilo base of exon per million)1 (Li and Dewey 2011), a common expression matrix for Piz-t-KJ201 and NPB-KJ201 at 0, 24, 48, and 72 hpi was constructed, with three replicates from each condition resulting in a total of 24 samples. Co-expression analysis was performed using R package. The similarity between gene-pairs was computed using a signed Pearson's correlation matrix, scaled to power β = 11, based on the approximate scale free-topology criterion.

Real-Time PCR Validation
Real-time PCR (qRT-PCR) was used to validate DEGs obtained from RNA-Seq. Primers were designed using Primer5 and summarized in Supplementary File 2. The ProtoScript M-MuLV First Strand cDNA Synthesis Kit (NEB) was used to synthesize cDNA. The reaction was performed in the LightCyclerR 480 II PCR system (Roche) with a volume of 10 μl, containing 5 μl of SYBR Green I Master mix, 0.8 μM of forward primer, 0.8 μM reverse primer, 2 μl of 1:5 diluted cDNA template (1-