BioMed Research International
Volume 2016 (2016), Article ID 3702789, 10 pages
http://dx.doi.org/10.1155/2016/3702789

Research Article

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

Liangbin Zeng,¹ Airong Shen,² Jia Chen,¹ Zhun Yan,¹ Touming Liu,¹ Zhaodong Xue,¹ and Yongting Yu¹

¹Institute of Bast Fiber Crops, Chinese Academy of Agricultural Sciences, Changsha, Hunan 410205, China
²Hunan Academy of Forestry Science, Changsha, Hunan 410004, China

Received 5 October 2015; Revised 13 December 2015; Accepted 1 February 2016

Academic Editor: Pulugurtha B. Kirti

Copyright © 2016 Liangbin Zeng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

1. 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.

(a)

(b)

(c)

(d)

Figure 1: Pictures of Cocytodes coerulea Guenée larvae and their resulting damage to ramie leaves. (a) J2 stage larvae feeding on ramie leaves. (b) The resulting net-like structure for the ramie leaves after J2 stage herbivory of C. coerulea. (c) A J5 stage larva eating a ramie leaf. (d) The remaining stems and main veins after infestation of the J5 stage larvae.

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–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.

2. Materials and Methods

2.1. 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, 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.

2.2. 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.

2.3. 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⁻⁵) 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.

2.4. 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 . Sequences were regarded to be significantly differentially expressed if the adjusted value was <0.01, and the absolute value of the log₂ (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.

2.5. 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 .

2.6. 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 method [21].

3. Results

3.1. 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 base-calling 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.

Table 1: Summary of sequencing and assembly results.

CK
CH
Total

Total clean reads
40,187,616
62,784,364

Total clean nucleotides (nt)
4,149,401,262
6,340,667,793

GC percentage (%)
51.29
50.73

percentage (%)
0.04
0.04

Q20 percentage (%)
90.39
90.35

Contigs

Total number
2,704,882
4,239,322
3,582,696

Total length (nt)
129,385,682
185,382,731
208,724,373

N50 length (nt)
47
43
49

Mean length (nt)
48
44
58

Transcripts

Total number
52,165
58,093
75,553

Total length (nt)
66,692,367
72,503,398
88,128,254

N50 length (nt)
1,846
1,849
1933

Mean length (nt)
1,278
1,248
1166

Unigenes

Total number
26,759
19,988
46,533

Total length (nt)
27,761,126
29,962,009
39,341,510

N50 length (nt)
1,708
1,678
1,585

Mean length (nt)
1,037
982
845

Figure 2: The proportion of leaf consumed per larva (J2 stage) at different time after infestation.

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).

3.2. 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.

Table 2: Number of unigenes annotated in five public databases.

Database
Annotated number
300 ≤ length < 1000
Length ≥ 1000

Nr
24,217
8,771
11,840

Swiss-Prot
17,948
6,299
9,198

GO
20,736
7,115
10,781

COG
9,089
2,608
5,394

KEGG
5,929
2,008
2,914

Total
24,327
8,808
11,842

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” (1,106, 8.75%), “transcription” (1,085, 8.59%), and “signal transduction mechanisms” (938, 7.42%). Collectively, 3,863 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.”

3.3. 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 RM-infested plants (Table S2).

Table 3: Summary of differentially expressed genes annotated as transcription factors.

Gene family
Upregulated
Downregulated
DEG number

AP2
1
0
1

ARF
1
1
2

bHLH
5
14
19

bZip
1
0
1

C2H2
0
2
2

CAMTA
0
1
1

COL
0
1
1

EFR
5
9
14

GATA
0
2
2

GLK
0
1
1

TRY
1
0
1

HSF
1
1
2

HY5
0
1
1

MYB
3
5
8

NAC
5
1
6

ORG
0
1
1

RF2
0
3
3

TGA3
0
1
1

TT2
0
1
1

WRKY
2
2
4

YABBY
0
1
1

ZFP
4
10
14

Total
29
58
87

Figure 3: Numbers of differentially expressed genes (DEGs) annotated in five public databases.

3.4. 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.

Figure 4: Validation of the gene expression results obtained from Illumina sequencing by qRT-PCR. The left and right vertical ordinates indicate normalized gene expression determined by RNA sequencing and qRT-PCR, respectively.

3.5. 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 (): 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 (), which included 84 upregulated and two downregulated DEGs. The second most significantly affected pathway was that of alpha-linolenic acid (ALA) metabolism (), 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.

Table 4: List of pathways significantly enriched in differentially expressed genes ().

Pathway term
DEGs tested
value
Pathway ID

Pathways for both upregulated and downregulated DEGs

Ribosome
84 (23.46%)

ko03010

Glycosphingolipid biosynthesis-globo series
4 (1.12%)

ko00603

Zeatin biosynthesis
6 (1.68%)

ko00908

Ubiquinone and other terpenoid-quinone biosyntheses
7 (1.96%)

ko00130

Pathways for upregulated DEGs

Alpha-linolenic acid metabolism
5 (1.4%)

ko00592

Valine, leucine, and isoleucine degradation
13 (3.63%)

ko00280

Arachidonic acid metabolism
4 (1.12%)

ko00590

Isoquinoline alkaloid biosynthesis
4 (1.12%)

ko00950

Pathways for downregulated DEGs

Carotenoid biosynthesis
8 (2.23%)

ko00906

Photosynthesis
21 (5.87%)

ko00195

Photosynthesis: antenna proteins
9 (2.51%)

ko00196

Sulfur metabolism
7 (1.96%)

ko00920

Pentose and glucuronate interconversions
10 (2.79%)

ko00040

4. Discussion

4.1. 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 RM-infested 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.

4.2. 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].

4.3. 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–33]. TFs can activate plant defense responses or be involved directly in different defense signaling pathways [34–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 signaling-mediated 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₂O₂ 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.

4.4. Protease Inhibitors Responding to RM

Insects feed directly on plant tissue, and the plant endogenous enzyme inhibitors obstruct this digestion by influencing the digestive amylases and proteases present in the insect guts [41]. Transgenic research showed that some endogenous protease inhibitors (PIs) of plants displayed antiherbivore characteristics. For example, bean α-amylase inhibitor 1 gives transgenic cowpeas an increased resistance to two storage pests, Callosobruchus maculatus and C. chinensis, which cause severe damage to cowpea seeds during storage [42]. Furthermore, the barley trypsin inhibitor (CMe) increases the resistance of transgenic rice to Sitophilus oryzae [43], and trypsin inhibitors from other plants have shown to increase resistance to insects [44–48]. In our study, 14 genes encoding for different protease inhibitors were found to respond to RM larvae infestation. These 14 genes included trypsin inhibitors (comp15649_c0, comp27345_c0, comp30853_c0, and comp38186_c0), alpha-amylase/subtilisin inhibitors (comp28610_c0 and comp29690_c0), polygalacturonase inhibitors (comp24833_c0 and comp27076_c0), and a xyloglucanase inhibitor (comp36158_c1). Three of these genes (comp15649_c0, comp39830_c0, and comp36636_c0) were also annotated by GO analysis for involvement in the insect defense response in the biological process group. Additionally, several genes were upregulated following RM infestation. This included four genes encoding trypsin inhibitors, such as comp15649_c0 (upregulated 18.6-fold), comp27345_c0 (4.7-fold), comp30853_c0 (4.4-fold), and comp38186_c0 (7.3-fold), as well as two alpha-amylase/subtilisin inhibitor encoding genes, being comp28610_c0 (8.7-fold) and comp29690_c0 (2.6-fold). These results imply that these genes play active roles in ramie-pest defense.

4.5. 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₂O₂ 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].

4.6. 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-isomethyl-jasmonate. 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.

5. 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.

Conflict of Interests

All authors have declared that no competing interests exist.

Acknowledgments

This study was supported financially by the Earmarked Fund for China Agriculture Research System (CARS-19-E19) and the National Natural Science Foundation of China (31201494 and 31501653). The authors thank Editage Company for its linguistic assistance during the preparation of this paper.

References

T. M. Liu, S. Y. Zhu, Q. M. Tang, P. Chen, Y. Yu, and S. Tang, “De novo assembly and characterization of transcriptome using Illumina paired-end sequencing and identification of CesA gene in ramie (Boehmeria nivea L. Gaud),” BMC Genomics, vol. 14, article 125, 2013.

C. M. Yu, J. R. Chen, Y. Z. Wang et al., “Research progresses of molecular breeding in ramie and forage Ramie,” Plant Fiber Sciences in China, vol. 29, pp. 389–392, 2007.

T. Guo, W. She, C. X. Xiao et al., “Research progress of ramie for forage,” Crop Research, vol. 6, pp. 730–733, 2012.

J.-Y. Ide, “Inter- and intra-shoot distributions of the ramie moth caterpillar, Arcte coerulea (Lepidoptera: Noctuidae), in ramie shrubs,” Applied Entomology and Zoology, vol. 41, no. 1, pp. 49–55, 2006.

L. B. Zeng, Z. D. Xue, Y. T. Yu et al., “Changes of occurrence and control of Cocytodes caerulea Guenée,” Hunan Agricultural Sciences, vol. 10, pp. 23–24, 2013.

L.-B. Zeng, Z.-D. Xue, Z.-Y. Yan et al., “Relationship between main chemical compositions in leaves of ramie and resistance to Cocytodes coerulea,” Hunan Agricultural Science, vol. 19, pp. 77–78, 2011.

M. Schena, D. Shalon, R. W. Davis, and P. O. Brown, “Quantitative monitoring of gene expression patterns with a complementary DNA microarray,” Science, vol. 270, no. 5235, pp. 467–470, 1995.

C. W. B. Bachem, R. S. van der Hoeven, S. M. de Bruijn, D. Vreugdenhil, M. Zabeau, and R. G. F. Visser, “Visualization of differential gene expression using a novel method of RNA fingerprinting based on AFLP: analysis of gene expression during potato tuber development,” Plant Journal, vol. 9, no. 5, pp. 745–753, 1996.

L. Diatchenko, Y.-F. C. Lau, A. P. Campbell et al., “Suppression subtractive hybridization: a method for generating differentially regulated or tissue-specific cDNA probes and libraries,” Proceedings of the National Academy of Sciences of the United States of America, vol. 93, no. 12, pp. 6025–6030, 1996.

X. Wang, W. Liu, X. Chen et al., “Differential gene expression in incompatible interaction between wheat and stripe rust fungus revealed by cDNA-AFLP and comparison to compatible interaction,” BMC Plant Biology, vol. 10, article 9, 2010.

A. Moumeni, K. Satoh, H. Kondoh et al., “Comparative analysis of root transcriptome profiles of two pairs of drought-tolerant and susceptible rice near-isogenic lines under different drought stress,” BMC Plant Biology, vol. 11, article 174, 2011.

N. K. Dubey, R. Goel, A. Ranjan et al., “Comparative transcriptome analysis of Gossypium hirsutum L. in response to sap sucking insects: aphid and whitefly,” BMC Genomics, vol. 14, article 241, 2013.

S. Artico, M. Ribeiro-Alves, O. Oliveira-Neto et al., “Transcriptome analysis of Gossypium hirsutum flower buds infested by cotton boll weevil (Anthonomus grandis) larvae,” BMC Genomics, vol. 15, article 854, 2014.

X. Wei, X. Zhang, D. Shen et al., “Transcriptome analysis of Barbarea vulgaris infested with diamondback moth (Plutella xylostella) larvae,” PLoS ONE, vol. 8, no. 5, Article ID e64481, 2013.

M. G. Grabherr, B. J. Haas, M. Yassour et al., “Full-length transcriptome assembly from RNA-Seq data without a reference genome,” Nature Biotechnology, vol. 29, no. 7, pp. 644–652, 2011.

A. Conesa, S. Götz, J. M. García-Gómez, J. Terol, M. Talón, and M. Robles, “Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research,” Bioinformatics, vol. 21, no. 18, pp. 3674–3676, 2005.

J. Ye, L. Fang, H. Zheng et al., “WEGO: a web tool for plotting GO annotations,” Nucleic Acids Research, vol. 34, pp. W293–W297, 2006.

A. Mortazavi, B. A. Williams, K. McCue, L. Schaeffer, and B. Wold, “Mapping and quantifying mammalian transcriptomes by RNA-Seq,” Nature Methods, vol. 5, no. 7, pp. 621–628, 2008.

S. Anders and W. Huber, “Differential expression analysis for sequence count data,” Genome Biology, vol. 11, article R106, 2010.

T. M. Liu, S. Y. Zhu, Q. M. Tang, and S. Tang, “Identification of 32 full-length NAC transcription factors in ramie (Boehmeria nivea L. Gaud) and characterization of the expression pattern of these genes,” Molecular Genetics and Genomics, vol. 289, no. 4, pp. 675–684, 2014.

K. J. Livak and T. D. Schmittgen, “Analysis of relative gene expression data using real-time quantitative PCR and the method,” Methods, vol. 25, no. 4, pp. 402–408, 2001.

S. Audic and J.-M. Claverie, “The significance of digital gene expression profiles,” Genome Research, vol. 7, no. 10, pp. 986–995, 1997.

S. Y. Zhu, S. W. Tang, Q. M. Tang, and T. Liu, “Genome-wide transcriptional changes of ramie (Boehmeria nivea L. Gaud) in response to root-lesion nematode infection,” Gene, vol. 552, no. 1, pp. 67–74, 2014.

J. R. Chen, F. Liu, Y. H. Tang, Y. Yuan, and Q. Guo, “Transcriptome sequencing and profiling of expressed genes in phloem and xylem of ramie (Boehmeria nivea L. Gaud),” PLoS ONE, vol. 9, no. 10, Article ID e110623, 2014.

J. Chen, Z. Pei, L. Dai et al., “Transcriptome profiling using pyrosequencing shows genes associated with bast fiber development in ramie (Boehmeria nivea L.),” BMC Genomics, vol. 15, article 919, 2014.

X. X. Wang, J. Chen, B. Wang et al., “Characterization by suppression subtractive hybridization of transcripts that are differentially expressed in leaves of anthracnose-resistant ramie cultivar,” Plant Molecular Biology Reporter, vol. 30, no. 3, pp. 547–555, 2012.

T. M. Liu, S. Zhu, Q. Tang, Y. Yu, and S. Tang, “Identification of drought stress-responsive transcription factors in ramie (Boehmeria nivea L. Gaud),” BMC Plant Biology, vol. 13, no. 1, article 130, 2013.

T. M. Liu, S. Y. Zhu, Q. M. Tang et al., “Genome-wide transcriptomic profiling of ramie (Boehmeria nivea L. Gaud) in response to cadmium stress,” Gene, vol. 558, no. 1, pp. 131–137, 2015.

W. Chon, N. J. Provart, J. Glazebrook et al., “Expression profile matrix of Arabidopsis transcription factor genes suggests their putative functions in response to environmental stresses,” Plant Cell, vol. 14, no. 3, pp. 559–574, 2002.

K. B. Singh, R. C. Foley, and L. Oñate-Sánchez, “Transcription factors in plant defense and stress responses,” Current Opinion in Plant Biology, vol. 5, no. 5, pp. 430–436, 2002.

M. Jakoby, B. Weisshaar, W. Dröge-Laser et al., “bZIP transcription factors in Arabidopsis,” Trends in Plant Science, vol. 7, no. 3, pp. 106–111, 2002.

A. N. Olsen, H. A. Ernst, L. L. Leggio, and K. Skriver, “NAC transcription factors: structurally distinct, functionally diverse,” Trends in Plant Science, vol. 10, no. 2, pp. 79–87, 2005.

P. J. Rushton, I. E. Somssich, P. Ringler, and Q. J. Shen, “WRKY transcription factors,” Trends in Plant Science, vol. 15, no. 5, pp. 247–258, 2010.

Y.-Q. Gu, M. C. Wildermuth, S. Chakravarthy et al., “Tomato transcription factors Pti4, Pti5, and Pti6 activate defense responses when expressed in Arabidopsis,” Plant Cell, vol. 14, no. 4, pp. 817–831, 2002.

M. Kalde, M. Barth, I. E. Somssich, and B. Lippok, “Members of the Arabidopsis WRKY group III transcription factors are part of different plant defense signaling pathways,” Molecular Plant-Microbe Interactions, vol. 16, no. 4, pp. 295–305, 2003.

Q. Bu, H. Jiang, C.-B. Li et al., “Role of the Arabidopsis thaliana NAC transcription factors ANAC019 and ANAC055 in regulating jasmonic acid-signaled defense responses,” Cell Research, vol. 18, no. 7, pp. 756–767, 2008.

M. Skibbe, N. Qu, I. Galis, and I. T. Baldwin, “Induced plant defenses in the natural environment: Nicotiana attenuata WRKY3 and WRKY6 coordinate responses to herbivory,” Plant Cell, vol. 20, no. 7, pp. 1984–2000, 2008.

J. Lu, H. Ju, G. Zhou et al., “An EAR-motif-containing ERF transcription factor affects herbivore-induced signaling, defense and resistance in rice,” Plant Journal, vol. 68, no. 4, pp. 583–596, 2011.

P. Misra, A. Pandey, M. Tiwari et al., “Modulation of transcriptome and metabolome of tobacco by Arabidopsis transcription factor, AtMYB12, leads to insect resistance,” Plant Physiology, vol. 152, no. 4, pp. 2258–2268, 2010.

B.-B. Lü, X.-J. Li, W.-W. Sun et al., “AtMYB44 regulates resistance to the green peach aphid and diamondback moth by activating EIN2-affected defences in Arabidopsis,” Plant Biology, vol. 15, no. 5, pp. 841–850, 2013.

O. L. Franco, D. J. Rigden, F. R. Melo, and M. F. Grossi-de-Sá, “Plant α-amylase inhibitors and their interaction with insect α-amylases,” European Journal of Biochemistry, vol. 269, no. 2, pp. 397–412, 2002.

S. K. Solleti, S. Bakshi, J. Purkayastha, S. K. Panda, and L. Sahoo, “Transgenic cowpea (Vigna unguiculata) seeds expressing a bean α-amylase inhibitor 1 confer resistance to storage pests, bruchid beetles,” Plant Cell Reports, vol. 27, no. 12, pp. 1841–1850, 2008.

J. Alfonso-Rubí, F. Ortego, P. Castañera, P. Carbonero, and I. Díaz, “Transgenic expression of trypsin inhibitor CMe from barley in indica and japonica rice, confers resistance to the rice weevil Sitophilus oryzae,” Transgenic Research, vol. 12, no. 1, pp. 23–31, 2003.

F. De Leo, M. Bonadé-Bottino, L. R. Ceci, R. Gallerani, and L. Jouanin, “Effects of a mustard trypsin inhibitor expressed in different plants on three lepidopteran pests,” Insect Biochemistry and Molecular Biology, vol. 31, no. 6-7, pp. 593–602, 2001.

T. Srinivasan, K. R. R. Kumar, and P. B. Kirti, “Constitutive expression of a trypsin protease inhibitor confers multiple stress tolerance in transgenic tobacco,” Plant and Cell Physiology, vol. 50, no. 3, pp. 541–553, 2009.

M. L. Rodrigues Macedo, R. A. Durigan, D. S. da Silva, S. Marangoni, M. D. G.-A. M. Freire, and J. R. P. Parra, “Adenanthera pavonina trypsin inhibitor retard growth of Anagasta kuehniella (Lepidoptera: Pyralidae),” Archives of Insect Biochemistry and Physiology, vol. 73, no. 4, pp. 213–231, 2010.

C. F. R. de Oliveira, S. Marangoni, and M. L. R. Macedo, “The trypsin inhibitor from Entada acaciifolia seeds affects negatively the development of Mediterranean flour moth, Anagasta kuehniella,” Pesticide Biochemistry and Physiology, vol. 108, no. 1, pp. 74–79, 2014.

A. Mittal, R. Kansal, V. Kalia, M. Tripathi, and V. K. Gupta, “A kidney bean trypsin inhibitor with an insecticidal potential against Helicoverpa armigera and Spodoptera litura,” Acta Physiologiae Plantarum, vol. 36, no. 2, pp. 525–539, 2014.

C. Lamb and R. A. Dixon, “The oxidative burst in plant disease resistance,” Annual Review of Plant Physiology and Plant Molecular Biology, vol. 48, pp. 251–275, 1997.

A. Kessler and I. T. Baldwin, “Plant responses to insect herbivory: the emerging molecular analysis,” Annual Review of Plant Biology, vol. 53, no. 1, pp. 299–328, 2002.

F. J. Corpas, J. B. Barroso, and L. A. del Río, “Peroxisomes as a source of reactive oxygen species and nitric oxide signal molecules in plant cells,” Trends in Plant Science, vol. 6, no. 4, pp. 145–150, 2001.

R. Mittler, S. Vanderauwera, M. Gollery, and F. Van Breusegem, “Reactive oxygen gene network of plants,” Trends in Plant Science, vol. 9, no. 10, pp. 490–498, 2004.

T. Kawano, “Roles of the reactive oxygen species-generating peroxidase reactions in plant defense and growth induction,” Plant Cell Reports, vol. 21, no. 9, pp. 829–837, 2003.

G. K. Taggar, R. S. Gill, A. K. Gupta, and J. S. Sandhu, “Fluctuations in peroxidase and catalase activities of resistant and susceptible black gram (Vigna mungo (L.) Hepper) genotypes elicited by Bemisia tabaci (Gennadius) feeding,” Plant Signaling and Behavior, vol. 7, no. 10, pp. 1321–1329, 2012.

Z. Wei, W. Hu, Q. Lin et al., “Understanding rice plant resistance to the Brown Planthopper (Nilaparvata lugens): a proteomic approach,” Proteomics, vol. 9, no. 10, pp. 2798–2808, 2009.

Y. T. Yu, L. B. Zeng, Z. Yan et al., “Identification of ramie genes in response to Pratylenchus coffeae infection challenge by digital gene expression analysis,” International Journal of Molecular Sciences, vol. 16, no. 9, pp. 21989–22007, 2015.

C. C. C. Chang, I. Šlesak, L. Jordá et al., “Arabidopsis chloroplastic glutathione peroxidases play a role in cross talk between photooxidative stress and immune ,” Plant Physiology, vol. 150, no. 2, pp. 670–683, 2009.

K. Yoshimura, K. Miyao, A. Gaber et al., “Enhancement of stress tolerance in transgenic tobacco plants overexpressing Chlamydomonas glutathione peroxidase in chloroplasts or cytosol,” The Plant Journal, vol. 37, no. 1, pp. 21–33, 2004.

H. Steiling, B. Munz, S. Werner, and M. Brauchle, “Different types of ROS-scavenging enzymes are expressed during cutaneous wound repair,” Experimental Cell Research, vol. 247, no. 2, pp. 484–494, 1999.

K. Zhu-Salzman, R. A. Salzman, J.-E. Ahn, and H. Koiwa, “Transcriptional regulation of sorghum defense determinants against a phloem-feeding aphid,” Plant Physiology, vol. 134, no. 1, pp. 420–431, 2004.

G. A. Thompson and F. L. Goggin, “Transcriptomics and functional genomics of plant defence induction by phloem-feeding insects,” Journal of Experimental Botany, vol. 57, no. 4, pp. 755–766, 2006.

K. A. Aliferis, D. Faubert, and S. Jabaji, “A metabolic profiling strategy for the dissection of plant defense against fungal pathogens,” PLoS ONE, vol. 9, no. 11, Article ID e111930, 2014.

G. A. Howe and G. Jander, “Plant immunity to insect herbivores,” Annual Review of Plant Biology, vol. 59, pp. 41–66, 2008.

J. S. Thaler, “Jasmonate-inducible plant defences cause increased parasitism of herbivores,” Nature, vol. 399, no. 6737, pp. 686–688, 1999.