Comparative Transcriptome Analysis of Two Rice Varieties in Response to Rice Stripe Virus and Small Brown Planthoppers during Early Interaction

Rice stripe, a virus disease, transmitted by a small brown planthopper (SBPH), has greatly reduced production of japonica rice in East Asia, especially in China. Although we have made great progress in mapping resistance genes, little is known about the mechanism of resistance. By de novo transcriptome assembling, we gained sufficient transcript data to analyze changes in gene expression of early interaction in response to SBPH and RSV infection in rice. Respectively 648 and 937 DEGs were detected from the disease-resistant (Liaonong 979) and the susceptible (Fengjin) varieties, most of which were up-regulated. We found 37 genes related to insect resistance, which mainly included genes for jasmonate-induced protein, TIFY protein, lipoxygenase, as well as trypsin inhibitor genes and transcription factor genes. In the interaction process between RSV and rice, 87 genes were thought to be related to RSV resistance; these primarily included 12 peroxidase biosynthesis genes, 12 LRR receptor-like protein kinase genes, 6 genes coding pathogenesis-related proteins, 4 glycine-rich cell wall structural protein genes, 2 xyloglucan hydrolase genes and a cellulose synthase. The results indicate that the rice-pathogen interaction happened both in disease-resistant and susceptible varieties, and some genes related to JA biosynthesis played key roles in the interaction between SBPHs and rice. When rice was infected by RSV a hypersensitive reaction (HR) in the disease-resistant variety was suppressed, which resulted from an increase in peroxidase expression and down-regulation of LRR receptor-like protein kinase and pathogenesis-related proteins, while, the changes of peroxidase biosynthesis, glycine-rich cell wall structural protein, cellulose synthase and xyloglucan endotransglucosylase/hydrolase could lead to the strengthening of physical barriers of rice, which may be an important resistance mechanism to RSV in rice.


Introduction
Rice stripe virus (RSV), transmitted by a small brown planthopper (SBPH, Laodelphax striatellus Fallen) [1], has caused a disastrous disease of rice in East Asia, particularly in China, Japan, Korea and North Korea [2]. When infected with RSV at the seeding stage, rice grows poorly and often develops folded and twisted leaves, with the central leaves yellowing and withering; growth may terminate and ultimately the plant will die. When infection occurs after tillering, upper leaves become discolored or yellow-green and the number of grains per ear is reduced significantly [3]. In China, rice stripe is increasing in severity, particularly in Jiangsu province, where approximately 0.6 M ha per year of rice were infected by RSV in the period 2000 to 2003, rising to 1 M ha in 2004. Rice yield is reduced by 30-50% in heavily infected fields, and in some of the most severely infected fields, no harvest is possible [4]. Rice stripe is still one of the most damaging virus threats to Japonica production in China. It is difficult to control this disease because the virus is transmitted in a circulative, transovarial propagative manner by viruliferous female planthoppers to their offspring. The most economical and effective way for controlling the disease is to plant resistant varieties.
The resistance gene to RSV that has been frequently utilized in breeding is Stvb i , which originated in the indica cultivar 'Modan'. Hayanoo et al constructed a physical map spanning 1.8-cM distance between flanking markers, consisting of 18 bacterial artificial chromosome (BAC) clones, around the Stvb-i locus on rice chromosome 11 [5]. Despite considerable progress in mapping RSV resistance genes [6,7], there have been only a few studies of the resistance mechanism [8,9] Understanding the responses of rice to RSV infection is important for developing strategies for disease control. Since rice responses to viral infection are complex and relate to many kinds of physiological processes, system-level transcriptomic studies are required to fully understand the responses. Recently, nextgeneration deep-sequencing techniques such as Illumina RNA-Seq and digital gene expression have provided new approaches for studying the transcriptome. RNA-Seq is a whole transcriptome sequencing method, and many transcriptome studies have greatly extended our knowledge of mechanisms of resistance to plant pathogens [10,11,12].
In this study, we analyzed the early response of two rice cultivars to infection by RSV and its carrier at the transcriptome level using next-generation deep-sequencing techniques. We investigated the alteration in gene expression between a disease-resistant cultivar and a susceptible cultivar before and after inoculation with RSV by co-culturing with Laodelphax striatellus for 48 h. Our study provides insight at the molecular level into the mechanism of development of rice stripe disease, which contributes to our understanding of the rice-RSV interaction.

Results
Illumina sequencing and aligning to the reference genome Our transcriptome sequence data included gene expression profiling following different treatments of two rice varieties. The RNA-Seq method generates absolute information rather than relative gene expression measurements; thus, it avoids many of the inherent limitations of microarray analysis. This method was used to analyze differences in gene expression between the diseaseresistant and the sensitive hosts. We sequenced four cDNA libraries, R1 (Liaonong 979 CK), R2 (Fengjin CK), R3 (Liaonong 979 inoculated treatment), R4 (Fengjin inoculated treatment), and generated 45,763,250 sequence reads, each of which was 50-101 bp in length, encompassing 7.91 Gb of sequence data ( Table 1). GC% of sequence data from the four libraries were all approximately 50%, and CycleQ20% were all 100%, which showed that the accuracy and quality of the sequencing data were sufficient for analysis. Each treatment was represented by approximately 10 million reads, the tag density sufficient for the quantitative analysis of gene expression. Of the total reads, 70-80% matched to genomic locations. Sequence data from the four libraries were combined, and 33420 unigenes were finally obtained with an average depth of .70x and 55-60% coverage. The length distribution of total unigenes had similar patterns among the four samples, suggesting there was no bias in the construction of the cDNA libraries ( Figure 1).
Difference of sequencing data between the reference genome (Nipponbare) and two rice cultivars As is shown in Table 2, many new genes, new exons,UTRs, rtIntroes, skippedexons and alternativesplices were detected from the four cDNA libraries (Table S1-6). Annotations of new genes in eight kinds of database were listed in Table S7.
We detected a total of 3918 new genes by RNA-Seq, annotated results of the different lengths of the new genes are shown in Table 2. 982 new genes were categorized into 25 functional groups with a COG classification. Among these COG categories, the cluster of 'general function prediction' occupied the highest number, followed by 'transcription', 'replication, recombination and repair' and 'signal transduction mechanisms'. It's worth mentioning that 48 and 45 new genes were categorized into the defense mechanisms functional group (V) and cell wall/membrane/envelope biogenesis group (M), respectively ( Figure 2). Analysis of these genes will contribute to our knowledge of the resistance mechanism to RSV and its transmission vector.

Differentially Expressed Genes(DEGs)among the four rice samples
A total of 3332 DEGs were detected by RNAseq. There were 1011 differentially expressed genes between Liaonong 979 CK (R1) and Fengjin CK (R2), of which 575 genes were up-regulated and 436 genes were down-regulated. The results show that there are wide discrepancies in gene expression in seedlings of the rice cultivars. After inoculation by SBPHs that carried RSV, gene expression of the two rice varieties changed significantly. 648 differentially expressed genes were detected in Liaonong 979, among which expression levels of 465 genes were up-regulated;  only 183 genes were down-regulated. In Fengjin, even though it is sensitive to RSV, 937 differentially expressed genes were detected, among which 78.7% were up-regulated. There were 736 differentially expressed genes between R3 and R4 ( Figure 3). Inoculating rice cultivars, Liaonong 979 and Fengjin, respectively resistant and sensitive to RSV, with SBPHs carried RSV, resulted in interaction with pathogens. In this process, the two cultivars shared 267 DEGs in response to piercing and sucking by SBPHs according to the Venn diagram ( Figure 4). In addition, the number of DEGs between R1_vs_R3 and R2_vs_R4 was 381, and there were 399 DEGs between R3_vs_R4 and R1_vs_R2 (Table 3). These genes were predicted to be relevant to interaction between rice and SBPHs or RSV. The data illustrate the difference between the resistant and sensitive cultivars, and deep analysis of these genes may shed light on the resistance mechanism.
Differential gene expression profiles in response to SBPH feeding After being inoculated by SBPHs, Liaonong 979 and Fengjin shared 267 DEGs, in which 221 genes could be annotated in the Swiss-Prot protein sequence database (Table S8). There were 27 genes that, when up-regulated were expressed in one cultivar and down regulated in the other, while the remaining genes were all either up-regulated or down-regulated in each of the two cultivars.
Analysis of these 221 DEGs by annotation in the Swiss-Prot database showed that 37 DEGs related to interaction between rice and pathogens were all up-regulated in these two varieties. Because Liaonong 979 and Fengjin were all moderately resistant to SBPHs, the 37 DEGs having similar expression patterns in these two cultivars were considered to be related to SBPHs response. As was determined from the Swiss-Prot protein sequence database, these up-regulated DEGs included 1 jasmonate-induced protein gene, 5 TIFY protein genes, 4 lipoxygenase genes, 4 trypsin inhibitor genes, 6 transcription factor genes, 3 CBL-interacting protein kinase genes, 3 protein phosphatase 2C genes, 1 (E)-betafarnesene synthase gene ( Figure 5).
Jasmonic acid (JA) is a signal molecule extensively involved in the regulation of plant growth and response to stress [13,14,15]. In this study, piercing and sucking by the SBPHs induced expression of defense genes, of which JA was one of the most important factors. As is shown in Figure 5, RPKM of a gene coding jasmonate-induced protein changed most significantly, increasing from 4 to 7841 in Fengjin, and 10 to 2770 in Liaonong 979. In addition, 9 DEGs coding lipoxygenase and TIFY protein were related to JA biosynthesis. Lipoxygenase is a key enzyme of the octadecanoid pathway (also called LOX pathway) [16,17], and the TIFY protein family is related to JAZ protein. Wounding by chewing insects resulted in the localized and systemic accumulation of high levels of proteinase inhibitor proteins in leaves of several plant families [18]. In this study, after inoculating with SBPHs, 5 genes coding proteinase inhibitors were all markedly upregulated, indicating that proteinase inhibitors play critical roles in resistance to SBPHs resistance in addition to JA.
Analysis of key genes participating in the early interaction between rice and RSV 381 DEGs were only detected in RSV resistant variety (Liaonong 979), of which 289 could be annotated in the Swiss- Prot database (Table S9). We found that 87 genes were related to plant defense responses. As is shown in Figure 6, several DEGs were significantly up-regulated; of these 12 were related to peroxidase biosynthesis, 4 to coding glycine-rich cell wall structural protein, and 1 gene coding cellulose synthase. Several other genes were markedly down-regulated, including 12 genes coding LRR receptor-like protein kinase and 6 genes coding pathogenesis-related protein, which were all related to disease resistance, signal transduction and the hypersensitive necrosis reaction (HR). Down-regulated expression of these genes showed that HR was suppressed in Liaonong 979 during the early interaction between rice and RSV. The changes in peroxidase biosynthesis, glycine-rich cell wall structural protein and cellulose synthase could result in lignification and strengthening of cell walls [19,20,21], which indicates that RSV resistance may be related to physical barriers of cells.
We found that expression patterns of 24 genes were different in Liaonong 979 than in Fengjin (Table S10). These genes included 4 pathogenesis-related protein genes, 2 xyloglucan hydrolase genes, 3 chlorophyll a-b binding protein genes, 1 proline-rich protein gene. The results showed that RSV infection may trigger a series of biochemical reactions, and resistance to RSV may be related to changes in photosynthesis, xyloglucan hydrolysis, and proline enrichment.

New exons related to disease resistance
Induction of expression of 152 new exons was detected in R3 (Table S11), in which 32 and 82 new exons could be annotated respectively in COG and the Swiss-Prot databases. COG annotatation revealed 7 exons related to cell wall/membrane/ envelope biogenesis. In addition, it is reported that an h-type thioredoxin functions in tobacco defense responses to two species of viruses [22], we detected a new exon encoding TPR repeatcontaining thioredoxin. Whether these new exons are related to resistance to RSV is unknown.

Comparison of DGE tag data with qRT-PCR expression patterns
In order to validate our DGE data, ten annotated unigenes having annotations were selected for qRT-PCR analysis (Figure 7).   Table 3. Accounting of the DEGs shown in Figure 4. The results showed that the qRT-PCR data of these genes were consistent with the DGE results. For example, both qRT-PCR and DGE analyses showed that genes encoding jasmonate-induced protein, Bowman-Birk type bran trypsin inhibitor, lipoxygenase 2.1 and protein TIFY were more highly expressed in two inoculated rice samples (R3, R4) than in the two non-inoculated samples (R1, R2). Likewise, expression patterns of two genes encoding peroxidase in the two rice varieties were confirmed to be similar by DGE analysis and qRT-PCR analysis. We further analyzed expression of three genes encoding LRR receptor-like serine/threonine-protein kinase and disease resistance protein and xyloglucan hydrolase protein respectively, by DGE analysis; these genes were significantly down-regulated in the resistant variety but were little changed in the susceptible variety, which is also in agreement with our qRT-PCR results.

Discussion
HR may not be the mechanism for resistance to RSV in Liaonong 979 Plants are often attacked by pathogens, such us bacteria, fungi, stramenopiles, viruses, protozoa and nematodes. To combat the pathogens, plants have evolved sophisticated and effective defense mechanisms. The hypersensitive response (HR) is a common initial response to attack by pathogens. This type of resistance often occurs in the incompatible interactions between plants and pathogens, and is characterized by rapid necrosis of infected and neighboring host cells [23,24]. Previous transcriptome analysis using Illumina sequencing-based DGE led to the conclusion that HR could be considered as the mechanism of resistance to TMV and CMV [25]. However, our data indicated that HR may not be the mechanism for resistance to RSV in Liaonong 979. As can be seen in Figure 6, which shows plant disease resistance signals in the disease-resistant cultivar, the expression of 22 genes, encoding receptor-like protein kinase and pathogenesis-related protein did not increase, rather these genes were down-regulated significantly after infesting rice with SBPHs that carried RSV; Thus HR was not triggered by RSV attack in Liaonong 979.
Lignification and cell wall strengthening may play a critical part in resistance to RSV In addition to HR mediated by pathogenesis-related proteins, lignification was considered to be an important mechanism of resistance to pathogens in plants [26]. Lignin and callose can function as physical barriers to prevent propagation and movement of pathogens [27]. In this study, after inoculation, some DEGs including 12 peroxidase genes,4 genes coding glycine-rich cell wall structural protein and a cellulose synthase gene were markedly up-regulated in the disease-resistant host (Liaonong 979) and down-regulated or unchanged in the susceptible cultivar Fengjin. Peroxidases, glycine-rich cell wall structure protein and cellulose synthase have a role in defense by strengthening cell physical barriers and lignification, which indicates that strengthening physical barriers may be one of the factors for RSV resistance in rice. Previous studies have used different methods to characterize changes in expression of some genes in response to RSV infection. By microarray, Satoh et al studied changes in gene expression at 5 stages, respectively 0 d, 3 d, 6 d, 9d and 12d after inoculation with SBPHs carrying RSV. At 12 days postinoculation, two families of genes, lignification and cellulose synthesis were all significantly down-regulated in Nipponbare (susceptible to RSV), which is consistent with the results in this test.

Signal pathway induced by SBPH feeding might be similar to BPH in rice
In this study, a few genes including the jasmonate protein gene, TIFY protein genes and lipoxygenase genes were up-regulated significantly after these two rice varieties were co-cultured with SBPHs for 48 h, which indicated that the JA signaling pathway was activated by SBPH piercing and sucking. Similar results were obtained from a study of gene expression during feeding by brown planthoppers (BPH) [28]. Using microarray, expression of some genes that take part in a jasmonic acid-independent pathway were detected to be up-regulated when rice plants were co-cultured with BPH for 72 h,. Plant defense responses to attack by insects include the activation of pathways dependent on SA and JA/ethylene signaling molecules [29,30,31]. Results of microarray and transcriptome sequencing showed that signal pathways induced by BPH and SBPH were all JA signaling.

Conclusion
By using next generation sequencing, abundant and high quality data were gained for the analysis of early interaction of rice in response to SBPH and RSV infection. The results indicated that expression of many genes associated with JA signaling, pathogenrelated molecular patterns and cell wall structure would be activated or suppressed before and after being inoculated with SBPH carried RSV. Further analysis shows that HR may not be the mechanism for resistance to RSV in a disease resistant variety, and lignification and cell wall strengthening may play a critical part in resistance to RSV. But the molecular function of some genes remains largely unknown and requires further study. Even so, the identification of DEGs involved in rice stripe disease resistance will extend our understanding of the complex molecular and cellular events in this process and provide a foundation for future studies on breeding for resistance to rice stripe disease.

Plant materials
Liaonong979 is a highly resistant rice variety to RSV, and Fengjin is a highly susceptible cultivar to rice stripe disease, both of these varieties are moderately resistant to SBPH [32].

Obtaining samples for transcriptomic analysis
Two rice cultivars, Fengjin and Liaonong979, were grown in 1000 ml beakers, 20 seeds per beaker, in a green house (2561uC, natural sunlight). When the most vigorous seedlings reached the three-leaf stage, 100 small brown planthopper (SBPH) nymphs in the second or third instars were scattered into each beaker, 90% of which were viruliferous, and the virus-carrying rates were calculated by Enzyme-linked Immuno sorbent Assay (ELISA). SBPHs were scattered around the beaker three times daily to prevent them from remaining in one area. After 48 h, the insects were completely removed from the plants and parts of the leaves of every seedling were mixed for RNA extraction. RNA was extracted by the same method from two uninoculated controls (the control experiment without insects feeding, R1, R3).

RNA extraction and quality determination
Total RNA was isolated using the modified CTAB method performed as described by Chang et al [33], and the RNA samples were treated with DNase I (TaKaRa, Japan) for 4 h. The integrity of the RNA samples was examined with an Agilent 2100 Bioanalyzer.
cDNA library preparation and Illumina sequencing cDNA library preparation and sequencing reactions were conducted by the Biomarker Technology Company, Beijing, China. Preparation of the paired-end libraries and sequencing were performed following standard Illumina methods and protocols. The cDNA library was sequenced on the Illumina Cluster Station and Illumina Genome Analyzer system.

De novo assembly
Reads from each library were assembled separately. The Trinity method [34] was used for de novo assembly of Illumina reads of hawthorn (Crataegus sp.). Trinity consists of three software modules: Inchworm, Chrysalis and Butterfly, applied sequentially to process large volumes of RNA-Seq reads. In the first step in Trinity, reads are assembled into the contigs by the Inchworm program. The minimally overlapping contigs were clustered into sets of connected components by the Chrysalis program, and then the transcripts were constructed by the Butterfly program. Finally, the transcripts were clustered by similarity of correct match length beyond the 80% of longer transcripts or 90% of shorter transcripts using the multiple sequence alignment tool BLAT [35]. The raw sequence data of four samples in this test have been uploaded to NCBI(http://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?), and the accession number is SRP028592.

Functional annotation
We annotated unigenes based on a set of sequential BLAST searches [36] designed to find the most descriptive annotation for each sequence. The assembled unigenes were compared with sequences in the National Center for Biotechnology Information (NCBI) non-redundant (Nr) protein and nucleotide (Nt) databases, the Swiss-Prot protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database, the Cluster of Orthologous Groups (COG) database, the Translated EMBL Nucleotide Sequence Database (TrEMBL) and InterPro database. The Blast2GO program [37] was used to obtain GO annotation of the unigenes. The WEGO software was then used to perform GO functional classification of all unigenes to view the distribution of gene functions.

Digital gene expression analysis
Gene expression levels were measured in RNA-Seq (Invitrogen) analyses as numbers of reads and were normalized with RPKM [38]. IDEG6 software [39] was used to identify differentially expressed genes in pair-wise comparison, and the results of all statistical tests were corrected for multiple testing with the Benjamini-Hochberg false discovery rate (FDR,0.01). Sequences were deemed to be significantly differentially expressed if the adjusted P value obtained by this method was ,0.001 and there was at least a two fold change (.1 or ,21 in log 2 ratio value) in RPKM between two libraries.

Quantitative RT-PCR (qRT-PCR) analysis
To confirm the results of pyrosequencing analysis, we determined the expression levels of 10 DEGs by qRT-PCR. Total RNAs from each sample were extracted using TRIzol reagent (Invitrogen) and treated with DNase I (Invitrogen) according to the manufacturer's protocol. The concentration of each RNA sample was adjusted to 1 mg/ml with nuclease-free water and total RNA was reverse transcribed in a 20 ml reaction system using the AMV RNA PCR Kit (TaKaRa). qPR-PCR was carried out on the LightCycler 480@ II with LightCycler 480@ SYBR I Master (Roche Applied Science, Basel, Switzerland) under the following conditions: 95 mC for 5 min; and 40 cycles of 95 mC for 10 s, 60 mC for 15 s, and 72 mC for 20 s, followed by melting curve generation (68uCto95uC). Primers used in qRT-PCR for validation of differentially expressed genes were shown in Table 4.

Supporting Information
Table S1 New genes detected from the four cDNA libraries. These genes are those which have not been found in a japonica cultivar Nipponbare (below). (XLS)       Table S10 24 DEGs whose expression patterns were adverse in Liaonong979 and Fengjin. These genes were considered to be related to RSV resistance. (XLS)