Transcriptome Analysis of Ramie (Boehmeria nivea L. Gaud.) in Response to Ramie Moth (Cocytodes coerulea Guenée) Infestation

The ramie moth Cocytodes coerulea Guenée (RM) is an economically important pest that seriously impairs the yield of ramie, an important natural fiber crop. The molecular mechanisms that underlie the ramie-pest interactions are unclear up to date. Therefore, a transcriptome profiling analysis would aid in understanding the ramie defense mechanisms against RM. In this study, we first constructed two cDNA libraries derived from RM-challenged (CH) and unchallenged (CK) ramie leaves. The subsequent sequencing of the CH and CK libraries yielded 40.2 and 62.8 million reads, respectively. Furthermore, de novo assembling of these reads generated 26,759 and 29,988 unigenes, respectively. An integrated assembly of data from these two libraries resulted in 46,533 unigenes, with an average length of 845 bp per unigene. Among these genes, 24,327 (52.28%) were functionally annotated by predicted protein function. A comparative analysis of the CK and CH transcriptome profiles revealed 1,980 differentially expressed genes (DEGs), of which 750 were upregulated and 1,230 were downregulated. A quantitative real-time PCR (qRT-PCR) analysis of 13 random selected genes confirmed the gene expression patterns that were determined by Illumina sequencing. Among the DEGs, the expression patterns of transcription factors, protease inhibitors, and antioxidant enzymes were studied. Overall, these results provide useful insights into the defense mechanism of ramie against RM.


Introduction
Ramie (Boehmeria nivea L. Gaud.), also called China grass, is a perennial herbaceous plant that belongs to the family Urticaceae. It is an important natural fiber crop, mainly grown in China, India, and other Southeast Asian and Pacific Rim countries [1]. Traditionally, ramie was planted as a fiber crop, and only the best fibers were harvested (constituting approximately 4% of the total dry matter). However, studies over the past decades have shown that ramie is rich in protein [2], and the shoots and leaves can be used as fodder for cattle and geese [3].
Cocytodes coerulea Guenée, also known as the ramie moth (RM), is a destructive plant pest. It can cause severe ramie yield reductions, as it feeds on ramie leaves and new shoots [4]. This pest is widely distributed in China, Japan, India, and Southeast Asian countries. In central and south China, up to four generations per year can occur during the vegetative growth phase of ramie [5]. Female RMs lay eggs (approximately 400 eggs each) on the abaxial leaf surface of ramie, which causes the leaves to yellow. After the eggs' hatching, the RM larvae begin to feed on the leaves, which will result in a net-like pattern of damage (Figure 1), and will seriously impair the plant's photosynthetic capacity, resulting in a yield reduction.
Chemicals such as dichlorvos, Trinox, and pyrethrum ester insecticides are effective in controlling RM larvae. However, the widespread use of agrochemicals can result in serious problems referred to as "3R-problem": residue, resistance, and resurgence. Selecting RM-resistant cultivars is therefore considered to be the most economical, effective, and environment-friendly control measure for reducing RM damage. A previous study showed that different ramie cultivars displayed varied levels of resistance to RM, which were correlated to the phenol, tannin, and oxadate contents in the leaves [6]. However, the underlying molecular mechanism of this resistance to RM is unclear. Methods such as suppression subtractive hybridization, cDNA-amplified fragment length polymorphism (cDNA-AFLP), and microarrays have been used to detect differentially expressed genes (DEGs) and hence elucidate the molecular mechanisms of plant responses to stress [7][8][9][10][11]. However, these methods are time consuming. Recently, high-throughput sequencing (Illumina HiSeq, Roche/454, and ABI SOLiD platform) has become a highly effective tool to identify DEGs, owing to its speed, cost-efficiency, and high-throughput ability. For example, high-throughput sequencing was employed to study the transcriptome profiles of Gossypium hirsutum [12,13] and Barbarea vulgaris [14] after infestation with herbivorous larvae or insects. B. vulgaris can resist the diamond back moth (DBM) and other insects by producing feeding-deterrent triterpenoid saponins. By investigating the DBM-induced changes in the transcriptome, triterpenoid saponin biosynthetic pathways and regulatory networks were analyzed, and the genes involved in these pathways were consequently identified [14]. Furthermore, by comparing transcriptome changes of cotton before and after aphid and whitefly infestation, the expression of some marker genes involved in phytohormonal-mediated plant resistance was found to be suppressed after insect infestation [12]. The suppressed marker genes included cationic peroxidase 3, lipoxygenase I, and nonspecific lipase, which suggest that insects suppress plant resistance in order to facilitate their feeding. In this study, we used Illumina sequencing to compare the levels of differentially expressed genes in leaves of ramie plants with and without RM infestation. The underlying mechanisms involved in ramie defense against insects are subsequently discussed.

Plant Growth, Pest Inoculation, and RNA Preparation.
Ramie Chuanzhu 8, which is moderately resistant to RM, was used in this study. Ramie seedlings were prepared with the shoot-cutting propagation method. The seedlings were cultured in a climate chamber at 26 ∘ C ± 1 ∘ C, 75% ± 1% relative humidity, with a photoperiod of 14 : 10 (L : D). RM egg masses attached to the back of ramie leaves were collected from infested ramie fields at the Institute of Bast Fiber Crops (112.11 E, 28.51 N). For hatching, the eggs were subjected to the same conditions as the ramie seedlings, and J2 larvae were used as inoculums. To prepare the challenged plants (CH), two J2 larvae were transferred to the ramie seedlings onto the fourth leaf from the shoot apex. To stimulate their appetite, BioMed Research International 3 larvae were kept on fasting for 12 h prior to inoculation. At 12,24,48, and 72 h after inoculation, the two topmost undamaged leaves were sampled. Unchallenged control (CK) plants were sampled simultaneously. There were five plants per treatment, and their sampled leaves were pooled per treatment. Overall, eight pooled samples were obtained (CK12, CK24, CK48, CK72, CH12, CH24, CH48, and CH72). Total RNA from each pooled sample was extracted using a commercially available EASYspin plus Total RNA kit (Aidlab, Beijing, China), following the manufacturer's protocol. The obtained RNA was subsequently stored at −80 ∘ C.

Transcriptome Library Preparation and Sequencing.
For transcriptome sequencing of the CK and CH samples, equal amounts of extracted RNA from all four sampling time points were mixed. From these two mixtures, 5 g mixed RNA was used to construct two sequencing libraries, each corresponding to either the CK or the CH sample. From the mixed RNA samples, mRNA was purified by briefly allowing it to bind to magnetic oligo (dT) beads, and it was subsequently broken into short fragments by the addition of a fragmentation buffer (Ambion). The short mRNA fragments were used as templates for synthesizing first-strand cDNA with a random hexamer-primer, dNTPs, RNase H, DNA polymerase I, and GEX. A second-strand buffer was then added to synthesize second-strand cDNA. The resulting cDNA fragments were purified with a QiaQuick PCR extraction kit and subsequently eluted in EB buffer from the kit for end repair and the addition of poly(A) tails. The fragments were then subjected to agarose gel electrophoresis. Suitable fragments were recovered and used as templates for PCR amplification. Finally, the amplified library was sequenced using an Illumina HiSeq™ 2000 platform at Biomarker Technologies Co., LTD, Beijing, China. The resulting data (i.e., the clean reads, counts, and RPKM values) have been submitted to the Gene Expression Omnibus (GEO) at the NCBI, with the GEO accession number GSE66447.

Assembly and Functional Annotation.
In order to yield clean reads, adaptor-only reads, reads containing more than 5% unknown nucleotides, and low-quality reads (reads containing more than 50% of bases with a Q-value of ≤20%) were filtered out of the results after sequencing. To generate nonredundant unigenes, clean reads were subsequently assembled de novo, using the Trinity method with an optimized K-mer length of 25 [15].
For annotation, the unigenes were first searched against the NCBI nonredundant (Nr) database (http://www.ncbi.nlm .nih.gov/), the Swiss-Prot database (http://www.expasy.ch/ sprot), the Clusters of Orthologous Groups (COG) database (http://www.ncbi.nlm.nih.gov/COG/), the Gene Ontology (GO) database (http://geneontology.org/), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) protein database (http://www.genome.jp/kegg) using local BLASTx (with an value < 10 −5 ) to obtain homologous protein information. With Nr annotation, we used the Blast2GO program [16] to obtain GO annotation according to the molecular function, biological process, and cellular component ontology. WEGO software [17] was subsequently used to obtain GO functional classification of all unigenes. In addition, all unigene sequences were aligned to the COG database to predict and classify their possible functions.

Identification of DEGs.
The transcript level of each expressed gene was calculated and normalized to reads per kilobase of exon model per million mapped read (RPKM) [18]. DESeq software (http://www-huber.embl.de/users/anders/ DESeq/) [19] was used to find differentially expressed genes using pairwise comparisons, and the results of all statistical tests were corrected for multiple testing with the Benjamini-Hochberg false discovery rate (FDR) of < 0.01. Sequences were regarded to be significantly differentially expressed if the adjusted value was <0.01, and the absolute value of the log 2 (fold change) was ≥1. Here, the fold change was calculated using the RPKM value of the CH library divided by that of the CK library.

Pathways Enrichment of DEGs.
Pathway enrichment analysis based on the KEGG pathway database (http://www .genome.jp/kegg) was used to identify markedly enriched metabolic pathways or signal transduction pathways in differentially expressed genes, compared with the whole genome background. The following equation was used for the calculations: where is the number of all genes with a KEGG annotation, is the number of DEGs in , is the number of all genes annotated to specific pathways, and is the number of DEGs in .

Quantitative Real-Time PCR (qRT-PCR) Analysis.
To verify genes that were differentially expressed in RM-challenged samples compared with unchallenged ones, qRT-PCR was performed, using an iQ SYBR Green Super Mix kit (Bio-Rad) on an iCycler iQ system (Bio-Rad, Hercules, CA, USA). Gene-specific primers of 13 candidate genes (Table S1 in Supplementary Material available online at http://dx.doi.org/ 10.1155/2016/3702789) were designed using the Primer Premier 5.0 software. The ramie gene encoding actin, which displays a stable expression under different stress conditions [20], was used as an internal control for data normalization. For each sample, first-strand cDNA was synthesized from 1 g of the pooled RNA sample of the CK or CH plants, using a RevertAid First-Strand cDNA Synthesis Kit (Thermo Scientific, Fermentas, Vilnius, Lithuania), according to the manufacturer's instructions. All reactions were performed in triplicate. Expression levels of each gene are presented as the fold change relative to that of the control gene, calculated with the 2 −ΔΔCt method [21].

Sequencing and De Novo Assembly.
To study differently expressed genes in ramie leaves after RM infestation, ramie leaves were subjected to the feeding larvae of C. coerulea Guenée (J2). At 12, 24, 48, and 72 h after inoculation, the  proportion of consumed leaves area by one larva were approximately 5%, 10%, 23%, and 35%, respectively, for each infested leaf (Figure 2). Leaves of CH (challenged) and CK (control) plants were sampled at these time points. RNA extracted from CH and CK samples was used to construct cDNA libraries with a fragment length of 200 bp. Fragments were then sequenced using Illumina paired-end sequencing technology. Raw sequencing data were obtained by basecalling transformation from the sequencing-received raw image data. After filtering, 40,187,616 and 62,784,364 clean reads of 90 bp in length were obtained for the CK and CH libraries, respectively ( Table 1). The GC content of sequence data from the two libraries was 51.29% and 50.73%, respectively. Furthermore, the CycleQ20% was 100% for both, which indicates that the accuracy and quality of the sequencing data were sufficient for further analysis.
De novo assembly was performed using the Trinity method. Sequences that were not extended on either end were defined as being a unigene. Clean reads of CK and CH libraries were assembled into 26,759 (mean length 1,037) and 29,988 unigenes (mean length 982), respectively. To obtain integrated information, reads of the two libraries were merged as well. Finally, 46,533 unigenes for ramie leaves (Chuanzhu 8) were obtained, with an average unigene length of 845 bp (Table 1, Figure S1).

Functional Annotation and Classification.
All of the unigenes assembled were aligned to five public protein databases ( Table 2). For these databases, the alignment resulted in 24,217 homologues (52.04% of all unigenes) with the Nr database, 17,948 (38.57%) with Swiss-Prot, 20,736 (44.56%) with GO, 9,089 (19.53%) with COG, and 5,929 (12.74%) with KEGG. Overall, 24,327 (52.27%) unigenes were functionally annotated, and 11,842 of them were larger than 1,000 bp in length. Approximately half of all unigenes (22,206,47.73%) had no database hits, suggesting that these genes are novel to the databases or specific for ramie.
Based on the functional annotation information, 20,736 unigenes were categorized to 56 different GO terms that followed three main ontology groups: biological process, cellular component, and molecular function ( Figure S2). Because of the matched homologous proteins in databases, some genes were assigned to more than one term. Overall, the 56 GO terms included 221,515 genes. Among these, 113,484 (51.23%) genes were classified into the biological process group and accounted for the greatest proportion, followed by 79,114 (35.71%) in the cellular component and 28,920 (13.06%) in the molecular function group. A total of 136 unigenes were annotated to be involved in biological processes associated with defense responses to insects (Table S2).
All unigenes were aligned to the COG database in order to predict and classify their possible functions. From the 24,217 unigenes (with Nr hits), 9,089 could be categorized according to the COG classifications ( Figure S3). From the 25 categories, the cluster "general function prediction only" (2,314, 18.32%) represented the largest group, followed by "replication, recombination, and repair" (1,216, 9.62%), "translation, ribosomal structure, and biogenesis"  unigenes were assigned to various metabolic processes, and 144 unigenes were assigned to "defense mechanisms." In contrast, only a few unigenes were assigned to "cell motility" (17, 0.13%) and "nuclear structure" (1, 0.01%), and none were assigned to "extracellular structures."

Identification of 1980 DEGs Induced by RM Infestation.
To identify the genes involved in RM responses, the transcriptomes of RM-infested and uninfested samples were compared. The number of clean reads for each gene was calculated. Individual sets of reads were then mapped back to the previously assembled transcript and counted as a proxy for gene expression. Differences in transcript expression between the CH and CK samples were identified with an algorithm developed by Audic and Claverie [22]. As a result, 1,980 genes (750 upregulated and 1,230 downregulated) with at least a twofold difference between the CK and CH treatment were identified as DEGs. According to the annotation information of the unigenes described above, 1,799 of these DEGs were functionally annotated ( Figure 3). Among these DEGs, 87 transcription factors (TFs) belonging to 22 families were found for plants infested with RM (Table 3, Table S3). Furthermore, 16 DEGs involved in the defense response to insects (by GO annotation) were identified for the RMinfested plants (Table S2).

Validation of DEGs by qRT-PCR.
To confirm the gene expression profiles obtained through Illumina sequencing, the expression patterns of 13 candidate genes (Table S1) were further analyzed with qRT-PCR. For the CH samples, the values are presented as the fold change in gene expression, normalized to the reference gene (actin), relative to the CK samples. Overall, eight genes were upregulated, whereas five genes were downregulated in the CH sample ( Figure 4). Although the fold changes in gene expression detected by qRT-PCR were smaller than those detected by Illumina sequencing, the qRT-PCR analysis validated the trend obtained through Illumina sequencing.

Pathway Enrichment Analysis of DEGs.
The RM-affected biological pathways were evaluated by enrichment analysis of DEGs. A total of 97 pathways in ramie were probably affected  by RM infestation (Table S4). Thirteen pathways were found to be significantly affected by DEGs ( < 0.05): four through upregulated DEGs, five by downregulated DEGs, and four enriched by both up-and downregulated DEGs ( Table 4). The most significantly affected pathway in ramie following RM infestation was ribosome ( < 0.0001), which included 84 upregulated and two downregulated DEGs. The second most significantly affected pathway was that of alpha-linolenic acid (ALA) metabolism ( < 0.0001), for which five DEGs were upregulated ( Figure S4). Finally, genes involved in photosynthesis and photosynthesis-antenna proteins pathways were also found to be significantly affected by RM infestation, with 30 of these genes being downregulated.

De Novo Assembly of the Ramie Transcriptome.
Genomewide gene expression profiles can help elucidate the molecular interactions between pests and hosts. In this study, transcriptome changes were analyzed with Illumina sequencing, in order to investigate the defense mechanisms of ramie (the medium-resistant Chuanzhu 8 hybrid) to RM. This resulted in two libraries corresponding to ramie leaves that were RMinfested and uninfested. De novo assembly yielded 46,533 unigenes with an average length of 845 bp. In 2013, the whole-plant transcriptome of ramie Zhongzhu 1 was de novo assembled, using data yielded from Illumina sequencing. A total of 43,990 unigenes with an average length of 824 bp were obtained from approximately 53 million reads [1]. Subsequently, the transcriptome of this ramie was reassembled by adding new data from a ramie-Pratylenchus coffeae interaction study. In that study, a total of 50,486 unigenes with an average length of 853 bp were obtained [23]. More recently, in order to uncover the molecular mechanism behind ramie fiber biosynthesis and development, the transcriptomes of ramie bast, phloem, and xylem were sequenced, and genes involved in cellulose synthesis were studied [24,25]. Our study focused on insect-induced transcriptome changes in ramie. The resulting data add to the available genetic resources for exploring the developmental biology and defense mechanisms of ramie in response to various stresses.

RM Infestation Induced Transcriptome Changes of Ramie.
Recently, there have been studies on changes in ramie transcriptome in response to biotic stresses from P. coffeae [23] and Colletotrichum gloeosporioides [26], or from abiotic stresses such as drought [27], or cadmium addition [28]. These studies provided insights into the molecular basis of ramie tolerance and resistance to stress. However, transcriptome-related information of ramie-RM interactions is still lacking up to date. In our study, by comparing transcriptome profiling of RM-challenged and unchallenged ramie leaves, 1,980 genes were identified as DEGs, of which 750 were upregulated and 1,230 were downregulated. By searching against a public database, 657 upregulated and 1,142 downregulated genes were functionally annotated. Furthermore, the expression patterns of 13 candidate genes, including three TFs genes, nine defense related protein genes, and one tubulin gene, were confirmed using qRT-PCR. Upon RM infestation, more genes in the ramie leaves were found to be downregulated than upregulated, which is similar to observations in B. vulgaris-Plutella xylostella interactions [14].

Transcription Factor Responding to RM.
A transcriptional control of stress, or responsive gene expression, by TFs is crucial in a plant's response to various biotic and abiotic stresses [29,30]. Well-known TFs involved in plant defense responses to pathogens or pest attacks include WRKY, ERF, MYB, bZIP, and NAC [30][31][32][33]. TFs can activate plant defense responses or be involved directly in different defense signaling pathways [34][35][36]. WRKY3 and WRKY6 of tobacco are two insect-responsive WRKY genes; they directly regulate the accumulation of jasmonate (JA) and trypsin proteinase inhibitors and indirectly regulate JA signalingmediated defenses. Silencing of these genes resulted in increased susceptibility to herbivores, indicating that WRKY3 and WRKY6 play important roles in the herbivore resistance of plants [37]. OsERF3, an ERF type TF of rice, was found to mediate resistance to the striped stem borer, mostly likely by suppressing MAPK repressors and modulating JA, SA, ethylene, and H 2 O 2 pathways [38]. Arabidopsis transcription factor AtMYB12 can regulate a number of pathways (including the phenylpropanoid pathway) in transgenic tobacco. This results in an increased accumulation of rutin, and enhanced resistance to Spodoptera litura and Helicoverpa armigera [39]. Furthermore, AtMYB44 can regulate resistance of Arabidopsis to the green peach aphid and diamondback moth by activating EIN2-affected defenses [40]. Previously, comparative transcriptome analysis revealed that, for ramie, drought stress and root lesion nematode infection regulated 24 TF genes (containing three NACs, two MYBs, and one EFR) and 40 TFs genes (containing 10 bHLHs, five MYBs, two NACs, and one ERF), respectively [20,27]. In this study, the expression of 87 TFs was affected by RM larvae infestation (29 were upregulated and 58 were downregulated) and most of them belonged to bHLH (19), EFR (14), ZFP (14), MYB (8), NAC (6), and WRKY (4) families. Overall, we assume that the insect defense mechanisms that involve these TFs are complicated.

Antioxidant Enzymes Responding to RM.
Plants respond to many forms of biotic stress by generating reactive oxygen species (ROS) that participate in defensive signaling and potentiate a hypersensitive response (HR) at the infection site [49]. Immediately after tissue damage caused by herbivore feeding, plants transiently produce ROS (such as the superoxide anion) in the damaged tissue. They additionally produce H 2 O 2 both locally and systemically throughout the plant [50]. After RM infestation, gene expression of two lignin-forming anionic peroxidases (comp30740 c0 and comp32955 c0), two peroxidases (comp24676 c0 and comp33647 c0), one glutathione peroxidase (comp36311 c0), and one catalase (comp191149 c0) was upregulated. Furthermore, gene expressions of two peroxidases (comp30251 c0 and comp31888 c0) were downregulated. Peroxidase (POD), catalase (CAT), and glutathione peroxidase (GPX) are well-known oxidative stress-related proteins that participate in ROS metabolism [51,52]. POD is important in ROS generation, whereas CAT and GPX are important in ROS scavenging [52,53]. A number of studies suggested that POD and CAT play important roles in plant resistance to herbivory [54,55]. These enzymes were also reported to be involved in ramie resistance to nematodes and anthracnose fungi [26,56]. In addition, studies also showed that GPX increases plant tolerance to different stress types and regulates cellular immune responses [57,58]. Moderate amounts of ROS are beneficial, but increased levels of ROS can inhibit cell migration and proliferation. It can even cause severe tissue damage. Therefore, cells must develop strategies for the detoxification of these molecules [59]. In general, plants appear to balance generating ROS as a defensive mechanism and producing ROS-detoxifying enzymes to cope with their own oxidative damage [60]. Hence, oxidative stress-related genes are not uniformly regulated (some upregulated and others downregulated) in leaf tissue in response to insect infestation [61].

ALA Metabolism Influenced by RM Infestation.
ALA metabolism is one of the pathways that was significantly affected by RM infestation; the end products of this pathway are (−)-methyl-jasmonate (MeJA) or (+)-7-isomethyljasmonate. Upregulation of five genes in this pathway leads to accumulation of JA. JA accumulation usually triggers the biosynthesis of phytoalexins and affects the expression of PR genes synergistically, or by antagonizing the action of other plant hormones [62]. In higher plants, JA is synthesized mainly via the octadecanoid pathway, and lipoxygenase (LOX) is a key enzyme that oxidizes alpha-linolenic acid in the first step of this pathway [63]. In RM-infested ramie, two potential LOX encoding genes (comp29491 c1 and comp29491 c2) were found to be upregulated, and two genes (comp37093 c0 and comp25237 c0) were downregulated. The octadecanoid pathway, however, was not significantly affected. Additionally, the efficiency of the JA mediated defense response varied in different herbivore-plant interaction systems and either decreased [64] or increased resistance [37]. Therefore, the role of JA and JA regulation in ramie-RM interactions needs to be investigated in future studies.

Conclusions
In summary, to our knowledge, this is the first study on transcriptome changes in ramie leaves induced by insects (RM larvae). The transcriptomes of RM-infested and uninfested ramie leaves were sequenced using Illumina sequencing, and 29,988 and 26,759 unigenes were obtained, respectively. Comparative transcriptome analysis showed that the expression of 1980 genes changed substantially in ramie leaves in response to infestation with RM larvae. Differential expression of some genes was confirmed by qRT-PCR. Many genes encoding for TFs, PIs, and antioxidant enzymes, which are thought to be involved in insect resistance, were regulated by RM larvae infestation. Further characterization of the targets of these genes will help us understand the details of ramie insect resistance.