Plant endogenous xylanase inhibitors (XIs), discovered in 1997 by DEBYSER et al.^[1] in wheat, are plant-produced proteinaceous inhibitors that inhibit xylanases. To date, three structurally different classes of XIs have been identified in plants, Triticum aestivum xylanase inhibitor (TAXI)-type^[2], xylanase-inhibiting protein (XIP)-type^[3] and thaumatin-like xylanase inhibitor (TLXI)-type^[4]. TAXI, XIP and TLXI genes have already been characterized in wheat, rice, maize, barley, sorghum and rye^[5].

XIs are only active against xylanases of microbial origin and cannot inhibit the plant endogenous xylanases; therefore, it has been hypothesized that XIs function only in plant defense rather than in plant growth and development^[6]. Additional features that support possibly involvement of XIs in plant defense include the findings that XIs share sequence and structure homology with some pathogenesis-related (PR) proteins, and their expression is induced under stress conditions, and they inhibit the growth and germination of fungi. For example, XIP and TLXI can be classified as PR proteins of classes PR-5 and PR-8, respectively, based on their homology to thaumatin-like proteins and chitinases. Some members of TAXI and XIP gene families are significantly induced by wounding, pathogen infection and treatment with jasmonic acid (JA) and methyl jasmonate^[7-8]. A chitinase-like XIP from coffee, Coffea arabica, inhibited the germination of spores of soybean Asian rust Phakopsora pachyrhizi^[9]; and XIP and TAXI inhibited the growth of Rhizoctonia solani and Fusarium graminearum, respectively^[10].

Rice was reported to contain only the XIP type XIs, including RIXI, OsXIP, riceXIP and OsHI-XIP^{[8, 11-13]}. The RIXI gene (LOC_Os11g47580) was first cloned from rice in 2005, and it inhibited fungal GH11 xylanases, but did not inhibit GH11 xylanases from bacteria and GH10 xylanases from both fungi and bacteria^[11]; and RIXI was upregulated in rice seedlings infected by Magnaporthe grisea^[14].

Although most studies agree that XIs are involved in plant defense, little is known about additional functions of XIs in rice growth, XI-related genes and their expression pattern.

The RNA-seq method generates absolute information, and is more sensitive for transcripts with low expression. In this study, to gain more insight into the function of XIs in plant development, deep RNA sequencing was employed in combination with digital gene expression (DGE) profiles to analyze the difference between the transcriptomes of R7 (RIXI-overexpression transgenic line 7) and WT (wild-type). Additionally, the biological functions of RIXI, XI, and the effects of overexpressing RIXI on expressions of other rice genes as well as the mechanism of its role in plant defense were illustrated.

1 Materials and methods 1.1 Plant materials and growth conditions

Rice (Oryza sativa L. ssp. japonica cv. Nipponbare) was used in this study. R7 plants were RIXI overexpressing transgenic rice plants with single copy integration^[15], which were generated by placing RIXI under control of the strong cauliflower mosaic virus (CaMV) 35S promoter in our laboratory. WT and R7 were grown in normal culture solution at pH 5.0-5.5 in a glasshouse at 28 °C with 16 h of light and 8 h of darkness. The leaves of 2-week-old plants were collected for sequencing. The experiments were repeated three times. The experimental materials for sequencing were obtained by mixing 50 leaves together in each group.

1.2 Library construction and sequencing

At least 3 μg of total RNA per sample was sent to the Beijing Genomics Institute for Solexa sequencing. Sequencing libraries were generated using Illumina TruSeq RNA sample preparation kit (Illumina, San Diego, CA, USA) following the manufacturer’s recommendations. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE cluster kit v3-cBot-HS (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2000 platform and 100 bp single-end reads were generated.

1.3 Data analysis

Raw data (raw reads) of FASTQ format were first processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapters, those containing poly-N and low-quality reads from raw data. At the same time, Q20, Q30, GC content and the clean data were calculated. All downstream analyses were based on the clean high-quality data.

Reference genome and gene model annotation files were directly downloaded from the genome website (ftp://ftp.ensemblgenomes.Org/pub/release-17/plants/fasta/oryza_sativa), and index of the reference genome was built using Bowtie v2.0.6 , and single-end clean reads were aligned to the reference genome using TopHat v2.0.9. HTSeq v0.5.4p3 was used to count the read numbers mapped to each gene. Then RPKM (reads per kilobase of exon model per million mapped reads) of each gene was calculated based on the length of the gene and read counts mapped to this gene^[16]. To obtain statistical confirmation of the differences in gene expression between the WT and R7 plants, we used the RPKM to normalize the expression level of genes. All the uniquely mapped reads were used for calculating RPKM values. The differentially expressed genes between the WT and R7 libraries were identified using DEG-Seq software^[17], and P-values were corrected using the method of BENJAMINI et al.^[18]. Genes with an adjusted P-value ＜ 0.05 found by DEG-Seq were assigned as differentially expressed.

1.4 GO and KEGG enrichment analysis of differentially expressed genes

Gene ontology (GO) enrichment analysis of differentially expressed genes was implemented by the GOseq R package, in which gene length bias was corrected. GO terms with corrected P-value ＜ 0.05 were considered significantly enriched by differentially expressed genes.

Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, organism and ecosystem from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS software to test the statistical enrichment of differentially expressed genes in KEGG pathways.

1.5 Quantitative real-time PCR analysis

For quantitative real-time polymerase chain reaction (qPCR), total RNAs from seedlings of different genotypes were extracted using RNA prep pure plant kit (TianGen Biotech Co. Ltd., Qiagen, http://www.qiagen.com), following the manufacturer’s protocols. Total RNAs were converted into cDNAs using PrimeScript RT reagent kit (TaKaRa, now Clontech, http://www.clontech.com/takara). Real-time PCR analysis used SYBR Premix Ex TaqTM Ⅱ (TaKaRa) on a Roche LightCycler480 real-time PCR system, following the manufacturer’s instructions (Roche, http://www.roche.com). Transcript levels of each mRNA were determined and normalized with the level of actin mRNAs using the ΔC_t method^[19]. Gene-specific primers are listed in Table 1.

1.6 Measurement of agronomic traits

Ten plants from each genotype were randomly selected and evaluated for five growth and development related agronomic traits^[20]. Plant height and panicle length were measured as the average of the ten plants, and number of tiller per plant and panicles per plant were the average which evaluated for all tillers and panicles on the ten plants (panicles with less than ten seeds were not counted)^[21-22]. Thousand-grain mass (TGM) was the average of a sample of 1 000 filled grains per plant. All data including at least three measurements were analyzed using SPSS 16.0 (SPSS Inc., Chicago, IL, USA) for analysis of variance (ANOVA) and Duncan’s test. P＜0.05 was considered statistically significant.

2 Results 2.1 Comparison of gene expression level between the two libraries

To investigate the expression patterns of genes in WT and R7 plants, the RNA extracted from leaves was sequenced by Illumina Hiseq 2000 platform. Comparing the data between WT and R7 plants, each of two replicate experiments showed highly correlated expression values (R²＞0.92), indicating that the results from these experiments were reliable.

Compared with WT, a total of 1 296 significantly changed genes were detected in R7 plants, in which 391 genes were upregulated and 905 genes were downregulated. The results also revealed that a total of 96 genes were upregulated by at least threefolds, and 431 genes were downregulated by at least threefolds in R7 plants. Overexpression of RIXI gene may have an inhibitory effect on expression of some genes, because RIXI itself is an inhibitory protein, and overexpressing one gene in rice can influence the expression levels of some other genes.

2.2 Functional categorization of differentially expressed genes

Gene expression was compared between R7 and WT rice plants, and GO (http://geneontology.org) assignments were used to classify the functions of the differentially expressed genes. The GO annotation of these genes was presented in Fig. 1. These genes were mainly classified into the following categories: biological process and molecular function with 19 and 11 functional groups, respectively (Fig. 1A, B). Among the 30 groups, five groups had more differentially expressed genes, such as single-organism metabolic process (GO: 0044710) with 260 genes (62 upregulated genes and 198 downregulated genes) in the main category of biological process. Anion binding (GO: 0043168) with 214 genes, small molecule binding (GO: 0036094) with 212 genes in the main categories of molecular function. These results showed that RIXI may be involved in the molecular function and the regulation of biological process, and plays a vital role.

2.3 Pathway enrichment analysis of the differentially expressed genes

The biological pathways affected by RIXI expression were evaluated by mapping the detected genes to reference canonical pathways in the KEGG. Based on the expression pattern in WT and R7 plants, the differentially expressed genes were assigned to 98 KEGG pathways. The pathways with the greatest representation by unique genes were the metabolic pathways (KO: osa01100) with 87 members (6.1% of background genes in metabolic pathway), biosynthesis of secondary metabolites (KO: osa01110) with 62 members (8.5%), plant-pathogen interaction (KO: osa04626) with 25 members (26.0%) and plant hormone signal transduction (KO: osa04075) with 20 members (13.7%). These annotations provide a valuable resource for investigating the function of RIXI in rice growth and pathogen defense.

2.4 Overexpression of RIXI modulates plant hormone signaling

The plant hormone signal transduction pathway contained a total of 146 genes, with 20 significantly upregulated or downregulated genes in R7 plants compared to WT. It should be pointed that seven jasmonate-signaling genes in R7 plants were downregulated (Table 2). These genes (log-ratio range of -6.302 8 to -1.281 5) included rice jasmonate ZIM-domain (JAZ) proteins OsJAZ6, OsJAZ7, OsJAZ8, OsJAZ9, OsJAZ10, OsJAZ12 and OsJAZ13.

Compared with WT plants, two regulatory protein NPR1 genes (LOC_Os01g09800 and LOC_Os03g46440, log-ratio range of -0.995 1 and -1.661 6) in transgenic plants were downregulated (Table 2). Ethylene (ET) signaling-related genes, such as ET receptor-like protein and EIN3-binding F-box protein 1, were upregulated in R7 plants. Overexpression of RIXI also influenced abscisic acid (ABA) signaling, and enhanced the expression of one gibberellin (GA) 20 oxidase gene (LOC_Os01g66100); while a chitin-inducible gibberellin-responsive protein 1 gene (CIGR1, LOC_Os07g36170) was downregulated, and two indole acetic acid (IAA) biosynthesis genes, indole-3-acetic acid-amido synthetases OsGH3.8 and OsGH3.2, were highly downregulated (Table 2).

2.5 RIXI expression and transcription factors (TFs)

Overexpression of RIXI induced changes in the expression of 95 TFs belonging to 23 TF families (Table 3), including AP2/EREBP, AUX/IAA, BES1, bHLH, bZIP, C3H, HSF, MYB, NAC, TCP, Trihelix, WRKY and TIFY, all of which play important roles in plant development and immunity. Of the 95 identified TF genes, 30 were upregulated, including two ethylene-responsive factor transcriptional regulator genes (LOC_Os09g11460 and LOC_Os04g46220) belonging to the AP2/EREBP family, two Aux/IAA TFs, seven bHLH genes, three C3H genes and one WRKY gene (LOC_Os08g38990). Meantime, 65 TF genes were downregulated, including six MYB TFs, seven TIFY TFs and 13 WRKY TFs (Table 3).

2.6 RIXI was involved in plant defense responses

The KEGG pathway enrichment analysis also revealed that the “plant-pathogen interaction” pathway was also significantly changed in R7 plants, indicating that RIXI overexpression influenced the expression of many genes involved in plant defense. This included such signal transduction components as Ca²⁺ and G-protein signaling, protein phosphatase, protein kinases, and defense response proteins including cyclic nucleotide-gated ion channel 2, calcium-binding protein and defense proteins against fungi or bacteria (Table 4).

Overexpression of RIXI also influenced the expression of genes involved in reactive oxygen species (ROS) production and antioxidants (Table 4). Of the responsive ROS-production genes (including amine oxidase, germin-like protein, respiratory burst oxidase homolog protein F and NADH dehydrogenase), three were upregulated and six were downregulated. Of the responsive genes involved in antioxidants, four were upregulated and 13 were downregulated (Table 4).

2.7 Confirmation of tag-mapped genes by qPCR analyses

To evaluate the reproducibility of the profile, 10 genes were selected for qPCR assays, including PR genes, TF genes, JA and salicylic acid (SA) biosynthesis and signaling genes. Expression levels of these genes (Cht, OsPR1a, OsPR1b, PR4, AOS, LOX, PAL, NH1, Pid2 and Myb) had similar patterns with qPCR to those of Tag-seq analysis (Fig. 2), indicating the validity of the profile.

2.8 Characteristics of transgenic rice plants

The transgenic lines R7 appeared to grow and develop normally. Phenotype parameters were measured at the stage of plant maturity. There were no significant differences in plant height, panicle length, number of tillers and panicles compared with non-transgenic plants (WT) (Table 5).The transgenic lines R7 had slight higher values for TGM. These results indicate that the overexpression of RIXI has little influence during the development stages of transgenic rice plants and don’t reduce the production of rice in total; moreover, it might have a positive effect on the yield of rice to a certain extent.

3 Discussion

This study is the first to study transcriptomic changes of transgenic rice (R7) and WT plants using deep RNA sequencing combined with DEG profile analysis. Compared with WT, the transcript levels of 1 296 genes in overexpressing transgenic R7 plants were differentially expressed. The annotated differentially expressed genes could be classified into two categories: biological process and molecular function with 19 and 11 functional groups, respectively, among which single-organism metabolic process, anion binding, small molecule binding, small molecule binding and nucleotide binding have more different expression genes. These results showed that RIXI may be involved in the molecular function and regulation of biological process by controlling the small molecular binding which involved in the metabolic process. According to KEGG pathway enrichment analysis, most differentially expressed genes were involved in metabolism, biosynthesis of secondary metabolites, plant-pathogen interaction and plant hormone signal transduction. By altering metabolism, plants could make more energy available for synthesis of defense materials and elicitation of the defense response^[23]. Changes in metabolism might be a general phenomenon in plant responses to stresses.

As it has been hypothesized that RIXI plays an important role in plant defense^[14], we focused on defense response genes, such as plant hormone signaling, TFs, signal transduction components and defense response protein genes.

Phytohormones are chemical messengers that coordinate cellular activities. The phytohormones of plants such as Auxin/IAAs, CK, GA, ABA, ET, BR and JA are chemical messengers that play distinct but overlapping roles in the regulation of defense response, growth and development of rice plants. In the present study, compared with WT, plants overexpressing RIXI downregulated JAZ proteins. JAZ transcriptional proteins acted as transcriptional repressors of JA responses in Arabidopsis and rice^[24-25]. This confirmed observations that XIs respond to stress via a JA-mediated signaling pathway^[8]. One of the major roles of ABA is to mediate adaptive responses to various environmental stresses during the plant growth. It has a critical role in the plant adaptation or acclimation to abiotic stresses, such as high salinity, drought, low temperature and mechanical wounding^[26]. The overexpression of RIXI also influenced ABA signaling, resulting in upregulation of two genes (LOC_Os02g15640 and LOC_Os01g59760) and downregulation of three genes (LOC_Os07g42940, LOC_Os03g16170 and LOC_Os02g52780). The LOC_Os01g59760 and LOC_Os02g52780 encode OsbZIP09 and OsbZIP23, respectively, which are members of the basic leucine zipper (bZIP) TF family in rice and regulate the expression of a wide spectrum of stress-related genes in response to abiotic stresses through an ABA-dependent regulation pathway^[26-27]. In addition, we found that the overexpression of RIXI elevated the mRNA levels of many genes involved in ET biosynthesis and signal transduction. ET plays a highly pleiotropic role in plant growth and development, and participates in responses to various biotic and abiotic stresses. Most genes involved in JA, SA, ABA, and ET biosynthesis or signaling pathways were altered, indicating that RIXI-induced plant defense responses are mainly dependent on the JA, SA, ABA and ET pathways.

All major processes of life depend on differential gene expression, which is largely controlled by TFs’ activity. Many TFs, such as MYB, NAC, WRKY, and TIFY family TFs, were downregulated or upregulated in R7 plants. Few genes of TFs, for example WRKY gene (OsWRKY30, LOC_Os08g38990) was highly increased in R7 plants. These downregulated TFs serve as negative regulators of defense signaling, for instance, OsWRKY71 has features characteristic of a transcriptional repressor of GA signaling in aleurone cells^[28], and OsWRKY45-1 negatively regulates ABA signaling and, in addition, OsWRKY45-2 negatively regulates rice response to salt stress^[29]. TIFY 10C/OsJAZ8 negatively regulated the JA-induced resistance to pathogens in rice and acts as a repressor of JA signaling in rice^[30].

Many genes involved in the biosynthesis pathways of secondary signaling compounds were regulated, including Ca²⁺ signaling, ROS, G-protein, phosphatidylinositol signaling, and protein kinase. Compared with transcript levels in WT plants, the expression levels of many genes encoding protein kinases, such as receptor-like protein kinase, mitogen-activated protein kinase (MAPK) and calmodulin dependent protein kinase, were changed in RIXI-overexpressing rice plants. Protein kinases represent an important mechanism in rice defense signal transduction, and MAPK signaling cascades are conserved in signaling pathways and have been implicated in a wide variety of plant biotic and abiotic stress responses^[31]. This indicates that the constitutive expression of RIXI in rice plants activated a complex signal transduction network and enhanced rice resistance to stress.

In conclusion, RIXI, a rice XI, involved in biological process and molecular function, especially activated a complex signal transduction network in response to various biotic and abiotic stresses. Meanwhile, the overexpression of RIXI didn’t have a negative influence on the growth and development of rice plants. More work should be done to understand the deep biological functions and regulatory mechanisms of this transgenic plant in future.