Comparative transcriptome analysis by Tilletia horrida infection in resistant and susceptible rice (Oryza sativa L.) male sterile lines reveals candidate genes and resistant mechanisms

Background Rice kernel smut (RKS), caused by the basidiomycete fungus Tilletia horrida , is one of the most devastating diseases affecting the production of male sterile lines of rice ( Oryza sativa ) worldwide. However, molecular mechanisms of resistance to T. horrida have not yet been explored.Results In the present study, analysis of the amount of T. horrida biomass in rice kernels and RNA sequencing (RNA-Seq) analysis of rice male sterile lines resistant and susceptible to RKS (Jiangcheng 3A and 9311A, respectively) were conducted after T. horrida infection. Transcriptomic analysis showed that 9, 840 and 7, 902 differentially expressed genes (DEGs) were observed in Jiangcheng3A and 9311A, respectively, after T. horrida inoculation, compared with the control. KEGG analyses of DEGs revealed that cutin, suberine and wax biosynthesis, avonoid biosynthesis, glutathione metabolism, tyrosine metabolism, and biosynthesis of unsaturated fatty acids were enriched by T. horrida inoculation in Jiangcheng 3A; however, not signicantly enriched in 9311A. Furthermore, the DEGs related to plant-pathogen interaction, plant hormone signal transduction, and transcriptional factor genes were identied in both rice male sterile lines.Conclusions This is the rst comparative transcriptome analysis of two rice genotypes with different responses to T. horrida infection, revealing DEGs with potentially important roles in defense against T. horrida infection. Ours results will serve as a strong foundation for developing effective strategies for T. horrida -resistance breeding. the resistance phenotype of Jiangcheng 3A against T. horrida. results of this will enhance our understanding of the basal differences in gene expression between these resistance and susceptible rice genotypes. This is the rst report of rice resistance mechanisms effective against T. horrida using comparative transcriptome proling.

technology for transcriptome sequencing, and its high-throughput ability has accelerated genomic research and gene function analysis [17]. Many genes related to pathogen resistance and stress tolerance have been identi ed in rice, wheat, maize, and Arabidopsis thaliana via transcriptome analysis [18,19,20,21]. Additionally, numerous genes imparting resistance to rice sheath blight have been identi ed using transcriptome sequencing, and the associated resistance mechanisms have been explored [22].
Transcriptome analysis has also revealed the molecular mechanisms of drought resistance in rice [23].
In this study, we used RNA-Seq to identify rice genes showing signi cant changes in expression in response to T. horrida infection during the early stages of infection. Rice genotypes used in this study included Jiangcheng 3A and 9311A, which are resistant and susceptible to T. horrida, respectively. Gene expression in these lines was analysed at 8, 12,24,48, and 72 h after T. horrida inoculation using Illumina RNA-Seq analysis, followed by quantitative real-time PCR (qRT-PCR) to identify differentially expressed genes (DEGs) potentially responsible for the resistance phenotype of Jiangcheng 3A against T. horrida. The results of this study will enhance our understanding of the basal differences in gene expression between these resistance and susceptible rice genotypes. This is the rst report of rice resistance mechanisms effective against T. horrida using comparative transcriptome pro ling.

Plant growth conditions and pathogen inoculation
The resistance phenotype of Jiangcheng 3A to T. horrida infection has been previously identi ed via inoculation tests from 2014-2016 [15], 9311A is a considered a highly susceptible standard. They were provided by the Department of Food Crop Research Institute of Hubei Academy of Agricultural Science.
The T. horrida isolate JY-521, a virulent strain with a sequenced the genome, was used for plant inoculations and separating and preserving by the Department of Rice Research Institute of Sichuan Agricultural University [14]. Monoconidial cultures of T. horrida were grown in potato dextrose broth in an incubator at 28℃ for 7 days on a shaker at 150 rpm. Approximately 1 mL of conidial suspension (10 6 conidia mL -1 ) was injected into the immature kernels of eld grown rice plants at the booting stage, i.e. 5-7 days before heading. Inoculations were performed in the late afternoon to facilitate infection. The biomass of T. horrida spores in rice kernels was quanti ed to determine the sampling time point. Panicles were harvested at 8, 12,24,48, and 72 h post inoculation. Kernels of each cultivar inoculated with 1 mL of sterilized distilled water for 8 h served as a control (infection 0 h). Three biological replicates of each sample were collected, immediately frozen in liquid nitrogen, and stored at -80℃ for RNA isolation.
qRT-PCR analysis First strand cDNA was synthesized from total RNA using the Transcriptor First-Strand cDNA Synthesis Kit (Roche, Indianapolis, IN, USA). The cDNA samples were then subjected to qRT-PCR on a Bio-Rad CFX96 Real-Time PCR System (Foster City, CA, USA), according to the manufacturer's instructions. The PCR reactions were prepared in a 20µL volume, containing 3 µL cDNA and 0.8 µL each of the forward and reverse gene speci c primers. Each PCR was replicated three times. The Ubiquitin (UBQ) gene was used as an internal control for data normalization. Gene expression levels were calculated using the 2 -∆∆Ct method. The biomass of T. horrida spores in rice kernels was quanti ed based on the abundance 18S rRNA of T. horrida relative to that of rice [24]. Primers used for qRT-PCR are listed in Additional le 1: Table  S1.
RNA-Seq library preparation, Illumina sequencing, and read-mapping RNA-Seq libraries were constructed using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA), according to the manufacturer's instructions, and unique index codes were added to each sample. RNA samples were sent to Beijing Novogene Biological Technology Co., Ltd for cDNA library construction and Illumina sequencing (HiSeq TM 2500, San Diego, NEB, USA). Sequencing was performed using Illumina Hiseq platform to generate 125 bp paired-end reads. Reads with low quality score and those containing adaptor sequences and stretches of -Ns were removed from the raw data. The reference genome of Nipponbare rice and gene model annotation les (Rice Annotation Project) were downloaded directly from the rice genome website (ftp://ftp.ensemblgenomes. org/ pub/ plants/ release_36/ fasta/ oryza_indica/ dna/ ). Index of the reference genome was built using Bowtie v2.2.3, and paired-end reads were aligned to the reference genome using TopHat v2.0.12 [25,26,27]. The number of reads mapped to each gene were counted using HTSeq v0.6.1 [28]. Then, the number of fragments per kilobase of transcript sequence per million (FPKM) of each gene was calculated based on the length of the gene and number of read counts mapped to that gene.

Differential expression analysis
Differential expression analysis of two conditions (three biological replicates per condition) was performed in R using DESeq package (v1.18.0) [29]. The DESeq method provides statistical routines for determining differential expression in digital gene expression data using a model based on negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). Genes with adjusted P <0.05 were designated as differentially expressed [30].

Series-clustered analysis
A total of 4729 DEGs co-expressed in the two rice male sterile lines were classi ed into 24 clustered pro les based on the trend observed in gene expression using the Short Time-series Expression Miner (STEM) software [31]. The clustered pro les with P≤ 0.05 were considered as statistically signi cant.
Then, DEGs in all or each pro le were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. Based on the hypothesis test of P-value calculation and FDR correction [30], the GO terms or pathways with q-value ≤ 0.05 were de ned as signi cantly enriched.
GO and KEGG enrichment analysis of DEGs GO enrichment analysis of the DEGs was performed in R using the GOseq package, after correcting for gene length bias [32,33]. GO terms with adjusted P < 0.05 were considered signi cantly enriched. The KEGG database enables understanding the high-level functions and utilities of biological systems, such as the cell, organism and ecosystem, using molecular information, especially large-scale molecular datasets generated by genome sequencing and other high throughput experimental technologies (http://www.genome.jp/kegg/). The statistical enrichment analysis of DEGs in KEGG pathways was performed The using the KOBAS software [34,35].

Results
Quanti cation of T. horrida biomass in resistant and susceptible rice male sterile lines We have previously shown that the rice male sterile line Jiangcheng 3A is resistant to T. horrida, whereas the rice male sterile cultivar 9311A is susceptible [15]. To con rm the response of Jiangcheng 3A and 9311A to T. horrida infection, we infected both these cultivars with T. horrida at the booting stage (3-5 days before heading) and quanti ed pathogen growth in immature kernels over an infection time point.
After 31 days of infection, 95% of 9311A panicles contained soft and hollow grains covered with a dark powdery mass of spores, and an average of 17.3 grains were infected per panicle (Fig. 1a, c). By contrast, only 3% of Jiangcheng 3A grains were diseased, and only 2.1 kernels were infected per panicle. Thus, lower disease incidence observed in the rice male sterile cultivar Jiangcheng 3A con rmed that it is resistant to T. horrida infection (Fig. 1b, c). To further reveal differences infection between the two rice male sterile cultivars in response to T. horrida infection and determine the ideal time points for transcriptome sequencing, we quanti ed the relative amount of fungal biomass in immature kernels of Jiangcheng 3A and 9311A using qRT-PCR during early stages of infection. Spores of T. horrida were detected in both cultivars at 8 h after infection. However, the biomass of T. horrida in immature kernels of 9311A continued to increase at a higher rate than in immature kernels of Jiangcheng 3A; at 48 h, the T. horrida biomass in 9311A kernels was ~ 2-fold greater than that in Jiangcheng 3A kernels (Fig. 1d).

RNA-Seq analysis of rice male sterile lines
To investigate changes in gene expression levels during the early infection stages, we conducted RNA-Seq analysis of Jiangcheng 3A and 9311A kernels at 8, 12, 24, 48, and 72 h post inoculation, with three biological replicates per time point. RNA-Seq analysis of all 36 samples generated 294.78 GB data, and 1, 965, 147, 576 clean reads were obtained after quality trimming and adaptor removal (Additional le 2: Table S2). Across all 36 samples, the ratio of total mapped reads varied from 76.68-88.62%, and that of unique mapped reads varied from 75.62-87.65% (Additional le 3: Table S3). We performed correlation analysis of expression levels among the three replicates of each sample to assess the reliability of RNA-Seq data. Heat map of the correlation coe cients (R 2 ) of all 36 transcriptome samples shown in Additional le 4: Fig. S1. Overall, a high correlation (R 2 >0.8) was detected among biological replicates of each RNA-Seq sample at each time point.
Differential expression of genes in T. horrida inoculated Jiangcheng 3A and 9311A kernels To identify genes involved in resistance to T. horrida, we compared changes in gene expression between resistance male sterile line Jiangcheng 3A and susceptible male sterile line 9311A at different time points.
We detected a total of 17, 742 DEGs in both the rice male sterile lines of which 9, 840 and 7, 902 DEGs were observed in Jiangcheng3A and 9311A, respectively (Fig. 2a). Among these DEGs, 4, 729 DEGs were shared by the two rice male sterile lines at different time points (Additional le 5: Table S4), and comparing these two groups, at 8, 12, 24, 48, and 72 h post inoculation, there are 332, 992, 82, 852, and 4, 246 DEGs were shared in both cultivars, respectively ( Fig. 2b-f). Furthermore, following T. horrida inoculation, the number of up-regulated genes in Jiangcheng 3A kernels was signi cantly higher than that in 9311A kernels at the early infection time points. For example, at 8 h, 2,289 DEGs (1,251 up-regulated, and 1,038 down-regulated) were detected in Jiangcheng 3A kernels, compared with 1,617 DEGs (52 upregulated, and 1,565 down-regulated) in 9311A kernels (Fig. 2a). The number of DEGs in Jiangcheng 3A kernels was the highest at 72h post inoculation (9,026 DEGs; 4,226 up-regulated, and 4,800 downregulated), which was higher than the number of DEGs in 9311A kernels at the same time point (5,737 DEGs, 2, 275 up-regulated, and 3, 462 down-regulated) (Fig.2a). We further performed the trend analysis of 4, 729 common DEGs in two male sterile lines, and these genes were signi cantly enriched in 24 and 22 pro les in 9311A and Jiangcheng 3A kernels, respectively (Fig. 3a, c). The number of DEGs in pro le 23 (sustained up regulation) and pro le 0 (sustained down regulation) was most in both male sterile lines ( Fig. 3a, c).

Annotation and function classi cation of DEGs in two male sterile lines
To identify whether these DEGs in both male sterile lines were signi cantly associated with a speci c biological process, molecular function, or cellular component, GO analyses of DEGs identi ed in each group (K8vsK0, K12vsK0, K24vsK0, K48vsK0, and K72vsK0; G8vsG0, G12vsG0, G24vsG0, G48vsG0, and G72vsG0) were performed, and the top 30 GO enrichment terms were listed (Additional le 6 and 7: Table   S5 and S6). For cellular component category, 35 GO terms were enriched at ve inoculation time points of Jiangcheng 3A, while 27 GO terms were enriched in 9311A. Amongst these, 21 GO terms were enriched in both the genotypes, 14 GO terms namely DNA packaging complex, nuclear chromatin, protein-DNA complex, nucleosome, nuclear chromosome part, cytosolic ribosome, ribosome, ribosomal subunit, cytosolic part, non-membrane-bounded organelle, intracellular non-membrane-bounded organelle, ribonucleoprotein complex, large ribosomal subunit, and cytosolic large ribosomal subunit were unique enriched in Jiangcheng 3A. The molecular function category was enriched with 17 and 24 GO terms in at least any of the inoculation time points of Jiangcheng 3A and 9311A, respectively. Both the genotypes showed 12 common GO terms enrichment during the T. horrida inoculation. Carbohydrate binding, metal ion binding, urea transmembrane transporter activity, structural constituent of ribosome, and structural molecule activity molecular function were only enriched during the inoculation period in Jiangcheng 3A. In biological process category, 20 GO terms, involve in polysaccharide metabolic process, cell wall modi cation, photosynthesis, pectin catabolic process, galacturonan metabolic process were enriched in both male sterile lines.
We further identi ed speci c function of 4,729 common DEGs using GO enrichment analysis. For DGEs in pro le 23, response to stimulus, developmental process, growth, and immune system process were enriched in both male sterile lines. These GO terms enriched suggest the general response or basic roles in the response to T. horrida. However, transcription factor activity was highly enriched in Jiangcheng 3A, but not in 9311A (Fig. 3b, d). For DEGs in Pro le 0, GO enrichment analysis showed that the genes enriched in pro le 0 were related to transporter activity and signal transducer activity in the two rice cultivars (Fig. 3b, d). These two GO terms enrichment showed the down-regulated genes might cut off the in ow of nutrition to T. horrida [36].

Analysis of metabolic pathways after T. horrida inoculation
Resistance to T. horrida is the result of a common effect of multiple genes, which play roles in certain biological processes by interacting with each other [37]. Pathway enrichment analysis of DEGs in both rice male sterile lines based on the KEGG database was performed, and the signi cantly enriched pathways are listed in Supplementary Table S7 and S8. Among these, 17 pathways, include photosynthesis-antenna proteins, phenylpropanoid biosynthesis, biosynthesis of secondary metabolites, phenylalanine metabolism, cysteine and methionine metabolism, terpenoid backbone biosynthesis, carotenoid biosynthesis, photosynthesis, plant-pathogen interaction, tyrosine metabolism, plant hormone signal transduction, pentose and glucuronate interconversions, diterpenoid biosynthesis, starch and sucrose metabolism, ribosome, carbon xation in photosynthetic organisms, and galactose metabolism were signi cantly enriched by T. horrida inoculation in both rice male sterile lines; However, pathways speci cally enriched in resistance male sterile line Jiangcheng 3A by T. horrida inoculation include cutin, suberine and wax biosynthesis, avonoid biosynthesis, glutathione metabolism, tyrosine metabolism, and biosynthesis of unsaturated fatty acids (Additional le 8 and 9: Table S7 and S8).

DEGs involve in plant-pathogen interaction
Plants employ two types of innate immune responses against invading pathogens: pathogen-associated molecule patterns (PAMPs) rst layer of immunity and effector-triggered immunity (ETI) second layer of immunity [38]. ETI is based on the recognition of pathogen-secreted effectors by plant-speci c resistance proteins [39]. We further analysis the DEGs involve in plant-pathogen interaction based on metabolic pathways analysis. In rice plants inoculated with T. horrida, 106 and 75 DEGs were found to be associated with plant-pathogen interaction in Jiangcheng 3A and 9311A, respectively (Additional le 10: Table S9).
Among them, 67 were common to both rice cultivars, and 39 were unique to Jiangcheng 3A (Supplementary Table S9

DEGs involved in plant hormone signal transduction
Once fungal effectors overcome the plant defence mechanism, plant resistance response becomes ineffective. In addition to the production of early signal molecules induced by pathogen infection, elicitor signals are frequently increased via the secretion of secondary signal molecules, such as ethylene (ET), salicylic acid (SA) and jasmonic acid (JA) [40]. Therefore, the DEGs related to plant hormone signal transduction were chosen for further analysis. A total of 92 and 81 DEGs were found to be associated with plant hormone signal transduction in Jiangcheng 3A and 9311A, respectively (Additional le 11: Table   S10), and more than half of all these genes were common to both rice cultivars, there are seven DEGs were involved in the JA signaling pathway, including six jasmonatezim-domain (JAZ) genes and one jasmonyl-L-isoleucine synthase gene (OsJAR2); two DEGs related to the SA signaling pathway, including a key regulation factor of SA mediated disease resistance NPR1-like gene OsNPR1 and TGA gene OsNIF1; three DEGs related to the ET signaling pathway, including one ET receptor gene, two ET-insensitive genes.
In addition, abscisic acid (ABA) signaling pathway-related genes, such as protein phosphatase 2C (PP2C) and stress-activated protein Kinase genes, pyrabactin resistance-like abscisic acid receptor genes, and bZIP transcription factor; auxin (AUX) signaling pathway-related genes, including small auxin upregulated RNA genes, indole-3-acetic acid-amido synthetase gene, and Aux/IAA genes were found in our results.
Furthermore, we also found 39 DGEs associated with plant hormone signal transduction were unique to Jiangcheng 3A (Fig. 4b). Among these genes, one gibberellin (GA) signaling pathway-related gene OsGAI and three cytokinin (CTK) signaling pathway-related genes OsRR6, OsRR4, and OsRR1 were found. This result showed that GA and CTK signaling pathway may play important role in regulated of T. horrida resistant.

Transcriptional factors (TFs) identi ed from DEGs
TFs play a key role in plant defence responses [41]. A total of 760 and 647 TFs were identi ed from the DEGs using the plant TF database PlantTFDB version 4.0 (http://planttfdb.cbi.pku.edu.cn/) in Jiangcheng 3A and 9311A, respectively (Fig. 5a). Interestingly, compared with the control conditions, 89 up-regulated and 74 down-regulated TF genes were identi ed in inoculated Jiangcheng 3A kernels, while 8 up-regulated and 90 down-regulated TF genes were identi ed in inoculated 9311A at 8 h (Fig. 5b). At 12 h, these numbers changed to 78 up-regulated and 71 down-regulated TF genes in Jiangcheng 3A, and 122 up-regulated, 142 down-regulated TF genes in 9311A (Fig. 5c). Although most of the TF genes exhibiting altered expression after T. horrida infection were common to both male sterile lines, some were unique to Jiangcheng 3A.

Discussion
RNA sequencing and DEGs in response to T. horrida RKS caused by the fungus T. horrida, is the most the most widespread disease occurred in many eld of hybrid seed produce, caused large yield lose. However, little information is known about the resistance pathways for T. horrida. As our previously reported, the rice male sterile line Jiangcheng 3A was highly resistant to T. horrida, and displayed consistent results between inoculation experiments carried out in 3 years [15]. We further chose two male sterile lines, Jiangcheng 3A and 9311A, as resistance and susceptible varieties, and performed comparative transcriptome analysis of different inoculation time points. Main aim was identify key pathway and genes involved in the resistant for T. horrida. Further, comparative analysis of DEGs of Jiangcheng 3A and 9311A revealed 4, 729 DEGs were commonly involved both rice male sterile lines, when we analysis the expression pro le of these common 4, 729 DEGs of both the genotypes, numbers of DEGs in sustained up regulation and sustained down regulation was more compared to others expression pro le. down-regulated, account for 42.58% and 57.42% of all DEGs, respectively. This result showed that resistance male sterile line Jiangcheng 3A has the higher percentage of genes up-regulated during the course of T. horrida inoculation, and suggests that these up-regulated genes may confer resistance to Jiangcheng 3A during T. horrida inoculation. Our results were different from Teqing and Lemont, in which the number of up-regulated DEGs were more in Lemont (susceptible rice cultivar) than in TeQing (moderately resistant rice cultivar) after Rhizoctonia solani inoculation, probably due to difference of rice response inoculation of necrotrophic (R. solani) and biotrophic (T. horrida) fungus [22]. Interestingly, upon exposure to low temperature stress, in cold tolerant variety were present higher percentage of genes upregulated than in cold susceptible variety [42]. In addition, at early (8 h) and late inoculation stages (72 h), more number of genes were participated the pathway of resistant to T. horrida in Jiangcheng 3A compared with 9311A, indicated that primary pathogen colonization may be occurred in 8 h, and the pathogen then established re-infection, causing the degree of infection to be enhanced at 72 h.

Functional characterization of DEGs
For GO enrichment analysis, some common terms were enriched in both male sterile lines; however, there are some terms, such as nuclear chromatin, protein-DNA complex, nucleosome, nuclear chromosome part, cytosolic ribosome, ribosome, ribosomal subunit, cytosolic part, non-membrane-bounded organelle, intracellular non-membrane-bounded organelle, ribonucleoprotein complex, carbohydrate binding, metal ion binding, urea transmembrane transporter activity, structural constituent of ribosome, and structural molecule activity molecular function, were unique enriched in resistant male sterile line Jiangcheng 3A.
We speculated that genotype that lacks these functions may be susceptible to T. horrida.
At 8 and 12 h as early infection times, a very high number of DEGs were enriched in hydrolase activity, hydrolyzing O-glycosyl compounds, hydrolase activity, acting on glycosyl bonds, cation homeostasis, photosynthetic electron transport chain, response to radiation, and plastoglobule in Jiangcheng 3A as compared to 9311A. During 24-72h of T. horrida inoculation as late infection times, similar to early response DEGs, more DEGs were enriched under T. horrida infection. In this response phase, there are 86 terms enrichment in Jiangcheng 3A; however, 9311A was enriched in 60 terms, such as carbohydrate metabolic process, ion transport, cation transport, urea transmembrane transporter activity, peptide biosynthetic process, amide biosynthetic process, and regulation of catalytic activity were unique enriched in Jiangcheng 3A; indicated that there are more GO terms involve in resistant variety and these terms may play an important role for ability of Jiangcheng 3A resistant to T. horrida.
Key pathway resistance to T. horrida in rice KEGG analysis of DEGs showed that cutin, suberine and wax biosynthesis, avonoid biosynthesis, tyrosine metabolism, glutathione metabolism, and biosynthesis of unsaturated fatty acids were only enriched in Jiangcheng 3A. Flavonoid and tyrosine ammonialyases play important roles in plant defence early during infection [43]. In addition, glutathione and some unsaturated fatty acids were also participates in signaling of plant defense against biotic and abiotic stresses [44]. These results indicated that these four pathways may play important role for ability of Jiangcheng 3A resistance T. horrida.
A total of 6 DEGs were found to be associated with avonoid biosynthesis. Among them, including two caffeoyl-CoA O-methyltransferase genes, two cytochrome P450 genes and one leucoanthocyanidin reductase genes, and one transferase family protein (Table 1). Former reports showed that caffeoyl-CoA O-methyltransferase gene may involve in regulation of lignin and other metabolites of the phenylpropanoid pathway and can enhance the multiple pathogens resistance of host plants [45]. In our results, two caffeoyl-CoA O-methyltransferase genes, Os08g0498100 and Os09g0481400 were sustained up-regulation in Jiangcheng 3A kernels after T. horrida inoculated; however in 9311A, only Os08g0498100 was up-regulation at 24 h after T. horrida inoculated (Table 1). Therefore, we speculate that these two genes were involved in pathway of T. horrida resistance. Proanthocyanidins are the polyphenolic compounds synthesized via the avonoid biosynthetic pathway and play a key role in defence against pathogens [46]. Leucoanthocyanidin reductase is a key enzyme of biosynthesis of proanthocyanidins, can converting leucocyanidin to catechin [46]. Leucoanthocyanidin reductase gene Os04g0630800, was sustained up-regulation in both male sterile lines after T. horrida inoculated; however, the gene expression levels (log 2 FPKM) was higher in Jiangcheng 3A than in 9311A at 8 h (Table 1).
Plant cutin and wax are important extracellular hydrophobic substances that prevent the entry of pathogens into host cells [47]. In this study, two cutin, suberine and wax biosynthesis-related genes ONI3 and DPW were identi ed with up-regulation in the Jiangcheng 3A (Table 1). Expression levels of ONI3 and DPW increased following T. horrida infection were happen in 8, suggesting that cutin and wax are required for protecting rice plants from the invading pathogens during the early stages of infection. In sum, we carefully concluded that Jiangcheng 3A kernels are capable of preventing T. horrida infection enhancing cell wall keratinization and by synthesising more secondary metabolites early during infection.
Moreover, there are six biosynthesis of unsaturated fatty acids-related genes were up-regulation in Jiangcheng 3A, including one oxidoreductase, two acyl-desaturase, one plastid ω-3 desaturase, and two microsomal Δ12-fatty acid desaturase genes (Table 1). Furthermore, four up-regulation genes were involved in the tyrosine metabolism, including one nicotianamine aminotransferase gene, one tyrosine decarboxylase gene, and two dehydrogenase genes. In addition, there were ten up-regulation genes related to glutathione metabolism, including one spermidine synthase gene, one RNRS1 homolog gene, one large subunit of ribonucleotide reductase gene and one glucose-6-phosphate 1-dehydrogenase gene, one dehydrogenase gene, one L-ascorbate peroxidase 4 gene, one glutathione reductase gene, and three glutathione S-transferase genes. Among these genes, eight were present signi cantly different expression between the both male sterile lines. Three Glutathione S-transferase (GSTs) encoded gene Os09g0467200, Os01g0369700, and Os01g0949800 were signi cantly up-regulated in Jiangcheng 3A, while not signi cantly different expression in 9311A. GSTs are ubiquitous and multifunctional enzymes in plant, and play a crucial role in plant defense though conjugation with glutathione, the attenuation of oxidative stress and the participation in hormone transport [48]. Notably, multiple GST genes were shown to be strongly inducible by the key defense hormone salicylic acid, which plays a key role in the activation of defence responses against biotrophic fungi [48]. Oxidoreductase and dehydrogenase as known to participate in generation of reactive oxygen species, which is one of several defence mechanisms used by plants against invading pathogens [49]. In our study, glucose-6-phosphate 1-dehydrogenase gene Os02g0600400 was up-regulated in both male sterile lines; however, the gene expression levels (log 2 FPKM) were higher in Jiangcheng 3A than in 9311A after T. horrida inoculation ( Table 1). The roles of these genes in resistant T. horrida need to be further examined.

DEGs involve in plant-pathogen interaction
We identi cation 39 plant-pathogen interaction-related DEGs were exclusively present in resistance male sterile line Jiangcheng 3A. Among these, 10 DEGs were up-regulated and 29 DEGs were down-regulated after T. horrida inoculation (Fig. 4a); 10 up-regulation DEGs including OsWRKY24, OsWRKY53, OsPBL1, NB-ARC/LRR disease resistance protein gene Os04g0514600, OsCML14, OsCML7, OsJAZ8, hsp82B, SPK, and Osrboh9. Intracellular calcium transients participate in some biological processes associated with resistance, such as oxidative burst, PR gene expression, hypersensitive cell death, and systemic acquired resistance during plant-pathogen interactions [50,51,52]. In this study, OsCML7 and OsCML14 were upregulated in Jiangcheng 3A kernels, we speculate that cellular calcium changes may participate in regulation pathway of T. horrida resistance, but the underlying mechanism should be further studied. Moreover, the highly conserved Heat shock proteins are involved in resistance to elevated temperatures and other environmental stresses [53], for example HSP90 interacts with a SA-induced protein kinase, which acts as a key component in signalling during plant defence [54]. Heat shock proteins hsp82B was induced in Jiangcheng 3A after infection with T. horrida, showed that it may also play important roles in rice tolerance to T. horrida, but its functions should be further studied. Osrboh9 encoding respiratory burst oxidase homologue protein, and participate in generation of reactive oxygen species [55]; OsPBL1 encoding serine/threonine protein kinase, and play important role in SA-mediated defense [56]; these two genes were up-regulation at 8 h and 72 h in Jiangcheng 3A, respectively. This result showed that ROS may regulation T. horrida resistant at early stage and SA play a role at late stage. In addition, two WRKY transcription factor genes association with rice resistant OsWRKY24 and OsWRKY53 were up-regulation at 8 h, their function are required for Jiangcheng 3A' ability to resistance T. horrida. For 29 down-regulation DEGs, several JAZ proteins encoding genes, such as Os04g0653000 and Os09g0401300 were downregulated at 8 h. Jasmonic acid (JA) is an important plant hormone involved in regulation of plant defence, and JAZ proteins function as transcriptional repressors during in JA signaling [57]. Therefore, we speculate that these genes negatively regulate resistance to T. horrida.

DEGs involved in plant hormone signal transduction
Plant hormones are involved in many biological processes in enhanced resistance to environmental stresses, diseases and pathogen infections [58]. In the present study, many DEGs were related to plant hormone (including SA, JA, ET, and IAA) signal transduction. Several lines of evidence show that SA plays a key role in the activation of defence responses against biotrophic fungi [59]. Notably, the key gene OsNPR1 of SA-mediated defense pathway was up-regulation in both male sterile lines after T. horrida inoculation [60]; however, the gene expression level was higher in Jiangcheng 3A than in 9311A. TGA gene OsNIF1 could negatively regulates resistance to Xanthomonas oryzae pv. oryzae in rice [61]; in our result, this gene present down-regulation in both male sterile lines after T. horrida inoculation. We speculate that these 2 SA signal-related genes may play key roles in T. horrida resistance, but the speci c functions of these genes should be further studied. ET and JA were also involved in resistance to pathogen infections [62]. In this study, several ET and JA signal-related genes were considered to be related to T. horrida resistance in rice. Some ethylene receptor, such as ETR2 were found to be involved in oral stage, overexpression of ETR2 in transgenic rice plants reduced ethylene sensitivity and delayed oral period [63]. ERS2 was up-regulated in Jiangcheng 3A after T. horrida infected but were not signi cantly changed in 9311A. We speculate that ERS2 play a role in T. horrida resistant though in uencing owering. The jasmonyl-L-isoleucine synthase gene OsJAR2 was induced by T. horrida inoculation in both male sterile lines, the expression levels of OsJAR1 and OsJAR2 were induced in response to wounding and the expression of OsJAR1 was up-regulation after blast infection [64]. Jasmonatezim-domain (JAZ) proteins function as transcriptional repressors during in JA signaling, several JAZ proteins encoding genes, such as Os03g0180900 and Os10g0392400 were down-regulated in both male sterile lines.
Furthermore, ABA responsive element binding factor OsABF1 was a suppressor of oral stage, inhibit its expression in transgenic rice plants was delayed oral period [65]. In this study, OsABF1 was downregulation in Jiangcheng 3A after T. horrida inoculation; however, not signi cantly changed in 9311A, so it may play an important role in T. horrida resistant though in uencing owering. Cytokinins (CTK) are plant hormones that play a key role in plant morphology, plant defence, leaf senescence and source-sink relationships [66]. OsRR1 was also involve in regulation of owering time, OsRR1-overexpressing plants show a late-owering phenotype [67]. In the present study, OsRR1 was up-regulation at 12 h in Jiangcheng 3A, but not signi cantly changed in 9311A. Therefore, we speculate that this gene may play similar roles as ERS2 and OsABF1, could hinder the T. horrida infection by delaying owering.

Transcriptional factors (TFs) identi ed from DEGs
More than half of the 103 WRKY TF genes present in the rice genome are differentially expressed under biotic and abiotic stress conditions [68,69]. In our study, 34 WRKY genes showed altered expression in Jiangcheng 3A kernels following T. horrida infection, of which 23 were up-regulated and 11 were downregulated at different time points. By contrast, 16 WRKY genes were showed altered expression in 9311A kernels following T. horrida infection, of which 10 were up-regulated and 6 were down-regulated (Additional le 13: Table S11). Out of the total DEGs WRKY genes detected, 13 were commonly regulated while 3 and 21 were exclusively regulated in 9311A and Jiangcheng 3A, respectively. Among 21 WRKY genes were unique to Jiangcheng 3A, OsWRKY12 and OsWRKY70 were up-regulated in Jiangcheng 3A at 8h post inoculation; suggesting that these gene play a positive role in resistance to T. horrida (Fig. 4c). Additionally, it has been shown that WRKY70 induces SA biosynthesis and inhibits JA biosynthesis. The expression level of WRKY70 is an indicator of the balance between SA and JA and determines which one of these hormones is more important [70]. OsWRKY12 is located upstream of OsNPR1 as a transcriptional activator in SA-dependent or JA-dependent defense signaling cascades [71]. OsWRKY10 was downregulation at 48 and 72 h, and this gene was increased by JA treatment. These result showed that SA may play a more important role than JA in defence against T. horrida, Balbi and Devoto [72] showed that SA is more important than JA in defence against biotrophic fungi, which is consistent with our results. In addition, OsWRKY22, OsWRKY24, OsWRKY15, OsWRKY23, OsWRKY44, OsWRKY93, OsWRKY69, OsWRKY50, and OsWRKY35 were also up-regulation in Jiangcheng 3A after T. horrida inoculation, but were not signi cantly changed in 9311A. The roles of these WRKY genes in resistant T. horrida need to be further identi ed.
MYB TFs are positive and negative regulators of pathogen resistance; for example, the Arabidopsis TF MYB30 promotes plant defence by positive regulation of HR [73]. In this study, we further found 60 MYB TF genes that were induced in Jiangcheng 3A following T. horrida infection, of which 25 were upregulated, and 35 were down-regulated at different times points. In 9311A, 50 MYB TF genes showing differential expressed were identi ed after infection with T. horrida, of which 21 were up-regulated and 29 were down-regulated at different infection time points (Additional le 14: Table S12). Out of the total DEGs MYB genes detected, 36 were commonly regulated while 14 and 24 were exclusively regulated in 9311A and Jiangcheng 3A, respectively. Among 24 MYB genes were unique to Jiangcheng 3A (Fig. 4d), Osmyb4 was up-regulated in Jiangcheng 3A at 8 and 12 h post inoculation, but were not signi cantly changed in 9311A after T. horrida inoculation. In Arabidopsis and cotton, overexpression Osmyb4 could increases biotic and abiotic stress tolerance [74]. We speculate that this gene play similar roles in T. horrida resistance regulation, but the speci c functions of this gene under T. horrida inoculation in rice should be further studied. MYB genes Os05G0543600 and Os11G0207600 were also up-regulated in Jiangcheng 3A at 8 h, Os04g0594100 at 24 h, Oa08g0433400 and Os11g0207600at 48 h, and Os01g0127450, Os01g0228901, Os01g0274800, Os01g0977300, Os04g0533250, Os05g0459000, OS06G0605750 at 72 h were up-regulated. Our results indicate that these MYB TFs may also play important roles in rice resistance to T. horrida; however, their functions should be further studied.

Conclusions
Transcriptome analysis showed that a resistant rice male sterile line exhibited more active responses to susceptible. DEGs in both male sterile lines were identi ed and enriched in many GO terms and pathways, suggesting that T. horrida resistance of rice is a complex mechanism. Glutathione metabolism-related genes, biosynthesis of unsaturated fatty acids-related genes, cutin, suberine and wax biosynthesis-related genes, avonoid biosynthesis-related genes, plant hormone signal transduction-related genes, plantpathogen interaction-related genes, and many TFs were found to be involved in T. horrida resistant in rice. The ndings reported herein increase our understanding of the molecular characteristics of rice resistant of T. horrida. The identi ed male sterile line Jiangcheng 3A will be useful for inbreeding new rice varieties with a high T. horrida resistance.