Understanding Host-Pathogen Interactions with Expression Profiling of NILs Carrying Rice-Blast Resistance Pi9 Gene

Magnaporthe oryzae infection causes rice blast, a destructive disease that is responsible for considerable decrease in rice yield. Development of resistant varieties via introgressing resistance genes with marker-assisted breeding can eliminate pesticide use and minimize crop losses. Here, resistant near-isogenic line (NIL) of Pusa Basmati-1(PB1) carrying broad spectrum rice blast resistance gene Pi9 was used to investigate Pi9-mediated resistance response. Infected and uninfected resistant NIL and susceptible control line were subjected to RNA-Seq. With the exception of one gene (Pi9), transcriptional signatures between the two lines were alike, reflecting basal similarities in their profiles. Resistant and susceptible lines possessed 1043 (727 up-regulated and 316 down-regulated) and 568 (341 up-regulated and 227 down-regulated) unique and significant differentially expressed loci (SDEL), respectively. Pathway analysis revealed higher transcriptional activation of kinases, WRKY, MYB, and ERF transcription factors, JA-ET hormones, chitinases, glycosyl hydrolases, lipid biosynthesis, pathogenesis and secondary metabolism related genes in resistant NIL than susceptible line. Singular enrichment analysis demonstrated that blast resistant NIL is significantly enriched with genes for primary and secondary metabolism, response to biotic stimulus and transcriptional regulation. The co-expression network showed proteins of genes in response to biotic stimulus interacted in a manner unique to resistant NIL upon M. oryzae infection. These data suggest that Pi9 modulates genome-wide transcriptional regulation in resistant NIL but not in susceptible PB1. We successfully used transcriptome profiling to understand the molecular basis of Pi9-mediated resistance mechanisms, identified potential candidate genes involved in early pathogen response and revealed the sophisticated transcriptional reprogramming during rice-M. oryzae interactions.


INTRODUCTION
Rice is among the most important staple foods, contributing 23% of total calories consumed globally. Over 600 million tons of rice is produced annually from 150 million hectares of rice paddies. Asia is responsible for 92% of that production, with three-quarters of the output from India and China (Maclean et al., 2002). Rice is the model system for monocotyledons because its genome has been fully sequenced (International Rice Genome Sequencing Project, 2005;430 Mb), the availability of high-density genetic maps, genome-wide microarrays (Jung et al., 2007), and genetic transformation methods. Both biotic and abiotic stresses affect rice growth. Among biotic stresses, infection by Magnaporthe oryzae (a hemibiotrophic fungus) causes rice blast disease, resulting in a 20-100% global crop loss, an amount that can feed 60 million people .
The development of new molecular tools and availability of model plant-fungal systems have increased our understanding of the mechanisms underlying fungal infections. Rice blast is one of the model diseases under study because annotated genome sequences are available for both organisms (M. oryzae: Dean et al., 2005;rice: Ohyanagi et al., 2006). The complete genome sequences of the cultivated rice strains Oryza sativa L. ssp. japonica and ssp. indica also allow for genome-wide transcriptome analyses (Jiao et al., 2009;Sato et al., 2011).
The gene-for-gene hypothesis operates in the rice-M. oryzae system (Flor, 1971). Infection on the rice plant starts with fungal spore attachment to the leaf, followed by germination and appressoria formation within 24 h to penetrate the leaf cuticle, invading epidermal cells (Talbot, 2003;Skamnioti and Gurr, 2009). Effector-triggered immunity (ETI) occurs when the M. oryzae Avr (avirulence) protein recognizes rice R (resistance) proteins, leading to a hypersensitive response that stops fungal growth. Pathogen triggered immunity (PTI) is triggered if fungal pathogen-associated molecular patterns (PAMP; e.g., chitins) is recognized by rice pattern-recognition receptors (PRR) and fungal hyphae spread inside plant causing disease symptoms (Liu et al., 2010;Chen and Ronald, 2011). Of the two defense responses, ETI is stronger and faster (Tao et al., 2003). Riceblast resistant varieties can eliminate pesticide use and minimize crop losses, but high pathogen variability and little mechanistic understanding of R-mediated resistance pathways mean that newly developed resistant cultivars are often susceptible after only a few years. Currently, over 100 rice-blast R genes (50% indica, 45% japonica, 4% wild species) have been mapped, but only 25 have been characterized and cloned (Sharma et al., 2016).
Studying transcriptome dynamics provides insight into functionally important genomic elements, their expression patterns, and their regulation in different developmental stages, tissues, and environmental stressors (Wang et al., 2011). Several methods are available for transcriptome isolation and quantification, including expressed sequence tag (EST) library sequencing, serial analysis of gene expression (SAGE), and SuperSAGE. However, these methods have numerous disadvantages, including low throughput, cloning bias, low sensitivity, and high cost. Microarray hybridization avoids some of these issues and thus sees common use, but the technique is prone to high background noise and requires known gene sequences, meaning it cannot identify novel transcribed regions. Fortunately, next-generation RNA sequencing technology (mRNA-Seq) is now available as a high-throughput method for simultaneously sequencing, mapping, and quantifying transcriptome reads. It is a robust tool for identifying rare or novel transcripts, alternative splice junctions and transcription start sites .
Monogenic or near-isogenic lines (NILs) that differ in a single rice-blast resistance gene are useful as differential varieties in pathogenicity tests and as genetic resources in rice breeding programs. However, because the development and phenotyping process is time-consuming and laborious, such lines exist only for a few genes. Among these is Pi9, a broad-spectrum riceblast resistance gene that encodes a nucleotide-binding siteleucine-rich repeat (NBS-LRR) protein and is part of a multigene family on chromosome 6. This gene is widely used in pyramiding programs for increasing broad-spectrum resistance to M. oryzae, and is more effective than 14 other NILs against several virulent M. oryzae strains (Imam et al., 2014). For example, Pi9+Pita is an effective combination for incorporation in Indian rice varieties (Khanna et al., 2015a). Numerous studies have examined the expression profiles of rice defense response genes (Vergne et al., 2008;Sharma et al., 2016). However, transcriptome profiling studies of rice NILs upon M. oryzae infection are few in number (Sharma et al., 2016). Earlier transcriptome profiling of NIL IRBL22 carrying Pi9 gene upon M. oryzae infection in the background of a japonica cultivar LTH was performed using microarray which provides information only about the known genes (Wei et al., 2013). Here, we have used novel monogenic lines containing Pi9 in the background of Pusa Basmati1 (PB1), a variety released in 1989 as the first high-yielding, semi-dwarf, photoperiod-insensitive, and superior quality scented rice line. This NIL (PB1+Pi9) shows broad-spectrum resistance against 100 diverse strains of M. oryzae from eastern (Variar et al., 2009;Imam et al., 2014) and northern (Khanna et al., 2015b) India. This is the first NIL carrying Pi9 gene in the background of any scented rice and serves an excellent biological material for understanding the molecular basis of rice-Magnaporthe interactions characterized by RNA-seq.
The objectives of this study were to identify early transcriptional changes in the PB1+Pi9 compared with PB1 after 24 h post-infection (hpi) with M. oryzae and to find the unique set of genes that were regulated only in PB1+Pi9 compared to PB1, in providing resistance against M. oryzae in PB1+Pi9. Our results will improve understanding of the molecular pathways and interactions during Pi9-meditated resistance.

Plant Material and Growth Condition
The experiment was performed using resistant NILs of PB1 (O. sativa L. ssp indica) containing Pi9 and its susceptible counterpart line PB1. The seeds of BC3F6 generation of PB1+Pi9 NIL was used. The seeds of both rice lines were surface-sterilized and germinated via soaking in water at 37 • C for 5 days. Seedlings were transferred to pots filled with sterile soil under standard growth conditions of 16 h light (115 µ Mol m −2 s −1 ) and 8 h dark at 25 ± 2 • C. Two-week-old healthy plants (four-leaf stage) were used for inoculation with highly virulent isolate Mo-nwi-53 of M. oryzae.

Fungal Inoculation
Magnaporthe oryzae was maintained on oatmeal agar and Mathur's media at 25 ± 1 • C for 15 days. Conidia of M. oryzae were collected from culture plates via rinsing with 0.25% gelatin. Conidia were filtered with two layers of gauze and their concentration adjusted to 10 5 spores/ml using a hemocytometer. The mock plants were sprayed only with 0.25% gelatin. An atomizer was used for fine spraying so that spore suspension was retained on the leaves. Both lines (PB1+Pi9 and PB1) were phenotyped against several M. oryzae strains collected from eastern and northern India. Finally, Magnaporthe oryzae strain Mo-nwi-53 (from northwest India) was used for inoculation. Three biological replicates were performed for M. oryzae and mock inoculation on both resistant and susceptible lines at fourleaf stage. The inoculated plants were placed in a darkened (covered) humid chamber for 24 h at 25 ± 1 • C and 90% relative humidity. Leaf tissues were collected from fully expanded leaves of each rice line 24 hpi and frozen in liquid nitrogen. Some leaves from each pot of PB1+Pi9 and PB1 were kept for disease development (Mackill and Bonman, 1992). On PB1 leaves, spindle-shaped lesions, characteristic of rice blast, and a disease rating of 5 (on a 0-5 scale) were observed after seven days of incubation. In PB1+Pi9, a hypersensitive response was observed within 2 days of infection ( Figure 1A).

RNA Isolation and Library Preparation for Illumina Hiseq 1000
Total RNA was extracted from leaf tissues using Spectrum TM Plant Total RNA Kit (Sigma) following manufacturer protocol. The quantity and quality of RNA was measured using Nanodrop-1000 (Thermo Fischer Scientific), and quality was further assessed in the Agilent 2100 Bioanalyzer (Agilent Technologies). Samples with RNA Integrity Number (RIN) > 8.5 were used for library preparation. Total extracted RNA samples (5 µg) from each biological replicate were used to isolate poly(A) mRNA and to prepare DNA libraries using TruSeq RNA preparation protocol v.2 (Low Throughput protocol, Illumina, Inc., USA). Library quality and size were assessed with the Agilent 2100 Bioanalyzer (Agilent Technologies) using a high sensitivity DNA Kit. Paired-end sequencing was performed with the TruSeq SBS Kit v3-HS (Illumina, Inc.) on Illumina HiSeq 1000. The CASAVA pipeline 1.7 was used for data processing, de-multiplexing and Bcl conversion. Illumina filters were kept for passing all samples reads. Read quality was checked in FastQC. Low-quality reads and adapters were removed using Trimmomatics (Bolger et al., 2014). The data were deposited as GSE81906.

Bioinformatics Analysis
Raw reads were filtered and checked for sequence contaminants using Trimmomatic and FastQC. To build index of the reference genome (O. sativa japonica) bowtie2-build was used. The splicedread mappers Tophat v2.0.9 (Trapnell et al., 2009) is built on the ultrafast short read mapping program bowtie. Tophat was used to map reads individually for each biological replicate in PB1+Pi9 and PB1 against O. sativa japonica group, cultivar Nipponbare; MSU release 7. During alignment of reads with tophat the following parameters were used read mismatch of 2, read gap length of 2, read edit-distance of 2, splice mismatches 0, minimum intron length of 50, maximum intron length of 500,000 maximum multihits 20, maximum insertion and deletion length of 3. The aligned reads obtained from Tophat were analyzed in Qualimap (García-Alcalde et al., 2012) and visualized on IGV. Qualimap was used to obtain reads mapped to exons (including splice-junction reads mapped to exon ends) for each biological replicate of both lines. Cufflink v2.1.1 was used for assembly into transcripts with reference annotation to guide assembly (Trapnell et al., 2010). The following parameters were used in cufflinks for abundance estimation, average fragment length of 200, fragment length standard deviation of 80 and unlimited alignment allowed per fragment. Cuffdiff was used to quantify transcripts with the merged transcript assembly as well as to detect differentially expressed genes (DEG). Quantification of transcript was done in terms of fragment per kilobase of transcript per Million mapped reads (FPKM), false discovery rate of 0.05 was used for testing and pooled dispersion method was used to estimate dispersion model. Cuffdiff output was analyzed in CummeRbund. The significant differentially expressed loci (SDEL) were identified from all expressed loci after applying multiple corrections (FDR adjusted p ≤ 0.05). The SDEL unique to PB1+Pi9 or PB1 with log 2 fold change ≥2 were used for further analysis (Figure 2). Heat maps were generated using R package, heatmap2. Pearson correlation coefficients were calculated in R. Venn diagrams were generated in InteractiVenn.

Pathway Analysis
Pathway analysis was performed in MapMan version 3.1.1 from the Max Planck Institute of Molecular Plant Physiology, Germany (Thimm et al., 2004). Non-redundant SDEL unique to PB1+Pi9 and PB1 were classified and their functional annotation was visualized via searching against the Oryza sativa TIGR7 database.

Singular Enrichment Analysis
Singular enrichment analysis (SEA) was performed on SDEL unique to PB1+Pi9 and PB1 using AgriGo (Du et al., 2010). Significant GO terms were identified in biological and molecular function category in PB1+Pi9 and PB1.

Co-expression Network Analysis
Co-expression network analysis of proteins corresponding to significant differentially expressed loci (SDEL; FDR adjusted p ≤ 0.05 & log 2 fold change ≥2) unique to PB1+Pi9 in response to biotic stimulus. In the co-expression network, circular node represents the protein encoded by the respective SDEL was performed with STRING 10 (Szklarczyk et al., 2015). Proteins were represented as a node with a specific number. Nodes were connected by interconnecting lines that represent the source by which protein interactions were derived. Source of protein interaction were represented with black, pink, green and blue lines that represent co-expression, experimental data, text mining, and homology respectively. Confidence of each interaction between proteins was obtained in terms of total score.

Real Time PCR Validation
Real time PCR (qRT-PCR) was used to validate SDEL obtained from RNA-Seq. Primers were designed using PrimerQuest (IDT) and were listed in Supplementary Table 1. The ProtoScript M-MuLV First Strand cDNA synthesis kit (NEB) was used to synthesize cDNA. The reaction was performed in the LightCycler R 480 II PCR system (Roche) with a volume of 10 µL, containing 5 µL of SYBR Green I Master, 0.2 µL of forward primer, 0.2 µL reverse primer, 2 µL of 1:5 diluted cDNA template and RNAase free water. Actin was used as internal control. The efficiency of RT-PCR was calculated in both control and target samples. Fold change was calculated using 2 − CT method.

RESULTS AND DISCUSSION
Differentially Expressed Genes in Rice-Blast Resistant NIL (PB1+Pi9) To study early defense mechanisms against M. oryzae in rice, RNA-Seq was performed on PB1+Pi9 NIL and susceptible control PB1 24 hpi. The PB1+Pi9 NIL was developed with marker-assisted backcross breeding (Khanna et al., 2015b). The functional marker NBS2-Pi9 195-1  and genelinked marker AP5659-5 (Fjellstrom et al., 2006) were used to confirm the presence of Pi9 in the NIL. Background selection was performed to minimize linkage drag and maximize recipient parent recovery, which was upto 95.6% (Khanna et al., 2015b). Pi9 is a constitutively expressing gene regardless of pathogen infection; however, post-transcriptional reprogramming allows the R gene to activate downstream processes involved in biotic stress only after M. oryzae infection . Here, the PB1+Pi9-M. oryzae interaction was observed incompatible, possibly mediated by ETI, whereas the PB1-M. oryzae interaction was compatible and likely represented PTI ( Figure 1A).
Among commonly expressed loci, 74% (17,926) of the loci were common between NIL PB1+Pi9 and PB1, with and without pathogen treatment, reflecting basal similarities between their transcriptional profiles, despite differing in Pi9 presence. Wei et al. (2013) reported similar results for the same model system, as did Tao et al. (2003) for Arabidopsis thaliana-Pseudomonas syringae interactions. Together, these data indicate that the O. sativa-M. oryzae interaction is quantitative, like other plantpathogen relationships. Total 8418 (4914 up-regulated and 3504 down-regulated) and 3336 (1837 up-regulated and 1494 downregulated) SDEL were obtained in the PB1+Pi9 and PB1 rice lines 24 hpi, respectively (Supplementary Figure 1). Among them 2027 SDEL were common between PB1+Pi9 and PB1. The number of SDEL with log 2 fold change ≥2 were 1254 (850 upregulated and 404 down-regulated) and 781 (466 up-regulated and 313 down-regulated) in PB1+Pi9 and PB1, respectively ( Figure 1D). Among them, 211 SDEL were found common between resistance and susceptible lines ( Figure 1C). The SDEL exclusively regulated in PB1+Pi9 were 1043 (727 up-regulated and 316 down-regulated) and in PB1 were 568 (341 up-regulated and 227 down-regulated) respectively. The proportion of up and down-regulated unique SDEL was higher in PB1+Pi9 than in PB1. The SDEL with log 2 fold change ≥2 unique to PB1+Pi9 or PB1 were used for further analysis.

Identification of Genes for Primary Defense Response through Respiratory Burst
Plants use ROS as stress-signal transduction molecules, and their accumulation plays a central role in plant stress response (Fujita et al., 2006;Ton et al., 2009), inducing unique ROS-responsive genes (Gadjev et al., 2006). In infected plants, respiratory burst is the earliest and most rapid defense response which involves generation of active ROS (primarily superoxide and H 2 O 2 ) to control pathogen spread (Grant et al., 2000). NADPH oxidase also plays a central role in ROS generation. ROS causes lipid peroxidation and membrane damage (Montillet et al., 2005). Major ROS scavenging enzymes (redoxin, glutathioredoxin, catalase, and peroxidase) have to restrict ROS dependent damage and fine tune ROS signaling (Mittler et al., 2004). Peroxidase catalyze ROS during last step of cell wall fortification to polymerize lignin (Wally and Punja, 2010) and cross-linking of cell-wall components. This makes the cell wall stronger to fight against the invading pathogen. Our results showed differential expression of three categories of respiratory-burst genes, namely, redox state, peroxidases, and glutathione S-transferases (Figure 3 and Supplementary Table 5). Five redoxin genes were found up-regulated in PB1+Pi9 while one is down-regulated. Among seven peroxidase precursors found in PB1+Pi9, six were up-regulated and one is down-regulated. All glutathione S transferase loci were found up-regulated in PB1+Pi9. Previous studies have demonstrated the up-regulation of 10 peroxidase genes in M. oryzae-infected rice (Sasaki et al., 2004) and a 16-fold up-regulation of class III peroxidases in the transgenic rice line TP-Pi54 using microarray . For example, M. grisea infection triggered the differential expression of 34 GST genes in a susceptible rice line (Ribot et al., 2008). Under pathogen attack, many plants differentially regulate GSTs (Alvarez et al., 1998;Wagner et al., 2002). The relative expression of peroxidase (LOC_Os08g02110), thioredoxin (LOC_Os07g29410), and cytochrome P450 (LOC_Os06g39780) was validated with real time PCR for infected and uninfected samples of PB1+Pi9 after 24 hpi (Figure 4, Supplementary Table 6, and Supplementary Figure 2). The qRT-PCR results also support the transcriptome data. Although PB1 shows compatible reaction, many ROS scavenging genes were differentially expressed after M. oryzae infection. Total five redoxin, five peroxidase precursors and three GSTs were up-regulated (Supplementary Figure 3). This may be a result of primary ROS burst after infection as stated by Lamb and Dixon (1997). In PB1+Pi9 and PB1, major loci of redoxins, peroxidases and glutathione S transferases were showing similar trend of up-regulation. This reflects that in both PB1+Pi9 (incompatible interaction) and PB1 (compatible interaction) upon M. oryzae attack, unique set of respiratory burst responsive loci were activated. However, during incompatible interaction FIGURE 4 | The qRT-PCR validation of significant differentially expressed loci (SDEL) in PB1+Pi9 upon M. oryzae infection. Blue bar represents the absolute fold change of PB1+Pi9. Actin was used for transcript normalization. Standard error bar show the standard deviation for three triplicate assay. Abbreviations: WRKY71, WRKY71 transcription factor (LOC_Os02g08440); WRKY76, WRKY76 transcription factor (LOC_Os09g25060); WRKY108, WRKY108 transcription factor (LOC_Os01g60600); bZIP, bZIP transcription factor domain (LOC_Os07g48660); TF HBP-1b, transcription factor HBP-1b (LOC_Os01g06560); AP2, AP2 domain containing protein (LOC_Os02g45420); POX, Peroxidase (LOC_Os08g02110); TRX, thioredoxin (LOC_Os07g29410); DR, disease resistance protein SlVe2 precursor (LOC_Os01g06836); NB/LRR DR, NB-ARC/LRR disease resistance protein (LOC_Os04g43440); PR, pathogenesis-related Bet v I (LOC_Os08g28670); rp3, rp3 protein (LOC_Os12g03080); Cf-2, Cf-2 (LOC_Os01g06876); CHIT, CHIT1-Chitinase (LOC_Os02g39330); OsWAK, OsWAK95-OsWAK receptor-like protein kinase (LOC_Os10g02250); CYP450, cytochrome P450 (LOC_Os06g39780).
ROS accumulation starts with lower amplitude followed by sustained phase of much higher amplitude, while in compatible interaction there is only transient low amplitude single phase ROS accumulation.

Host Hormonal Regulation during M. oryzae Infection
The burst of jasmonic acid (JA) and ethylene (ET) is well characterized in plant-pathogen interactions. Beside JA and ET, salicylic acid (SA) also acts as important signal molecules in biotic stress responses (Pieterse et al., 2009). SA primarily involved in resistance to biotrophic pathogens, whereas JA and ET influence necrotroph resistance (Roberts et al., 2011). In PB1+Pi9, one loci (LOC_Os06g20920) of SAM-dependent carboxyl methyl transferases (SAMT) is upregulated whereas other (LOC_Os06g20790) is down-regulated (Figure 3 and Supplementary Table 7). SAMT converts SA into methyl salicylate (MeSA). During pathogen attack MeSA acts as airborne signal involved in intra and inter plant communication (Shulaev et al., 1997). Earlier it was shown that exogenous SA treatment did not induce M. oryzae resistance in young rice plants (four-leaf stage), but it induces in adults eight-leaf stage (Iwai et al., 2007). Here, plants from PB1+Pi9 and PB1 were infected at the four-leaf stage; therefore, SA did not change much significantly in response to M. oryzae.
JA biosynthesis and signal transduction related six genes were up-regulated in PB1+Pi9 post infection (Figure 3 and Supplementary Table 7). Among six up-regulated loci, four were involved in JA biosynthesis, out of which, two were lipoxygenase (LOX) that converts the Linolenic acid to hydroperoxy derivatives and two were 12-oxophytodienoate reductase (similar to OPR) that converts 12-oxo phytodienoic acid to 3oxo-2-cyclopentane-1-octanoic acid during JA biosynthesis. The rest two up-regulated loci were zinc finger protein (LOC_Os07g42370) involved in JA signal transduction and CBS domain containing membrane protein (LOC_Os09g02710) similar to loss of the timing of ET and JA biosynthesis 2 gene (LEJ2). CBS domain containing protein is reported to regulate thioredoxin system (Yoo et al., 2011).
Ten ET-metabolism-related unique loci were up-regulated and two were down-regulated in infected PB1+Pi9 (Figure 3 and Supplementary Table 7). Among 10 up-regulated loci four were 1-aminocyclopropane-1-carboxylate oxidase proteins similar to ethylene forming enzyme (AT1GO5010), four were AP2 domain containing protein involved in ethylene signal transduction and two were ethylene responsive proteins. The down-regulated loci were AP2 domain containing protein and ethylene responsive protein in PB1+Pi9.
In PB1 only one locus related to JA & ET biosynthesis were found up-regulated (Supplementary Figure 3). The induction of several genes involved in JA and ET signaling and biosynthesis was quantitatively higher in PB1+Pi9 than in PB1. The above trend leads to enhance JA-ET signaling in PB1+Pi9. Previous several reports indicated that JA is a strong activator of resistance against hemibiotrophs M. oryzae and X. oryzae pv. oryzae (Deng et al., 2012;Yamada et al., 2012;Riemann et al., 2013). During transcriptome profiling of resistant (IRBL18, IRBL22) and susceptible (LTH) lines after M. oryzae infection, eight out of 34 enzymes involved in JA biosynthesis were up-regulated (Wei et al., 2013). It was observed that exogenous application of ET inhibitor and generator respectively induced and suppressed riceblast infection (Singh et al., 2004). Activation of ET emission was found little earlier in incompatible compared to compatible interaction (Iwai et al., 2006). Additionally, ET-overproducing rice transformants increased resistance to both the fungal pathogen Rhizoctonia solani and M. oryzae (Helliwell et al., 2013). Normally, ETI response involves redundant activities of both SA and JA-ET pathways (Tsuda et al., 2009). However, when SA signaling is not active, substantial resistance to pathogen is contributed by JA-ET pathway (Dodds and Rathjen, 2010). Our results reflect that the increased level of ET and JA production in PB1+Pi9 helps in providing resistance to M. oryzae compared to PB1.The induction in level of JA-ET hormone in PB1+Pi9 might be controlling the expression of defense genes by regulating the abundance of transcription factors.

Kinase-Mediated Signaling
Plant signal perception and activation of downstream processes is crucial for the innate immunity. Defense signaling involves molecules such as receptor-like kinases (RLKs), MAPK, and calmodulin-related calcium sensor protein. In PB1+Pi9 29 receptor kinases (RK) were found to be differentially expressed (Figure 3 and Supplementary Table 8). Previous research in IRBL18 and IRBL22 (resistant NILs), 103 receptor kinases (subfamily RLK and WAK) were up-regulated, while one was down-regulated (Wei et al., 2013). In Digu rice, 48 SDEL of receptor kinases were found upon M. oryzae infection (Li et al., 2016). Further, in a transgenic rice line carrying rice-blast resistance gene Pi54 , as well as resistant NILs IRBL18 and IRBL22 (Wei et al., 2013), signaling-related genes were highly up-regulated post M. oryzae infection, compared with susceptible controls.
During infection, WAK triggers innate immune response through detecting fungal cell wall-associated oligogalacturonides (Brutus et al., 2010). From OsWAK subfamily four OsWAK (LOC_Os02g56370, LOC_Os09g38850, LOC_Os09g38910, LOC_Os10g02250) and one OsMAPK (LOC_Os05g49140) were up-regulated in PB1+Pi9 (Figure 3 and Supplementary  Table 8). OsWAK receptor-like protein kinase (OsWAK95, LOC_Os10g02250) was validated with real time PCR in infected and uninfected samples of PB1+Pi9 (Figure 4, Supplementary Table 6, Supplementary Figure 2). Upon M. oryzae infection several calcium dependent kinase were also activated. In PB1+Pi9 eight calmodulin related calcium sensor protein and three G-proteins were up-regulated. The most highly up-regulated receptor kinase in PB1+Pi9 after M. oryzae infection is LOC_Os06g13320.1 with FPKM value of 0.6 (3.5-fold; Figure 3 and Supplementary Table 8). A previous microarray analysis of blast-infected rice revealed high upregulation of OsWAK71 and OsWAK25 (Wei et al., 2013). During transcriptional profiling of resistant vs. susceptible lines, the former exhibited more up-regulated WAKs (Bagnaresi et al., 2012). Mitogen-activated protein kinases are conserved signaling molecules that transduce extracellular stimuli into intra-cellular responses. Active MAPK further activates the downstream transcription factors (TF). Earlier up-regulation of four MAPK transcripts, as well as MAP3K.3 and MAP3K.1 isoforms, were observed only in a resistant rice line (Bagnaresi et al., 2012). In PB1 few receptor kinases were differentially expressed (Supplementary Figure 3). But these kinases were different from kinases differentially expressed in PB1+Pi9. So, the signaling mediated by kinases in PB1 is not able to provide resistance against M. oryzae. The expression trend of kinases in mediating signaling to provide resistance against M. oryzae in PB1+Pi9 reflects that up-regulation of several subfamilies of kinases, calcium sensor and G protein that positively regulate the activation of transcription factors. Few kinases being down-regulated negatively regulate the defense response in PB1+Pi9. This also shows that the positive and negative regulation by kinases being differentially expressed in PB1+Pi9 together mediated stronger ETI-responses in PB1+Pi9 compared to PB1, with a basal response similar to those of earlier studies.
The ERF transcription factors integrate signals received from JA and ET (Lorenzo et al., 2003).

Activation of Primary and Secondary Metabolites
Carbohydrate, lipid, and protein metabolism, as well as photosynthesis, produce primary metabolites in plants. In contrast, secondary metabolites include phenylpropanoids, lignin, phenolics, waxes, terpens, and flavanoids. Pathogen infection affects host lipid metabolism. We observed the upregulation of all six loci involved in phospholipid-biosynthesis in infected PB1+Pi9. Among six loci four (LOC_Os01g50030, LOC_Os01g50032, LOC_Os05g47540, LOC_Os05g47545) were similar to phosphoethanolamine-N-methyltransferase that converts S-adenosyl methionine to S-methionine homocysteine during phospholipid biosynthesis (Figure 5 and Supplementary  Table 10). Rest two up-regulated loci were phosphatidate cytidylyltransferase (LOC_Os10g17990) and diacylglycerol kinase (LOC_Os04g54200) involved in phospholipid biosynthesis (Li-Beisson et al., 2013). Phosphatidic acid (PA) formed during phopsholipid synthesis induce generation of ROS in Arabidopsis (Sang et al., 2001) and rice (Yamaguchi et al., 2003). PA activates MAPK and induces expression of defense related genes (Yamaguchi et al., 2005) and plays positive role in mobilizing defense response (Yamaguchi et al., 2009). A locus of sphingolipid delta desaturase (LOC_Os02g42660) was found two-fold up-regulated in PB1+Pi9. Sphingolipids induce ROS production that leads to programmed cell death (Brodersen et al., 2002;Shi et al., 2007). Genes involved in fatty acid (FA) synthesis and elongation were also highly up-regulated in PB1+Pi9 (Figure 5 and Supplementary Table 10). Four loci involved in FA synthesis and elongation were found up-regulated in PB1+Pi9. Among four up-regulated loci two were ketoacyl CoA synthase involved in both long chain FA and wax biosynthesis (Todd et al., 1999). Rest two were 3-oxoacyl synthase and acyl desturase. Fatty acid plays role in regulating enzyme activity that were involved in generation of signal molecules in plant defense (Shah, 2005). Two loci involved in fatty acid desaturation namely omega-3 and omega-6 fatty acid desaturase (LOC_Os07g23430, LOC_Os03g18070) were up-regulated in PB1+Pi9. These enzymes were responsible for the synthesis of polyunsaturated fatty acids (PUFA), which act as antimicrobial agents and signaling molecules (Iba, 2002;Turner et al., 2002;Weber, 2002;Yaeno et al., 2004). Oxylipins, another antimicrobial compound, were synthesized from PUFA through the action of lipoxygenase. Two loci of lipoxygenase (LOC_Os08g39840, LOC_Os04g37430) were up-regulated in PB1+Pi9 (Table 1). In plants, oxylipins were potent signaling molecule in defense response. Five loci of lipases were found up-regulated in PB1+Pi9, among which three loci (LOC_Os08g04800, LOC_Os04g56240, LOC_Os07g34420) were lipases and two were phospholipases (LOC_Os06g40170, LOC_Os05g07880) (Figure 5 and Supplementary Table 10). Phospholipases catalyze hydrolysis of phospholipid to release free fatty acid for synthesis of JA during plant defense response (Shah, 2005). Infection also alters protein degradation and modification. The consecutive action of three protein classes (E1s, E2s, E3s; Komander, 2009) mediates the process of ubiquitination. PB1+Pi9 exhibited differential regulation of genes related to F-Box and RING subcomplex of E3 ligase in ubiquitin-dependent degradation (Figure 3 and Supplementary Table 11). The present results were in agreement with previous data in which it has been found that a UPS ubiquitin proteasome system protein OsBBI1 with E3 ligase activity possesses broad-spectrum resistance to M. oryzae . JA-ET hormones induced in PB1+Pi9 may be regulating expression of defense genes by TF through regulated protein degradation. Photosynthesis in rice is compromised upon M. oryzae infection because the pathogen competes with the host for photosynthates (Bastiaans and Kropff, 1993). Several studies reported down-regulation of photosynthesis-related genes during M. oryzae attack (Vergne et al., 2007;Bagnaresi et al., 2012). Many genes from the light reaction of photosynthesis were differentially expressed in PB1+Pi9 (Figure 5). Similarly, an earlier study indicated that the down-regulation of photosynthesis-related genes reflects the usage of energy and resources to defend against invading pathogens (Bolton, 2009).
Secondary metabolites were the products of specialized metabolic pathways that were important in defense but not essential to plant existence. Several enzymes like phenyl ammonia lyase (PAL), coumarate CoA ligase (4CL), cinnamyl alcohol dehydrogenase (CAD) and caffeoyl CoA-o-methyl transferase (CCoAOMT) involved in phenyl propanoid (PP) pathway were found up-regulated in PB1+Pi9 upon M. oryzae attack ( Table 2). Two loci, each of PAL and CCoAOMT that convert phenylalanine to cinnamic acid and caffeoyl-CoA to feruloyl-CoA respectively were found up-regulated in PB1+Pi9. Two loci of AMP binding domain were found upregulated that is similar to 4-coumarate CoA ligase (4CL) and were involved in several conversions during PP pathway. In PB1+Pi9 two loci of dehydrogenase similar to cinnamyl alcohol dehydrogenase (CAD7) were found up-regulated. And CAD enzyme converts aldehyde to alcohol during PP pathway. Peroxidases found up-regulated in PB1+Pi9 were involved in conversion of alcohol to lignin during PP pathway (Figure 5 and Supplementary Table 12). The phenylpropanoid pathway is critical to plant defense because it is involved in synthesis of phytoalexin (include isoflavonoids, terpenoids, alkaloids etc), lignin, flavanoids, coumarins, phenylpropanoid esters, and cutin synthesis. Phytoalexin has antimicrobial activity and lignin acts as a physical barrier against pathogens Dixon et al., 2002). During pathogen attack, flavonoid biosynthesis and accumulation is enhanced (Treutter, 2005). In a previous study of resistant NIL, 19 out of 80 enzymes related to phenylpropanoid biosynthesis and 6 out of 14 genes related to shikimate biosynthesis were up-regulated (Wei et al., 2013). Additionally, secondary metabolites were highly enriched in resistant IRBL18 and IRBL22 compared with a susceptible control line after M. oryzae infection (Wei et al., 2013). Plant susceptibility increases with down-regulation of phenylpropanoid pathway genes (Bhuiyan et al., 2009;Naoumkina et al., 2010) that were regulated mainly by MYB transcription factors Zhao and Dixon, 2010). As already described, MYBs were highly up-regulated in PB1+Pi9, correlating with phenylpropanoid biosynthesis. In PB1+Pi9 peroxidase which catalyze phytoalexin and lignin biosynthesis were also found highly up-regulated. The upregulated loci involved in primary and secondary metabolism were different in PB1+Pi9 from PB1 with variable magnitudes. In PB1, very few loci involved in primary and secondary metabolite synthesis and degradation were altered upon M. oryzae infection (Supplementary Figure 4 and Table 2). Only two loci involved in PP pathway namely AMP domain binding protein (similar to 4CL enzyme) and dehydrogenase (similar to CAD9) were found up-regulated in PB1. This shows that during primary and secondary metabolism, metabolites formed in PB1 were not able to provide resistance while in PB1+Pi9 the higher activation of enzymes involved in primary and secondary metabolism for example phospholipid and fatty acid biosynthesis that induce defense genes expression by activating MAPK and higher accumulation of antimicrobial compounds like PUFA, oxylipin, phytoalexin and phenol also prevent growth of M. oryzae in PB1+Pi9.
Pathogenesis-related proteins were induced in response to pathogen infection and thus often used as markers to identify plant defense. The PR Bet v I family encodes genes homologous to PR10 (Radauer et al., 2008). In infected PB1+Pi9 two PR Bet v I family proteins (LOC_Os08g28670 and LOC_Os04g50700) were highly up-regulated. In addition, two disease resistance proteins (LOC_Os04g43440 & LOC_Os01g06836) were upregulated in PB1+Pi9 (Figure 3 and Table 1). Disease resistance protein SlVe2 precursor (LOC_Os01g06836), pathogenesisrelated Bet v I family protein (LOC_Os08g28670), NB-ARC/LRR disease resistance protein (LOC_Os04g43440), rp3 protein (LOC_Os12g03080) and Cf-2 (LOC_Os01g06876) were validated with real time PCR in infected and uninfected samples of PB1+Pi9 after M. oryzae infection (Figure 4, Supplementary Table 6, and Supplementary Figure 2). During infection, PR-10 proteins were induced in various plant species and exhibit ribonuclease activity (Mcgee et al., 2001;Kim et al., 2004Kim et al., , 2008. PR genes were shown to be expressed in the resistant Digu line compared with the susceptible LTH line (Li et al., 2016). Lipid transfer proteins (LTP) from the PR-14 family exhibit antifungal activity. PB1+Pi9 exhibited up-regulation of nine LTP loci 24 hpi ( Table 1). Earlier it was shown that differential expression of rice LTP in response to M. grisea was higher in incompatible compared with compatible interaction (Broekaert et al., 1997;van-Loon and van-Strien, 1999;Kim et al., 2006). In susceptible line (PB1) none of the chitinases, β-1,3-glucanases and LTP showed up-regulation (Supplementary Figure 3). In PB1+Pi9, the activation of PR and disease resistance proteins leads to synthesis of antimicrobial compounds (e.g., secondary metabolites) that arrest pathogen growth. These results showed expanded transcriptional activation of different metabolic pathways and PR genes during the early response of rice to M. oryzae. Furthermore, it indicates that downstream transcriptional regulation may be controlled by Pi9, which provides resistance to PB1 against M. oryzae.

Singular Enrichment Analysis
Singular enrichment analysis (SEA) was performed to understand the biological processes and molecular functions of genes which plays important role in Pi9-mediated resistance in PB1+Pi9. Enriched biological processes and molecular function were found for PB1+Pi9 and PB1. Eight biological processes (BP) and four molecular functions (MF) were found significantly (corrected p < 0.01) enriched among 1043 SDEL uniquely found in PB1+Pi9 (Figure 6 and Table 3). Among eight biological processes 151, 216, and 69 SDEL belonging to  the response to stress, response to stimulus and response to endogenous stimulus GO term, respectively. These results reflect that large numbers of loci were being differentially regulated in response to stress in PB1+Pi9 upon M. oryzae infection. During pathway analysis in PB1+Pi9, significant enrichment of up-regulated genes is observed in both primary metabolism and secondary metabolism. The similar trend was reflected in SE analysis in PB1+Pi9 with 489 and 30 SDEL falling in primary and secondary metabolic processes respectively. In PB1+Pi9, 16 SDEL were falling under photosynthesis GO term. This shows that photosynthesis related genes were affected upon M. oryzae attack because both host and pathogen were fighting for the photosynthates. Among three molecular functions in PB1+Pi9 77, 24, 77, and 317 SDEL were found in transcriptional regulation activity, oxygen binding, transcription factor activity and catalytic activity respectively. Overall, PB1+Pi9 have several unique SDEL involved in transcriptional regulation by different transcription factors suggesting that resistance involves the action of many genes regulating different biological processes. Pathway analysis demonstrated higher transcriptional regulation in PB1+Pi9 compared with PB1, corresponding results were observed in SE analysis (Supplementary Figure 5 and Supplementary Table 13).
In the present study, these data provide a higher degree of confidence regarding the regulatory role of a single riceblast resistance gene present in PB1+Pi9. The GO term directly associated with defense response to fungal infection was "response to biotic stimulus" (GO:0009607) was found significantly enriched only in PB1+Pi9. In PB1+Pi9, SDELs involved in biotic stimulus were highly up-regulated and they FIGURE 7 | Heat-map of significant differentially expressed loci (FDR adjusted p ≤ 0.05 & log 2 fold change ≥2) unique to PB1+Pi9 present in the gene ontology term, response to biotic stimulus (GO:0009607) and their respective log 2 fold change in both PB1+Pi9 and PB1. Red represents up-regulated loci and green represents down-regulated loci.
include transcriptions factors, glycosyl hydrolase, cytochrome P450, zinc finger proteins, lipoxygenase, disease resistance proteins and pathogensis related proteins (Figure 7 and Supplementary Table 14) These SDEL were important candidates for blast resistance in the rice-Magnaporthe pathosystem, and network analysis was thus performed to understand their co-regulation.
Co-expression Network of SDEL Unique to PB1+Pi9 in Response to Biotic Stimulus SDEL identified in response to biotic stimulus (GO:0009607) during enrichment analysis were significantly up-regulated only in PB1+Pi9. The co-expression network analysis of these important candidate SDEL helps in understanding Pi9 mediated resistance. In the network, each node represents a protein encoded by a SDEL (Figure 8 and Supplementary Table 15). FIGURE 8 | Co-expression network of proteins corresponding to significant differentially expressed loci (SDEL; FDR adjusted p ≤ 0.05 & log 2 fold change ≥2) unique to PB1+Pi9 in response to biotic stimulus. In the network, circular node represents the protein encoded by the respective SDEL and the interconnecting lines represent the source by which protein interactions are derived. The protein interactions are derived from coexpression, concurrence, experiments, neighborhood, database, text mining, and homology sources. All the interactions represented have score greater than 0.4 that shows medium confidence for the network.
Several interactions in the network were derived from coexpression and text mining source. The total score of each interaction in the network is above 0.4, so the overall network is of medium confidence. Several interactions in the network have total score above 0.7 so they represent the high confidence interactions (Supplementary Table 15).
The network revealed co-expression of three zinc finger proteins 4335698 (LOC_Os04g32480), 4347164 (LOC_Os09g26780), and 4331834 (LOC_Os03g08330) from the STRING database. These zinc finger proteins were involved in the KEGG plant-pathogen interaction pathway (osa04626) and plant hormone signal transduction pathway (osa04075) of rice (Table 4). In the network (supported by experimental data), three WRKY proteins 4328512 (OsWRKY28), 4347069 (OsWRKY71), 4341678 (OsWRKY76) were co-expressed. These WRKYs have zinc finger domains and show homology across related species. The orthologous genes for OsWRKY71 and OsWRKY76 in A. thaliana is AtWRKY40, whereas OsWRKY28 has two orthologous genes, AtWRKY60 and AtWRKY18. The STRING database report showed that WRKY40 and WRKY18 possessed confidence co-expression scores of 0.695; both also interact with the W box, an elicitor responsive cis-acting element. Functional and physical interaction between WRKY40, WRKY18, and WRKY60 has been reported in A. thaliana (Xu et al., 2006). These three WRKYs were structurally similar and induced by pathogen attack. Triple mutants of all three WRKYs (WRKY40, WRKY18, WRKY60) were resistant to P. syringae compared with wild type (Xu et al., 2006). Thus, network analysis showed that in PB1+Pi9 NIL, OsWRKY71, OsWRKY76, OsWRKY28 were co-expressed with zinc finger proteins (4335698,4347164,4331834). Proteins 4345717 (LOC_Os08g34790) and 4328485 (LOC_Os02g08100) involved ). This co-expression network of genes involved in biotic stimulus improves our understanding of genome-wide co-expression and suggest protein-level interactions among genes in PB1+Pi9 NIL to provide Pi9-mediated resistance.

Pi9-Mediated Resistance in PB1+Pi9
The present study confirms that a single functional blastresistant gene (Pi9) in PB1+Pi9 activates a cascade of defense response genes, leading to incompatible interaction between host and pathogen. The resistant NIL displays a broad spectrum of transcriptional changes upon M. oryzae attack. The pathway analysis of SDEL revealed that genes involved in cell wall fortification, respiratory burst, kinase signaling, hormone signaling, WRKY, MYB, and ERF transcription factors, defense response genes (peroxidases, glucanases, chitinases, laccases, lipoxygenase, phenyl ammonia lyase, PR proteins) were highly activated in PB1+Pi9 than in PB1. The activation of defense response genes results in the synthesis of antimicrobial secondary metabolites, inhibiting the spread of M. oryzae in PB1+Pi9. We also proposed a model showing how Pi9-mediated regulators controls incompatible interaction against M. oryzae infection (Figure 9).
The role of WRKYs in the O. sativa-M. oryzae interaction had been previously reported in several studies. In the present study, MYB, bZIP and ERF were also observed to mediate downstream resistance response 24 hpi, along with WRKY genes. Increased expression of JA/ET signaling genes in PB1+Pi9 showed that M. oryzae is regulated via a pathway between the typical SA regulation of biotrophs and JA regulation of necrotrophs. The pathway analysis reflects the increased complexity of the PB1+Pi9 defense response, compared with that of PB1, the former involving intricate mechanisms of cell wall fortification and ROS mediated signaling upon M. oryzae attack. The effector molecule released from M. oryzae bind to Pi9 gene and activates the downstream signaling by kinases and hormones that activates the transcription factors during the initial infection phase. These transcription factors activate the transcription of defense response genes and formation of secondary metabolites. The secondary metabolites like phytoalexin and phenols act as antimicrobial compounds and stops M. oryzae growth in PB1+Pi9, unlike PB1. Singular enrichment analysis of SDEL also supported the pathway analysis data, indicating that the genes involved in signaling, response to biotic stimulus, biological regulation, transcriptional regulation, as well as primary and secondary metabolism were enriched in PB1+Pi9 compared with PB1.
Earlier reports using rice NIL to study rice-M. oryzae interaction mainly used microarray technology (Sharma et al., 2016). The transcriptional profiling using microarray has limited range of expression which identifies only known genes, while RNA sequencing has greater dynamic range of detection along with identification of novel candidate genes (Wang et al., 2009;Agarwal et al., 2010;Xu et al., 2012). The level of capture in transcriptome data is clearly reflected by comparison of upregulated genes found in different studies on rice trancriptome 24 hpi with M. oryzae. It has been observed that numbers of differentially expressed genes were higher in both compatible and incompatible interactions using RNA-seq compared to microarray. In PB1+Pi9 and PB1, the number of SDEL were 1254 and 779 (log 2 fold change ≥2; FDR adjusted p ≤ 0.05), respectively, while the number of SDEL in IRBL22 and LTH 24 hpi upon M. oryzae infection were 649 and 131 (absolute fold change ≥2; p < 0.05), respectively (Wei et al., 2013). This is probably due to the sensitivity of technology to capture all transcripts, differences in the background of both the NILs and also the strain of M. oryzae used for infection. In our study in-depth transcriptional level changes covered signaling mediated by receptor kinases, wall associated kinases, calmodulin related calcium sensor proteins and Gproteins; hormone signaling mediated by JA/ET; MYB, bZIP, and ERF transcription factors along with WRKY; activation of enzymes involved in lipid metabolism like phospholipid and fatty acid biosynthesis; activation of defense response genes and accumulation of antimicrobial compounds like PUFA, oxylipin and phenol that prevent the growth of M. oryzae.

CONCLUSIONS
The present study deciphered the transcriptional snapshot of blast resistant NIL PB1+Pi9 at 24 hpi with M. oryzae. In this line, singular enrichment analysis showed that SDELs involved in biotic stimulus were highly up-regulated. The co-expression network of proteins involved in biotic stimulus (GO:0009607) played an important role in understanding broad-spectrum resistance against M. oryzae in PB1+Pi9, as these genes were down-regulated in PB1. Several proteins in the network were involved in both plant-pathogen interaction and plant hormone signal transduction pathways. The support from both transcriptional and co-expression interaction data indicates that these are prominent candidates for blast resistance in the rice-Magnapothe pathosystem. These important candidates can be manipulated to modify important pathway and enhance disease resistance. Thus, current study revealed additional transcriptional changes via modulating genome-wide transcriptional regulation using three important methods namely pathway, SEA and co-expression analysis, along with qRT-PCR validation of important candidates in rice NIL (PB1+Pi9), 24 hpi upon M. oryzae infection and helped in unraveling mechanism of broad-spectrum blast resistance mediated by Pi9 gene.