Skip to main content
  • Research article
  • Open access
  • Published:

Comparative transcriptome analysis of soybean response to bean pyralid larvae

Abstract

Background

Soybean is one of most important oilseed crop worldwide, however, its production is often limited by many insect pests. Bean pyralid is one of the major soybean leaf-feeding insects in China. To explore the defense mechanisms of soybean resistance to bean pyralid, the comparative transcriptome sequencing was completed between the leaves infested with bean pyralid larvae and no worm of soybean (Gantai-2-2 and Wan82–178) on the Illumina HiSeq™ 2000 platform.

Results

In total, we identified 1744 differentially expressed genes (DEGs) in the leaves of Gantai-2-2 (1064) and Wan82–178 (680) fed by bean pyralid for 48 h, compared to 0 h. Interestingly, 315 DEGs were shared by Gantai-2-2 and Wan82–178, while 749 and 365 DEGs specifically identified in Gantai-2-2 and Wan82–178, respectively. When comparing Gantai-2-2 with Wan82–178, 605 DEGs were identified at 0 h feeding, and 468 DEGs were identified at 48 h feeding. Gene Ontology (GO) annotation analysis revealed that the DEGs were mainly involved in the metabolic process, single-organism process, cellular process, responses to stimulus, catalytic activities and binding. Pathway analysis showed that most of the DEGs were associated with the plant-pathogen interaction, phenylpropanoid biosynthesis, phenylalanine metabolism, flavonoid biosynthesis, peroxisome, plant hormone signal transduction, terpenoid backbone biosynthesis, and so on. Finally, we used qRT-PCR to validate the expression patterns of several genes and the results showed an excellent agreement with deep sequencing.

Conclusions

According to the comparative transcriptome analysis results and related literature reports, we concluded that the response to bean pyralid feeding might be related to the disturbed functions and metabolism pathways of some key DEGs, such as DEGs involved in the ROS removal system, plant hormone metabolism, intracellular signal transduction pathways, secondary metabolism, transcription factors, biotic and abiotic stresses. We speculated that these genes may have played an important role in synthesizing substances to resist insect attacks in soybean. Our results provide a valuable resource of soybean defense genes that will benefit other studies in this field.

Background

Soybean (Glycine max (L.) Merr.) is the largest oil crop worldwide, and is widely used in the production of food, feed, industrial products and other sideline fields [1]. However, there have been large increases in soybean production costs due to pests. Bean pyralid (Lamprosema indicate (Fabricius)) is one of the major leaf-feeding insects that affects soybean crops in central and southern China, the larvae lurk inside soybean leaves, cause leaf curling and feed on leaf tissues, affecting photosynthesis. Therefore, the plants cannot grow normally [2]. So bean pyralid is different from other leaf-feeding insects with chewing mouthparts which cause holes or incisions by means of encroachment [3]. In the soybean-producing areas of southern China, many generations of bean pyralids may appear in 1 year. In serious pest-damaged years, only veins and petioles will be left, causing serious yield losses [4]. Sun et al. and Long et al. evaluated rolled leaflet number and larva number could be used as an evaluation index for bean pyralid in soybean, and screened the highly resistant line Gantai-2-2 and the highly sensitive line Wan 82–178 [4,5,6]. Two indicators of resistance to bean pyralid, rolled leaflet number and rolled leaflet percentage, were a significantly positively correlated with the pubescence angle, length on leaf blade, angle on petiole and a significantly negatively correlated with the pubescence density on leaf blade, but on correlation with pubescence tip shape was observed [7]. Xing et al. and Li et al. showed that soybean resistance to bean pyralid accords with two or three major genes and polygene, 81–92% of the phenotypic variation was accounted for by additive quantitative trait locus (QTL) (27–43%), epistatic QTL pairs (5–13%) and collective unmapped minor QTL (38–58%) [8,9,10]. The contents of soluble sugar, superoxide dismutase (SOD), polyphenol oxidase (PPO), jasmonate (JA) and abscisic acid (ABA) are significantly increased after bean pyralid feeding [11]. However, the results of comparative transcriptome research which has focused on soybeans’ resistance to bean pyralid has not yet been made available. This is the first study of soybean transcriptome in response to bean pyralid feeding.

Transcriptome sequencing has become an important method for gene expression analysis, differentially expressed genes (DEGs) selection, functional gene mining, and genetic evolution analysis. The soybean genome was released in 2010 [12]. Based on the soybean genome data and transcriptome technology, can be better to examine all the transcription reactions, structural functions, and transcriptional regulation of resistant soybean varieties at the overall level. In the present paper, we tried to find important DEGs and metabolism pathways might related to the soybean in response to bean pyralid larvae through the comparative transcrptome analysis between the leaves infested with bean pyralid larvae and no worm of soybean using the Illumina HiSeq™ 2000 platform. Our results provide a valuable resource of soybean defense genes that will benefit other studies in this field.

Results

Transcriptome sequencing and sequence alignment

An Illumina HiSeq™2000 sequencer was employed to analyze the comparative transcriptome of eight samples of Gantai-2-2 and Wan82–178 leaves that bean pyralid had been feeding on 0 h and 48 h. The original image data obtained by sequencing base-calling were the original sequence reads. Each read in the Solexa paired-end (PE) sequencing was 100 bp in length. There were 45.88 G original data sets produced during sequencing. The mean sequencing depth was 5.67. After the raw data were trimmed, 442,422,398 clean reads were obtained. The clean/raw read rates of the eight samples ranged from 95.30% to 97.11%. All clean reads were matched to the soybean reference genome by BWA software, allowing two base mismatches. The mapped genome reads ranged from 42,474,863 to 44,489,050 sets, genome map rates ranged from 77.23% to 78.95%, and unique match rates ranged from 72.86% to 75.22%. The expressed genes ranged from 50,283 to 53,739 (Table 1). To estimate whether the sequencing depth was sufficient for transcriptome coverage, the sequencing saturation in the eight cDNA libraries was analyzed. The results showed that most genes became saturated when the amount of PE reads was 20 M (200 × 100 kb) (Fig. 1), which indicated that the overall quality of sequencing saturation in the eight cDNA libraries was high and that the sequencing amount covered the vast majority of expressed genes.

Table 1 Number of reads sequenced and mapped to soybean genome
Fig. 1
figure 1

Analysis of sequencing saturation. a HRK0–1, bHRK0–2, c HRK48–1, d HRK48–2, e HSK0–1, f HSK0–2, g HSK48–1, h HSK48–2

Correlation analysis of samples

The correlation of gene expression levels among samples is a key criterion to test whether the experiments are reliable and whether the chosen samples are reasonable. If one sample is highly similar to another one, the correlation value between them is very close to 1. We calculated the correlation value between each of two samples based on the FPKM results. According to the standard that Encode plan recommends, the square of the correlation value (R2) should be ≥0.92 (under an ideal experimental environment with reasonable samples). Our results showed that the R2 of all repetitions were >0.95 (Fig. 2), which signified that our experimental samples and results were satisfactory and reliable.

Fig. 2
figure 2

Correlations value of each repetition. a HRK0–1 and HRK0–2. b HRK48–1 and HRK48–2. c HSK0–1 and HSK0–2. d HSK48–1 and HSK48–2

Screening of differentially expressed genes (DEGs)

Noiseq, DESeq2 and edgeR methods were used to screen DEGs (Fig. 3). The results showed that the Noiseq method can screen DEGs between two groups, showing a good performance when comparing it to other differential expression methods, such as edgeR and DESeq2. Noiseq maintains good True Positive and False Positive rates when the sequencing depth is increased, whereas most other methods show poor performance. Further more, Noiseq models the noise distribution from the actual data to better adapt to the size of the data set and be more effective in controlling the rate of false discoveries. Therefore, the Noiseq method was used to screen the DEGs.

Fig. 3
figure 3

The DEGs were screened by Noiseq, DESeq2 and edgeR. a HRK48/HRK0_UP In total, 894, 900 and 1050 up-regulated DEGs were identified by Noiseq, DESeq2 and edgeR, respectively. 460 DEGs were identified under the three methods, 62 DEGs were identified under both Noiseq and edgeR, 388 DEGs were identified under both DESeq2 and edgeR, 10 DEGs were identified under both Noiseq and DESeq2. b HRK48/HRK0_DOWN In total, 170, 991 and 1028 down-regulated DEGs were identified by Noiseq, DESeq2 and edgeR, respectively. 85 DEGs were identified under the three methods, 12 DEGs were identified under both Noiseq and edgeR, 771 DEGs were identified under both DESeq2 and edgeR, 1 DEGs were identified under both Noiseq and DESeq2. c HSK48/HSK0_UP In total, 495, 595 and 448 up-regulated DEGs were identified by Noiseq, DESeq2 and edgeR, respectively. 210 DEGs were identified under the three methods, 29 DEGs were identified under both Noiseq and edgeR, 196 DEGs were identified under both DESeq2 and edgeR, 15 DEGs were identified under both Noiseq and DESeq2. d HSK48/HSK0_DOWN In total, 185, 434 and 183 down-regulated DEGs were identified by Noiseq, DESeq2 and edgeR, respectively. 47 DEGs were identified under the three methods, 9 DEGs were identified under both Noiseq and edgeR, 122 DEGs were identified under both DESeq2 and edgeR, 2 DEGs were identified under both Noiseq and DESeq2. e HRK0/HSK0_UP In total, 192, 264 and 147 up-regulated DEGs were identified by Noiseq, DESeq2 and edgeR, respectively. 84 DEGs were identified under both Noiseq and DESeq2. f HRK0/HSK0_DOWN In total, 413, 241 and 116 down-regulated DEGs were identified by Noiseq, DESeq2 and edgeR, respectively. 120 DEGs were identified under both Noiseq and DESeq2. g HRK48/HSK48_UP In total, 202, 100 and 146 up-regulated DEGs were identified by Noiseq, DESeq2 and edgeR, respectively. 56 DEGs were identified under both Noiseq and DESeq2. h HRK48/HSK48_DOWN In total, 266, 131 and 121 down-regulated DEGs were identified by Noiseq, DESeq2 and edgeR, respectively. 74 DEGs were identified under both Noiseq and DESeq2

As a result, under bean pyralid larvae feeding for 48 h, 1064 DEGs were identified in the Gantai-2-2, of which 894 DEGs were up-regulated and 170 DEGs were down-regulated compared to 0 h (Additional file 1: Table S1). Additionally, 680 DEGs were identified in Wan82–178, of which 495 DEGs were up-regulated and 185 DEGs were down-regulated (Additional file 2: Table S2). After being induced with bean pyralid larvae, the number of up-regulated genes was significantly higher than that of down-regulated genes. These results indicated that most of the genes were activate and a few genes were inhibited after insect feeding.

To screen the constitutive defense genes, the highly resistant line and highly sensitive line were compared at different feeding times. The results showed that 605 DEGs were identified in Gantai-2-2 at 0 h of feeding, of which 192 DEGs were up-regulated and 413 DEGs were down-regulated (Additional file 3: Table S3), compared to Wan82–178. And at 48 h feeding, 468 DEGs were identified in Gantai-2-2, of which 202 DEGs were up-regulated and 266 DEGs were down-regulated (Additional file 4: Table S4), compared to Wan82–178.

The DEGs were further divided into three categories. The first category was the “DEGs with non-bean pyralid-induced genotype”, and there were 605 DEGs in total. This class of genes was the “DEGs identified in Gantai-2-2 compared to Wan 82-178 before bean pyralid feeding induction”, in which 52 DEGs were always up-regulated and 83 DEGs were always down-regulated at 0 h and 48 h. In addition, 9 DEGs were up-regulated at 0 h but down-regulated at 48 h, and 2 DEGs were down-regulated at 0 h but up-regulated at 48 h, whereas the other 459 DEGs displayed no changes at 48 h (Fig. 4a). The second category was the “bean pyralid-induced DEGs that appeared in both materials”. This category included 315 DEGs, which mainly displayed an up-regulated trend, with 274 DEGs were up-regulated and 31 DEGs were down-regulated in the two materials. A total of 8 DEGs were down-regulated in Gantai-2-2 but up-regulated in Wan82–178, 2 DEGs up-regulated in Gantai-2-2 but down-regulated in Wan82–178 (Fig. 4b). The third type was the “bean pyralid-induced genotype DEGs”, which consisted of a total of 1114 DEGs, of which 749 DEGs were only expressed in Gantai-2-2 and 365 DEGs were only expressed in Wan82–178.

Fig. 4
figure 4

Venn diagram of the distribution of DEGs. a HRK48/HRK0 and HSK48/HSK0. b HRK0/HSK0 and HRK48/HSK48. The circles are proportional to the number of genes identified in each treatment. The overlapping regions indicate the number of common genes. The ↑ indicate up-regulated, ↓ indicate down-regulated, ↑↓ indicate up-regulated in HRK48/HRK0 or HRK0/HSK0 but down-regulated in HSK48/HSK0 or HRK48/HSK48, ↓↑ indicate up-regulated in HSK48/HSK0 or HRK48/HSK48 but down-regulated in HRK48/HRK0 or HRK0/HSK0

Gene Ontology (GO) annotation analysis of the DEGs

To further analyze the cellular components, molecular functions and biological processes of the DEGs, GO annotation analysis was performed on all of the identified DEGs. The results showed that, under bean pyralid larvae feeding for 48 h, 572 DEGs (53.76%) of Gantai-2-2 were annotated to 42 functional groups, including 18 biological processes, 12 cellular components and 12 molecular functions compared to 0 h (Fig. 5a, Additional file 5: Table S5). Under bean pyralid larvae feeding for 48 h, 378 DEGs (55.59%) of Wan82–178 were annotated to 41 functional groups, including 18 biological processes, 11 cellular components and 12 molecular functions compared to 0 h (Fig. 5b, Additional file 6: Table S6). When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 285 DEGs (47.11%) were annotated to 39 functional groups, including 17 biological processes, 12 cellular components and 10 molecular functions (Fig. 5c, Additional file 7: Table S7). When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 240 DEGs (51.28%) were annotated to 34 functional groups, including 15 biological processes, 9 cellular components and 10 molecular functions (Fig. 5d, Additional file 8: Table S8).

Fig. 5
figure 5

GO function analysis of the DEGs. a HRK48/HRK0. b HSK48/HSK0. c HRK0/HSK0. d HRK48/HSK48

Among GO annotations of the DEGs in the above four groups, the largest common functional groups were metabolic process, cellular process, single-organism process, responses to stimuli, catalytic activities and binding. These results indicated that when the soybean was subjected to bean pyralid larvae feeding, the defense systems in the plants would immediately respond to the stimuli, appropriately increase metabolic activities in vivo and produce defense substances, such as defendant enzymes, and proteinase inhibitors, thereby enhancing the activities of various enzymes to promote defense.

Pathway analysis of the DEGs

Pathway-based analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. The results showed that 614 DEGs (57.71%) of Gantai-2-2 were assigned to 256 pathways under bean pyralid larvae feeding for 48 h, compared to 0 h. Which mainly included the metabolic pathways (263, 42.83%), biosynthesis of secondary metabolites (206, 33.55%), microbial metabolism in a diverse environment (85, 13.84%), flavonoid biosynthesis (52, 8.47%), phenylpropanoid biosynthesis (47, 7.65%) and phenylalanine metabolism (33, 5.37%) (Additional file 9: Table S9, Fig. 6a). And 380 DEGs (55.88%) of Wan82–178 were assigned to 224 pathways under bean pyralid larvae feeding for 48 h, compared to 0 h, which mainly included the metabolic pathways (147, 38.68%), biosynthesis of secondary metabolites (102, 26.84%), phenylpropanoid biosynthesis (27, 7.11%), phenylalanine metabolism (22, 5.79%) and flavonoid biosynthesis (14, 3.68%) (Additional file 10: Table S10, Fig. 6b). When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 367 DEGs (62.15%) were assigned to 208 pathways, which mainly included plant-pathogen interaction (84, 22.89%), insulin signaling pathway (14, 3.81%), phosphatidylinositol signaling system (13, 3.54%), and oocyte meiosis (13, 3.54%) (Additional file 11: Table S11, Fig. 6c). When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 209 DEGs (63.68%) were assigned to 209 pathways, which mainly included plant-pathogen interaction (42, 14.09%), aminoacyl-tRNA biosynthesis (12, 4.03%), ABC transporters (10, 3.36%) and huntington’s disease (9, 3.02%) (Additional file 12: Table S12, Fig. 6d).

Fig. 6
figure 6

Top 20 pathway entries of the DEGs.a HRK48/HRK0. b HSK48/HSK0. c HRK0/HSK0. d HRK48/HSK48

Pathway analysis showed that differentially expressed of these genes in the metabolic pathways after induced by bean pyralid might be related to the resistance, such as plant-pathogen interactions, phenylpropanoid biosynthesis, phenylalanine metabolism, flavonoid biosynthesis, peroxisome, plant hormone signal transduction, terpenoid backbone biosynthesis, and so on, and that they played a defensive role against insect stress.

R language was used to conduct a super geometric algorithm. The result showed that the DEGs belonged to the bean pyralid-induced DEGs that appeared in both materials, mainly involving the global and overview maps, biosynthesis of the other secondary metabolites, carbohydrate metabolism and amino acid metabolism (Fig. 7a). A total of 146 DEGs were identified in Gantai-2-2 compared to Wan 82–178 before and after bean pyralid feeding that were mainly involved environmental adaptation, translation, global and overview maps and signal transduction (Fig. 7b).

Fig. 7
figure 7

Pathway classification of the DEGs. a Pathway classification of “bean pyralid-induced DEGs which appeared in both materials”. b DEGs were identified in Gantai-2-2 compared to Wan 82–178 before and after bean pyralid feeding Note: The X axis represent the percent of genes (%), and the Y axis represent the metabolic process

Analysis of DEGs potentially related to anti-bean pyralid in soybean

Man pathway cluster analysis was employed to identify the pathway classification of the DEGs. The results showed that some important DEGs were categorized into ROS removal, hormone metabolism, signaling, stress, secondary metabolism and cell wall (Tables 2 and 3).

Table 2 Functional classification of DEGs
Table 3 Comparison of some DEGs of the Gantai-2-2 and Wan82–178 after bean pyralid larvae feeding

Under biotic and abiotic stresses, a great number of ROS in plants will be produced, and the cell structure is then destroyed. Plants often remove ROS in vivo through the production of a reactive oxygen scavenging agent, in order to alleviate damage to the plants caused by oxidative stress [13]. ROS removal system mainly includes enzymes, such as POD and PPO, and low molecular weight antioxidants, such as ferredoxin and thioredoxin (TRX) [14, 15]. After bean pyralid larvae feeding for 48 h, many genes related to the ROS removal system were identified that were significantly up-regulated (Tables 2 and 3). For example, there were 23 POD, 5 PPO, 9 glutathione S-transferase (GST) and 3 TRX1 identified in Gantai-2-2. Additionally, there were 20 POD, 4 PPO and 2 GST identified in Wan82–178. When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 1 GST was down-regulated. When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 1 TRX1 was up-regulated.

JA, ethylene (ET), and other plant hormone signaling pathways can be activated after pest feeding, which in turn causes a rise in the defense gene expression levels, an accumulation of defensive compounds, and an increase in the release of volatiles, finally, the plants showed resistance to the pests [16,17,18]. Our results showed that the genes related to plant hormone had changed after bean pyralid feeding, including JA, ET and auxin (Tables 2 and 3). Eight DEGs related to JA biosynthesis were all up-regulated in Gantai-2-2, including 1 lipoxygenase (LOX), 3 linoleate 9S–LOX, 1 alpha-dioxygenase (α-DOX), 1 hydroperoxide dehydratase and 2 12-oxophytodienoic acid reductase (OPDA). Additionally, 6 DEGs related to JA biosynthesis were all up-regulated in Wan 82–178, including 1 LOX, 2 linoleate 9S–LOX, 1 α-DOX, 1 hydroperoxide alpha dehydratase and 1 OPDA. When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 1 linoleate 9S–LOX and 1 OPDA were down-regulated. When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 1 LOX and 1 OPDA were down-regulated. There were 24 DEGs related to ET biosynthesis and signal transduction identified in Gantai-2-2, of which 22 DEGs were significantly up-regulated, including 7 aminocyclopropane carboxylate oxidase (ACC oxidase), 2 ethylene receptors, 3 ethylene responsive factor 1 (ERF1), 5 EREBP-like factor, and so on. Additionally, 14 DEGs related to ET biosynthesis and signal transduction identified in Wan 82–178 were all significantly up-regulated, including 5 ACC oxidase, 2 ethylene receptors, 2 EREBP-like factor, 3 ERF1, and so on. One ACC oxidase was found to be higher in Gantai-2-2 than in Wan82–178 at 0 h. In addition, 1 EREBP-like factor was found to be higher in Gantai-2-2 than in Wan82–178 at 48 h. There were 12 significantly up-regulated DEGs associated with auxin synthesis and signal transduction pathways in Gantai-2-2, including 4 IAA-amino acid hydrolase, 1 auxin responsive GH3 gene family, 3 SAUR family proteins, and so on. Six DEGs associated with auxin synthesis and signal transduction pathways identified in Wan82–178 were significantly up-regulated, including 1 IAA-amino acid hydrolase, 3 SAUR family protein, and so on. In addition, 1 SAUR family protein was found to be higher in Gantai-2-2 than in Wan82–178 at 48 h.

The results showed that many genes related to protein kinase could be induced by bean pyralid (Tables 2 and 3). After bean pyralid feeding for 48 h, 17 protein kinases were identified in Gantai-2-2, that were significantly up-regulated, including 2 protein kinase, 1 protein kinase A, 5 LRR receptor-like serine/threonine-protein kinase FLS2, 3 serine/threonine-protein kinase PBS1 (STK), 2 serine/threonine-protein kinase SRK2, 1 serine/threonine-protein kinase WNK1 and 3 serine/threonine-protein kinase CTR1. Additionally, 7 protein kinases were identified in Wan82–178, of which 6 DEGs were up-regulated, including 1 protein kinase A, 2 LRR receptor-like serine/threonine-protein kinase FLS2, 2 serine/threonine-protein kinase PBS1 (STK) and 1 serine/threonine-protein kinase WNK1. The results showed that many genes related to Ca2+ could be induced by bean pyralid (Table 3). After bean pyralid feeding for 48 h, 7 DEGs associated with Ca2+ signaling were identified in Gantai-2-2, of which 6 DEGs were up-regulated, including 1 calmodulin, 2 calcium-binding protein CML and 3 Ca2+-transporting ATPase. Additionally, 8 DEGs associated with Ca2+ signaling were identified in Wan82–178, of which 3 DEGs were up-regulated, including 1 calcium/calmodulin-dependent protein kinase, 1 calcium-binding protein CML and 1 Ca2+-transporting ATPase.

Some genes induced by abiotic stress may have been induced by insects as well, for example, the genes associated with harm and drought stresses are often induced by chewing insects [19,20,21]. Meanwhile, Our results showed that many genes related to biotic and abiotic stresses were induced by bean pyralid too. After bean pyralid feeding for 48 h, many genes related to biotic stress could be induced by bean pyralid (Tables 2 and 3). For example, 43 DEGs related to biotic stress were identified in Gantai-2-2, of which 39 DEGs were up-regulated, including 12 PR-proteins, 8 proteinase inhibitors, 9 chitinase, 1 lectin mannose-binding 2, and so on. Additionally, 22 DEGs were associated with biotic stress identified in Wan82–178, of which 20 DEGs were up-regulated, including 7 PR-proteins, 7 proteinase inhibitors, 3 chitinase, 1 lectin mannose-binding 2, and so on. In addition, 3 proteinase inhibitors were found to be higher in Gantai-2-2 than in Wan82–178 at 48 h. After bean pyralid feeding for 48 h, 23 DEGs associated with abiotic stress, such as pathogen infection, heat stress, cold stress and drought stress, were identified in Gantai-2-2, of which 17 DEGs were up-regulated. Additionally, 16 DEGs associated with abiotic stress were identified in Wan82–178, of which 9 DEGs were up-regulated (Tables 2 and 3). When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 3 DEGs associated with abiotic stress were up-regulated and 9 DEGs were down-regulated. When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 3 DEGs associated with abiotic stress were up-regulated and 3 DEGs were down-regulated.

It has widely been reported that secondary metabolism pathways, such as isoprenoid, phenylpropanoid, flavonoid, cytochrome P450 and simple phenol metabolism, were involved in plant resistance to insect feeding, or possibly functioned as direct defense compounds, antioxidants, signaling molecules or insect toxins [22,23,24]. Our results showed that many genes related to secondary metabolism pathways could be induced by bean pyralid (Tables 2 and 3). After bean pyralid feeding for 48 h, 6 DEGs related to isoprenoid biosynthesis pathway were identified in Gantai-2-2 that were all up-regulated, including 1 1-deoxy-D-xylulose-5-phosphate synthase (DXS), 2 homogentisate phytyltransferas (HPT), 1 prolycopene isomerase and 2 isoprene synthase. Additionally, 6 DEGs were identified in Wan82–178 that were all up-regulated, including 1 HPT, 2 isoprene synthase, 1 (E)-4-hydroxy-3-methylbut-2-enyl-diphosphate synthase, 1 acetyl-CoA C-acetyltransferase and 1 farnesyl diphosphate synthase. When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 1 prolycopene isomerase was down-regulated. When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 1 acetyl-CoA C-acetyltransferase was down-regulated. There were 13 DEGs associated with phenylpropanoid biosynthesis pathway were identified in Gantai-2-2 that were all up-regulated, including 3 phenylalanine ammonia-lyase (PAL), 3 caffeoyl-CoA 3-O-methyl transferase (CCOAOMT), 3 caffeic acid 3-O-methyltransferase (COMT), 1 cinnamyl-alcohol dehydrogenase (CAD), 2 trans-resveratrol di-O-methyltransferase and 1 shikimate O-hydroxycinnamoyl transferase. Additionally, 4 DEGs were identified in Wan82–178 were all up-regulated, including 2 COMT and 2 CAD. When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 1 CCOAOMT, 1 CAD and 2 shikimate O-hydroxycinnamoyl transferase were down-regulated, and 2 shikimate O-hydroxycinnamoyl transferase were up-regulated. When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 1 COMT, 1 CAD and 1 shikimate O-hydroxycinnamoyl transferase were down-regulated, and 2 shikimate O-hydroxycinnamoyl transferase was up-regulated. There were 22 DEGs associated with flavonoid biosynthesis pathway identified in Gantai-2-2 that were all up-regulated, including 4 chalcone isomerase (CHI), 8 chalcone synthase (CHS), 3 dihydroflavonol4-reductase (DFR) and 7 leucoanthocyanidin reductase (LAR). Additionally, 4 DEGs were all up-regulated identified in Wan82–178, including 1 CHI, 1 DFR and 2 LAR. When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 1 CHS was up-regulated. After bean pyralid feeding for 48 h, 26 DEGs related to cytochrome P450 (CPY) were identified in Gantai-2-2, of which 25 DEGs were up-regulated. Additionally, 8 DEGs were identified in Wan82–178, of which 7 DEGs were up-regulated. When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 1 CPY was up-regulated and 1 CPY was down-regulated. When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 1 CPY was up-regulated. After bean pyralid feeding for 48 h, 4 and 1 L-ascorbate oxidase (L-AO) were identified in Gantai-2-2 and Wan82–178, respectively, which were all up-regulated.

Bean pyralid-induced the transcription factor genes

Transcription factors, known as trans-acting factors, are a type of DNA binding protein that regulates the transcription of the target genes by combining with cis-acting elements in the gene promoter [25]. Research demonstrated that plant transcription factors, including NAC, MYB, WRKY, were involved in the plant defense responses [26]. Our results showed that 25 transcription factors were identified in Gantai-2-2, including 4 NAC, 7 AP2-EREBP, 3 MYB, 2 WRKY, 1 bHLH, 2 C3H, and so on, of which 21 were up-regulated and 4 were down-regulated, after bean pyralid larvae feeding for 48 h (Table 4). Fifteen transcription factors were identified in Wan82–178, including 3 NAC, 4 AP2-EREBP, 1 MYB, 1 WRKY, 2 bHLH, 3 C3H, and so on, which were all up-regulated, after bean pyralid larvae feeding for 48 h. When comparing Gantai-2-2 with Wan82–178 at 0 h feeding, 7 transcription factors were identified, of which 1 AP2-EREBP, 1 C2C2-YABBY and 1 FAR1 were up-regulated, and 1 AP2-EREBP, 2 MYB and 1 LIM were down-regulated. When comparing Gantai-2-2 with Wan82–178 at 48 h feeding, 5 transcription factors were identified, of which 1 FAR1 was up-regulated, and 1 AP2-EREBP, 2 MYB and 1 C3H were down-regulated (Table 4).

Table 4 Transcription factors (TF) related to the insect resistance identified in different alignment schemes

Analysis of DEGs by qRT-PCR

The quantitative real time-PCR (qRT-PCR) technology was used to verify the credibility of the RNA-Seq. Therefore, we used the qRT-PCR technology to verify 17 DEGs identified by RNA-Seq. It was found that the qRT-PCR expression patterns of 17 DEGs were all consistent with the RNA-Seq results (Fig. 8), which indicated that the RNA-Seq results were reliable in the present study.

Fig. 8
figure 8

DEGs confirmed by qRT-PCR using the same sample as that in RNA-Seq. X-axis represented gene name, the blue column represented qRT-PCR results in HRK48/HRK0, the red column represented RNA-Seq results in HRK48/HRK0, the green column represented qRT-PCR results in HSK48/HSK0, and the purple column represented RNA-Seq results in HSK48/HSK0; Y-axis represented the relative level of gene expression

Discussion

Influence of the genes related to ROS removal on bean pyralid response

As an electron acceptor of H2O2, POD oxidizes various materials in the process of secondary metabolites and can produce phenoxy and oxygen free radicals by interacting with some phenolic compounds, which may also directly interfere with the feeding of insects or reduce the nutrition level of leaves, thereby lowering the edibility of the plants [27]. Meanwhile, POD has been found to have significant insecticidal effects on Lepidoptera and Coleoptera [28]. Previous studies have shown that POD could be induced by insects in wheat [29], sorghum [30], cucumber [31] and rice [32, 33]. PPO is widely distributed in plants and can potentially catalyze the oxidation of polyphenols into a keto metal enzyme, it can directly oxidize acid into quinone using O2 as an oxidation substrate, then produce amino acids and proteins that are difficult to digest and toxic to herbivores. Meanwhile, as an anti-nutritional protein, PPO can produce protein affinity, and therefore reduce the edibility of plants, which plays an important role in the defense of insect feeding [34, 35]. Previous research studies have shown that the activity of PPO increased after tomatoes were eaten by Spodoptera exigua or were treated with oral secretions of Spodoptera exigua [36, 37]. GST is encoded by a large and diverse family of genes, that plays important roles in the responses of plants to oxidative damage induced by various environmental conditions, especially in the resistance to resist the toxic effects of ROS on cells [38,39,40]. In oxygen stress reactions, TRX would pass restoring forces to the reductase, which could potentially remove the lipid peroxide or repair the oxidized protein, resulting in an alleviation of oxygen stress [41]. TRXh5 and TRXh8 in Arabidopsis thaliana were strongly induced under biotic or abiotic stress [42, 43]. After bean pyralid larvae feeding for 48 h, we found that POD, PPO and GST were all significantly up-regulated in Gantai-2-2 and Wan82–178. These results suggested that these genes were involved in the defense reactions of soybean to bean pyralid. TRX1 was up-regulated in Gantai-2-2, which suggested that TRX1 involved in defense responses to bean pyralid in the highly resistant material only.

Influence of plant hormones on bean pyralid response

JA signaling pathway is implicated in insect-induced responses in plant; it was involved in the signal transduction of insect defense, regulated the expression of plant downstream defense genes, and significantly induced responses of defense systems in plant, thereby effectively reducing pests [44]. After bean pyralid feeding, 3 key genes involved in JA biosynthetic pathway were identified in Gantai-2-2 and Wan82–178, namely, LOX, α-DOX and OPDA, which were significantly up-regulated. LOX is the key enzyme in the synthesis of the JA [45], and plays an important signal factor in plant induction defense pathways [46]. Adversity stresses, such as insects and diseases, could induce single or multiple LOX genes in plants [47, 48]. Additionally, α-DOX could catalyze the oxidation of fatty acids and produce 2-hydrogen peroxide fatty acids, and it is also an enzyme related to plant stress resistance in JA biosynthetic pathway [32]. For example, α-DOX can be induced by insects in rice [32, 49]. OPDA is the intermediate product of JA biosynthesis pathway with biological activities that can regulate plant defense genes. For example, OPDA has been found to induce defense genes in Arabidopsis thaliana which resist attacks from Bradysia odoriphaga [50]. As a plant endogenous hormone, ET regulates the defense responses of plants to pests and diseases [51, 52]. Our results found that 3 types of genes related to ET biosynthesis and signal transduction pathway. For example, most of ACC oxidase, ethylene receptor and ERF1 genes were up-regulated after bean pyralid feeding. ACC oxidase is one of the key enzymes in ethylene biosynthesis pathway [52]. Ethylene receptor is an upstream component of ethylene signal transduction pathway, that plays an important role in plant growth and responses to adversity stresses [53]. ERF is a specific plant transcription factor that plays an important role in plant cell growth, hormone regulation, disease resistance and abiotic stresses [54,55,56]. Auxin can directly act on cell membranes or intracellular components, and affects some cellular responses. In addition, it can also indirectly regulate the expression of genes, and has a direct impact on plant stress [57, 58]. After bean pyralid feeding, some genes related to auxin synthesis and signal transduction pathways were identified, such as IAA-amino acid hydrolase, auxin responsive GH3 (Gretchen Hagen3) gene family and SAUR (small auxin up RNA) family protein. Early/primary auxin response genes were composed of three major gene families: Aux/IAA (auxin/indoleacetic acid), GH3 and SAUR transcription factor family [59]. GH3 gene family is a typical auxin-responsive gene family and can promote amino acids to be combined with IAA, JA and SA, then changed the concentration of biologically active forms within the cells. This aided in regulating plant growth, development and defense responses [60]. SAUR gene family is known to be the specific and largest family among the auxin response factors in plants and is related to environmental stimuli [61, 62]. These results indicated that JA, ET and auxin may be involved in elevating the basal resistance of soybean to herbivory.

Influence of the intracellular signal transduction pathway on bean pyralid response

Protein kinase is involved in the signal transduction pathway of plants and plays an important role in the defense responses of plants [63]. Protein kinases could be induced by bean pyralid, including protein kinase, protein kinase A, LRR receptor-like serine/threonine-protein kinase FLS2, serine/threonine-protein kinase PBS1 (STK), serine/threonine-protein kinase SRK, serine/threonine-protein kinase WNK1, and serine/threonine-protein kinase CTR1. LRR receptor-like kinase is the largest family of receptor protein kinase in plants; it plays an important role in the regulation of plant growth, development, biotic and abiotic stresses [64]. Serine/threonine-protein kinase (STK) is an important signal molecule; when plant suffers stimulation, such as pest feeding, salinity, drought stress, trauma, cytokines or hormones, STK is rapidly activated in the serine and threonine residue phosphorylation sites and is further activated in the downstream signal molecules through phosphorylation cascades, which activate specific signaling pathways. It is eventually transmitted by an outside signal to the nucleus, and activates or inhibits specific genes [64,65,66]. These results indicated that protein kinases played an important role in the defense against bean pyralid.

When the plants were stimulated by external environmental and received the signal, the receptor activated the calcium channels on the membrane through a series of phosphorylation reactions to cause calcium ions to be released from the calcium base into the cytoplasm. This response led to an increase in the concentration of Ca2+ in the cytoplasm, which then activated the plant defense response [67, 68]. Previous studies have shown that Ca2+ could be induced by Spodoptera littoralis in Phaseolus lunatus [69] and Ginkgo biloba [70]. After bean pyralid feeding, most of the DEGs related to Ca2+ signaling identified in Gantai-2-2 were up-regulated, but in Wan82–178 were down-regulated. It was suggested that Ca2+ signaling was involved in the defense responses to bean pyralid in the highly resistant material. After soybean was stimulated by bean pyralid, the concentrations of Ca2+ in the cytoplasm changed; therefore, calcium was sent to transfer the stimulation signal. However, Ca2+ signaling may have played a negative role in the regulation of defense responses in the highly susceptible material.

Influence of the genes related to biotic stress on bean pyralid response

When plants suffer from pest stress, a variety of defense proteins which inhibit the insect from producing digestive enzymes are produced to resist insect pests, including PR, proteinase inhibitors, chitinase and lectin [71, 72] and thereby destroy the insect’s normal digestive absorption functions, then disrupt their nutrient uptake and utilization, which could ultimately lead to malnutrition and inhibited growth of the insects [73].

PR proteins are generated and accumulated after pathogen invasions or abiotic stress and are inducible components in the self defense mechanisms of plants [74,75,76]. Uknes et al. found that the contents of PR-1, PR-2 and PR-5 were increased, when the Arabidopsis thaliana were treated with TCV and INA, and its disease resistance was increased [77]. After bean pyralid feeding, most of the DEGs related to PR were induced. Our results speculated that the PR proteins might be involved in the defense responses of soybean to bean pyralid.

Proteinase inhibitors can inhibit the activity of chymotrypsin and trypsin in the digestive tract of the Lepidoptera and Diptera, to inactivate the related digestive enzyme in intestines, and thereby interfering with the insects’ normal feeding and digestion. Therefore, the insects cannot attain enough nutrition, which in turn affects their growth and development, then leads to insect death [78, 79]. In other cases, it caused excessive production of proteases in insects, which resulted in a deficit of essential amino acids, then caused insect death [80]. Our results showed that the proteinase inhibitors could be induced by bean pyralid, so we assumed that the proteinase inhibitors were involved in the defense responses.

Chitinase could destroy the peritrophic membranes of insects [81]. When insect eats plants with sustained chitinase expression, its digestive tract becomes damaged and its epidermis cannot form normally [82]. Chitinase could be induced by aphid and spider mites in sorghum [83] and tomato [84], respectively. Our results showed that chitinase could be induced by bean pyralid as well. It was assumed that after the chitinase was absorbed into the insects’ bodies, chitin was hydrolyzed, which inhibited the growth and development of insects and achieved the purpose of defense for the plants [85].

Lectin has at least one non-catalytic domain-specific reversible binding to monosaccharides or oligosaccharides, and after being absorbed into the gut of an insect, lectin induces local or systematic toxic effects, which bond to glycoproteins on the peritrophic membrane and damage the structure of the peritrophic membrane. This response causes insect growth arrest, antifeedant and even death [86, 87]. Previous results have shown that lectin displays insecticidal activities against a variety of insects. For example, after Heliothis virescens [88], Myzus persicae [89] and Lacanobia oleracea [90] fed on lectin-transgenic plants, their survival rates and fecundity were greatly reduced. Our results showed that the lectin, mannose-binding 2 could be induced by bean pyralid. It was assumed that the lectin that was released from damaged insect cells of the soybean combined with the chitin of the peritrophic membranes, sugar compounds of the digestive tract epithelial cells, and the glycosylation digestive enzymes, after bean pyralid feeding. This affected the normal absorption of nutritional materials. It also induced disease in the digestive tracts and promoted gastrointestinal bacterial reproduction, thereby causing insect growth inhibition or death, which achieved the defensive purpose of killing the pests [91].

Influence of the genes related to the secondary metabolism on bean pyralid response

Isoprenoid biosynthesis pathway is an important secondary metabolism route that widely exists in plants. It has the ability to synthetize isoprenoid compounds, and plays an important role in plant resistance reactions, in addition to participating in the growth and development of plants [92]. The genes related to isoprenoid biosynthesis pathway could be induced by bean pyralid, which speculated that isoprenoid biosynthesis pathway was involved in the defense reactions of soybean to bean pyralid.

Phenylpropanoid and its derivatives play an important role in plant development and stress responses. Environmental stresses can promote plant carbon synthesis, which changes into benzene propane synthesis, causing the accumulations of large amounts of substances related to plant stress resistance, such as lignin, flavone, flavonol and alkaloids. These substances have been found to improve the resistance of plants to various biotic and abiotic stresses [93]. The genes related to phenylpropanoid biosynthesis pathway, including PAL, CCOAOMT, COMT and CAD, could be induced by bean pyralid. PAL is not only the enzyme that catalyzes the first step of phenylpropanoid biosynthesis pathway but also the key and rate limiting enzyme for the pathway. Various metabolite substances that are catalyzed by PAL, such as lignin, flavonoids and phytoalexin, were related to plant insect resistance [94]. PAL can be induced by many types of biotic and abiotic stresses, including pathogen infections, insect feeding, mechanical injuries and exogenous plant hormones [95,96,97]. CCOAOMT is a type of important methyltransferase in plant lignin biosynthesis [98, 99]. COMT catalyzes the reaction that converts caffeic acid into ferulic acid, which is an important step in the biosynthesis of lignin monomers [100]. CAD is the last step in lignin synthesis metabolic pathway and plays a key role in lignin biosynthesis [101, 102]. PAL, CCOAOMT, COMT and CAD were related to lignin, which suggested that soybean might have enhanced its cell mechanical strength through the synthesis of lignin, which reinforced the cell walls, strengthened the cell to blocked insect feeding, and defended the insect pests by nutrient limitation [103, 104].

Flavonoid plays an important role in insect-induced, disease, and other stress responses in plants [105]. The genes related to flavonoid biosynthesis pathway, including CHI, CHS, DFR and LAR, were found to be induced by bean pyralid. CHS and CHI are two key enzymes in flavonoid biosynthesis pathway; CHS catalyzes the condensation reaction of 4-coumaric acid-CoA and acyl-CoA, which causes chalcone to be formed, and flavonoid is then formed under the action of CHI, a necessary enzyme that syntheses flavanone, flavone, flavonol and anthocyanin substances [106, 107]. The amount of flavonoid metabolites was directly affected by the expression of CHS and CHI, and the loss of expression ability or loss of enzyme function [108]. DFR and LAR are key enzymes in flavonoid synthesis, and during the process of plant flavonoid biosynthesis, anthocyanins, proanthocyanidins and other common synthesis precursors are leucoanthocyanins. DFR was catalyzed by flavanonols to generate leucoanthocyanidin, LAR converted leucoanthocyanidin into 2,3-trans-flavan-3-ols [109]. Previous studies have shown that over-expressions of DFR could improve the content and antioxidant capacities of flavonoid [110]. Our results suggested that flavonoid biosynthesis pathway played an important role in the defense against bean pyralid and that flavonoid accumulation is an important characteristic of responses to abiotic stress in plants.

Cytochrome P450 (CYP) is a type of pheme containing oxidoreductases. It can catalyze some substances with a defensive function, such as sterol, isoflavone, alkaloid, terpenoid, and so on [111, 112]. CYP plays an important role in the defense of organisms against pests and abiotic stress [113]. CYP genes can be induced by aphid in soybean [114]. Our results showed that some genes related to CYP can also be induced by bean pyralid. Previous studies have shown that cyanogentic glycoside which was catalyzed and synthesized by the CYP79A and CYP71E1 genes in sorghum, was harmful to pests [115]. Therefore, it was speculated that soybean would use the CYP family to mitigate the threat of insect infestation.

AO is a copper-containing oxidase family that is wide spread in plants. It is a type of simple phenol and can catalyze the ascorbic acid (AA) of protoplasts in vitro and can oxygenate AA to generate monodehydroascorbate (MDHA), which regulates the overall oxidation state of the AA library of plant protoplasts in vitro [116, 117]. When activity of AO increased, the ratio of oxidation and reduction in extracellular AA was also increased. There was a significant redox gradient on both sides of the plasma membrane, which modulated plant defense responses against insect attack and other stressed [118,119,120]. Our results showed that some genes related to L-AO can also be induced by bean pyralid. Therefore, L-AO was involved in the response reactions of soybean to bean pyralid.

Transcription factor on bean pyralid response

WRKY is one of the largest families of transcription factors in plants and plays an important role in the regulation of plant growth and development as well as biotic and abiotic stresses [121,122,123]. Previous studies have shown that WRKY could be induced by insect stress [124]. For example, one WRKY23 gene was induced by Heterodera schachtii in Arabidopsis thaliana [125] and 20 WRKY transcription factors were induced by cotton boll weevil in Gossypium hirsutum [126]. NAC is known to be a second major transcription factor in plants and is also a specific transcription factor. NAC can activate downstream genes when plants are exposed to biotic or abiotic stress, so it is involved in the stress responses of plants [127,128,129]. Some NAC genes were induced by insects in rape [130], sugarcane [131] and Gossypium hirsutum [126]. MYB is one of the central regulators of development, metabolism, and response to abiotic and biotic stresses [132]. One of the known mechanisms is that MYB transcription factor contribute to the accumulation of flavonoids to protect plants from radiation injury and enhance their resistance to cold and insect pests in plants [133]. MYB could be induced by small cabbage white caterpillars [134] and cotton boll weevil [126]. AP2/EREBP is a type of specific transcription factor which is bound with ethylene responsive elements. It is able to combine with the GCC box of the promoter region of some resistance related genes, and carry out the function of activating the expression of these genes [135]. AP2/EREBP can potentially regulate the molecular responses of plants to hormones, pathogens, low temperatures, drought and high salt, which could improve the tolerance of crops to stress [136,137,138]. For example, the over expression of Tsi1 gene encoding EREBP/AP2 transcription factor increased the resistance to pathogens in tobacco, as well as the resistance ability to osmotic stress [138]. Our results showed that the transcription factor, such as WRKY, NAC, MYB and AP2/EREBP, could be induced by bean pyralid. It was confirmed that these transcription factors were involved in the defense reactions to bean pyralid.

Conclusions

To explore the defense mechanisms of soybean resistance to bean pyralid, the comparative transcriptome sequencing was completed between the leaves infested with bean pyralid larvae and no worm of soybean (Gantai-2-2 and Wan82–178) on the Illumina HiSeq™ 2000 platform. The results showed that there were 1064 and 680 DEGs were identified in the Gantai-2-2 and Wan82–178 after bean pyralid larvae feeding for 48 h, respectively, compared to 0 h. When comparing Gantai-2-2 with Wan82–178, 605 DEGs were identified at 0 h feeding, and 468 DEGs were identified at 48 h feeding. The DEGs were divided into three categories, including the DEGs with non-bean pyralid-induced genotype, bean pyralid-induced DEGs that appeared in both materials and bean pyralid-induced genotype DEGs. According to GO and KEGG functional and metabolic pathway analysis combined with the previously reported literatures, we concluded that the response to bean pyralid feeding might be related to the disturbed functions and metabolism pathways of some key DEGs, such as DEGs involved in the ROS removal system, mainly POD, PPO, GST and TRX; phyto-hormone signaling pathways, mainly JA, ET and auxin pathway; intracellular signal transduction pathways, including plant protein kinases and Ca2+ signaling; secondary metabolism, such as isoprenoid biosynthesis pathway, phenylpropanoid biosynthesis pathway, flavonoid biosynthesis pathway, CYP and L-AO; transcription factors, such as WRKY, NAC, MYB and AP2/EREBP. Meanwhile, bean pyralid activated a large number of genes related to biotic and abiotic stresses. These results will help to elucidate the molecular mechanism of response to bean pyralid feeding in soybean, and provide a valuable resource of soybean defense genes that will benefit other studies in this field. Future research will focus on the cloning and transgenic function validation of possible candidate genes associated with the anti-bean pyralid soybean.

Methods

Plant materials

The tested materials Gantai-2-2 (highly resistant line) [3] and Wan82–178 (highly susceptible line) [3] were planted inside insect net rooms in experimental fields at the Guangxi Academy of Agricultural Sciences. The plants were sown in three rows, with 10 strains per row. During the entire growth period of soybeans, pesticides and fertilizers were not used. When the plants grew to a level of 10 compound leaves, the seventh compound leaves on the left side were collected before the infestation. There was a total of five plants for each sample. Then, each of the samples was simultaneously artificially infested with 5 four-year-old bean pyralids on the right side of the seventh leaves, counted downward. There were two biological repetitions. In each sample, five leaves were mixed, and they were frozen with liquid nitrogen and stored at −80 °C for further use.

Total RNA extraction

The total RNA (5 μg) from the leaves (100 mg per sample) of Gantai-2-2 and Wan82–178 by using TRIzol kit (Invitrogen, Carlsbad, CA, USA) according to the manufacture’s protocol. Briefly 1.3 ml of TRIzol kit was added into 100 mg of leaf sample, the sample was homogenized using power homogenizer and centrifuged at 12,000×g for 10 min at 4 °C. After the fatty layer was removed and discarded, the cleared supernatant was transferred into a new tube and mixed with 0.2 ml of chloroform. The sample tube was shaken for 15 s, followed by an incubation for 3–5 min at room temperature. Next, the sample was centrifuged 12,000×g for 10 min at 4 °C and the aqueous phase was moved into a new tube for RNA precipitation. For precipitating RNA from each sample, 10 μg of RNase-free glycogen was added to the aqueous phase as a carrier, followed by 0.5 ml of 100% isopropanol, then samples were placed at −20 °C for 1 h and centrifuged at 13,600×g for 20 min at 4 °C and discarded the supernatant. To wash the RNA pellet, we added 1.0 ml of 75% ethanol into the tube, vortexed the tube gently, centrifuged the tube 13,600×g for 3 min at 4 °C and discarded the wash. The RNA pellet was air-dried, suspended in 50 μl Nuclease-free water, incubated at 55–60 °C for 10 min. RNA concentration, purity and integrity were determined using a NanoDrop2000 (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA).

cDNA library construction and transcriptome deep sequencing

Equal amount of total RNA (5.0 μg) was used for cDNA library construction using TruSeq™ RNA Sample Prep-aration Kit v2 (Illumina) (Illumina, SanDiego, CA, USA) and the cDNA library was sequenced on an Illumina HiSeq™ 2000 platform (Hiseq2000 Truseq SBS Kit v3-HS (200 cycles), Illumina) following the manufacturers’ protocols. Briefly, Dynal Oligo (dT) beads (Invitrogen) were obtained poly(A) mRNAs. Then mRNAs were chemically fragmented into ~200 nt fragments. First-strand cDNA was generated by using reverse transcriptase and random primers, the second strand cDNA synthesis using DNA Polymerase I (Invitrogen) and RNase H (Invitrogen) treatment. The cDNA fragments were end repaired by using End Repair Mix (Illumina) reagent, Finally, purified and enriched to create the final cDNA library. Eight cDNA libraries were sequenced by using pair-end (2 × 100 bp) sequencing technology with an Illumina HiSeq™ 2000 sequencer (The RNA-Seq test and its results were analyzed by the BGI Tech Solutions Co., Ltd. (BGI Tech) (BGI, Shenzhen, PR China).

Raw read filtering and mapping to the reference genome and gene sequences

SOAPfuse software was used to filter noise for the original sequencing reads [139]. Raw reads of eight libraries were performed by removing adapter sequences and low-quality reads, higher N rate sequences, and excessively short sequences. After passing the QC process of the alignment, the results were used for further analysis. The remaining high-quality reads were submitted for mapping analysis against the soybean reference genome (ftp://ftp.jgi-psf.org/pub/compgen/phytozome/v9.0/early_release/Gmax_275_Wm82.a2.v1/, version Glyma 2.0, 975 Mb), using BWA software and Bowtie software and allowing two base mismatches. We assembled the transcripts with reads using Cufflinks (http://cufflinks.cbcb.umd.edu/) [140].

DEG analysis

Genes and isoforms expression level are quantified by a software package: RSEMV1.2.12 (RNASeq by Expectation Maximization) [141]. RSEM computes maximum likelihood abundance estimates using the Expectation-Maximization (EM) algorithm for its statistical model, including the modeling of paired-end (PE) and variable-length reads, fragment length distributions, and quality scores, to determine which transcripts are isoforms of the same gene. The expression quantity of each gene (fragments per kilobase of exon model per million mapped fragments, FPKM) was used for the calculated expression level, the formula for which is as follows:

$$ \mathrm{FPKM}\ \left(\mathrm{A}\right)=\frac{10^6C}{N/{10}^3} $$

where FPKM (A) is the expression of gene A, C is the number of fragments that are uniquely aligned to gene A, N is the total number of fragments that are uniquely aligned to all genes, and L is the number of bases on gene A. The FPKM method can eliminate the influence of different gene lengths and sequencing discrepancies in the calculation of gene expression. Therefore, the calculated gene expression can be directly used for comparing the difference in gene expression among samples. According to the Noiseq method, an absolute value of |Log2FC (Fold Change)| ≥ 1 and diverge probability ≥0.8 were used as the threshold to screen DEGs. In the DESeq2 method, the absolute value of |Log2FC (Fold Change)| ≥ 1 and p-adj ≤ 0.001 were used as the threshold to judge the significance of the gene expression difference. For the EdgeR method, the absolute value of |Log2FC (Fold Change)| ≥ 1 and FDR (False Discovery Rate) ≤ 0.001 were used as the threshold to judge the significance of the gene expression difference.

Bioinformatics analysis

GO (GO, http://www.geneontology.org/) and functional enrichment analysis were conducted on all DEGs using TermFinder software (http://www.yeastgenome.org/help/analyze/go-term-finder). Then, all DEGs were mapped to a pathway in the KEGG database (http://www.genome.jp/kegg/pathway.html) using Blast_v2.2.26 software. A p-value was used as ≤0.05 was used as the threshold to judge the significance of the GO and KEGG pathway enrichment analyses.

The transcript version number of Wm82.a2.v1, which was obtained from the sequencing analysis was deleted from the website (http://www.soybaen.org/correspondence/) and only the names of the genes were retained. These names were converted into the differential gene transcription version number of Wm82.a2.v1. MapMan software (http://www.gabipd.de/projects/MapMan/) was used to assign the pathways and functional classifications to the DEGs. The latest Osa_MSU_v7 mapping file and pathways files were downloaded from the website.

Quantitative real time-PCR (qRT-PCR) analysis

qRT-PCR analysis was used to verify the RNA-Seq results. Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA) was used to design primers for qRT-PCR experiment (Additional file 13: Table S13). PrimeScript RT Reagent kits with gDNA Eraser (Takara, Dalian, China) were used to reverse transcriptase and synthesize the cDNA. The qRT-PCR reaction mixture (25 μl) contained 12.5 μl SybrGreen qPCR Master Mix (2 × concentration, Ruian Biotechnologies, Shanghai, China), 0.5 μl reverse and forward primers (10 μM), 2 μl cDNA and 9.5 μl ddH2O. Then, the qRT-PCR reaction was performed on an ABI7500FAST Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) as follows: 2 min at 95 °C, 40 cycles of heating at 95 °C for 10 s and annealing at 60 °C for 40 s. The β-actin gene was used as the internal control to calculate the relative expression. 2-ΔΔcq method was used to calculate the differential expression of a gene in two samples. Three technological and three biological reactions (n = 3 × 3) were performed for every gene in each sample.

Abbreviations

AA:

Ascorbic acid

ACC oxidase:

Aminocyclopropane carboxylate oxidase

CAD:

Cinnamyl-alcohol dehydrogenase

CCOAOMT:

Caffeoyl-coA 3-O-methyl transferase

CHI:

Chalcone isomerase

CHS:

Chalcone synthase

COMT:

Caffeic acid 3-O-methyltransferase

CYP:

Cytochrome P450

DEGs:

Differentially expressed genes

DFR:

Dihydroflavonol 4-reductase

DXS:

1-Deoxy-D-xylulose-5-phosphate synthas

ERF1:

Ethylene responsive factor 1

ET:

Ethylene

GO:

Gene Ontology

GST:

Glutathione S-transferase

HPT:

Homogenitisate phytyltransferas

JA:

Jasmonic acid

KEGG:

Kyoto Encyclopedia of Genes and Genomes

L-AO:

L-ascorbate oxidase

LAR:

Leucoanthocyanidin reductase

LOX:

Lipoxygenase

MDHA:

Monodehy-droaseorbate

OPDA:

12-oxophytodienoic acid reductase

PAL:

Phenylalanine ammonia-lyase

POD:

Peroxidase

PPO:

Polyphenol oxidase

qRT-PCR:

Quantitative real time-PCR

ROS:

Reactive oxygen species

TRX:

Thioredoxin

α-DOX:

alpha-dioxygenase

References

  1. Wlson RF. Chapter 1 soybean; market driven research needs. In: Stacey G, editor. Genetics and genomics of soybean. New York: Springer; 2008. p. 3–15.

    Chapter  Google Scholar 

  2. Editorial committee of plate of Chinese diseases and insects on crop. Plate of Chinese diseases and insects on crop, fifth fascicule, diseases and insects on oil crop (first). Beijing: Agricultural press; 1982. p. 136–7.

    Google Scholar 

  3. Cui ZL, Gai JY, Ji DF, Ren ZJ. A study on leaf-feeding insect species on soybeans in Nanjing area. Soybean Sci. 1997;16(1):12–20.

    Google Scholar 

  4. Sun ZD, Yang SZ, Chen HZ, Li CY, Long LP. Identification of soybean resistance to bean pyralid (Lamprosema indicate Fabricicus) and oviposition preference of bean pyralid on soybean varieties. Chin J Oil Crop Sci. 2005;27(4):69–71.

    Google Scholar 

  5. Sun ZD, Chen HZ, Wei DW. A study on leaf-feeding insect species on soybeans in Nanning. Guangxi Agric Sci. 2001;2:104–6.

    Google Scholar 

  6. Long LP, Yang SZ, Chen HZ, Qin JL, Li CY, Sun ZD. Effects of different genotypes of soybean varieties on the experimental population of bean pyralid (Lamprosema indicata Fabricius). Chin J oil. Crop Sci. 2004;26:67–70.

    Google Scholar 

  7. Xing GN, Tan LM, Liu ZXN, Yue H, Zhang HZ, Shi HF, et al. Morphological variation of pubescence on leaf blade and petiole and their correlation with resistance to bean pyralid (Lamprosema indicate Fabricius) in soybean landraces. Soybean Sci. 2012;31(5):691–6.

    Google Scholar 

  8. Xing GN, Zhao TJ, Gai JY. Inheritance of resistance to Lamprosema indicata Fabricius in soybean. Acta Agron Sin. 2008;34(1):8–16.

    Article  Google Scholar 

  9. Li GJ, Cheng LG, Zhang GZ, He XH, Zhi HJ, Zhang YM. Mixed major-gene plus polygenes inheritanceanalysis for resistance in soybean to bean pyralid (Lamprosema indicate Fabricius). Soybean Sci. 2008;27(1):33–6. 41

    Google Scholar 

  10. Xing GN, Zhou B, Wang YF, Zhao TJ, DY Y, Chen SY, et al. Genetic components and major QTL confer resistance to bean pyralid (Lamprosema indicate Fabricius) under multiple environments in four RIL populations of soybean. Theor Appl Genet. 2012;125(5):859–75.

    Article  PubMed  Google Scholar 

  11. Zeng WY, Cai ZY, Zhang ZP, Chen HZ, Yang SZ, Tang XM, et al. Studies on the physiological and biochemical of bean pyralid (Lamprosema indicata Fabricius) resistance in different soybean varieties. J Southern Agric. 2015;46(12):2112–6.

    Google Scholar 

  12. Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, et al. Genome sequence of the palaeopolyploid soybean. Nature. 2010;463(7278):178–83.

    Article  CAS  PubMed  Google Scholar 

  13. Foyer CH, Shigeoka S. Understanding oxidative stress and antioxidant functions to enhance photosynthesis. Plant Physiol. 2011;155(1):93–100.

    Article  CAS  PubMed  Google Scholar 

  14. Slesak I, Libik M, Karpinska B, Karpinski S, Miszalski Z. The role of hydrogen peroxide in regulation of plant metabolism and cellular signaling in response to environmental stresses. Acta Biochim Pol. 2007;54(1):39–50.

    CAS  PubMed  Google Scholar 

  15. Quan LJ, Zhang B, Shi WW, Li HY. Hydrogen peroxide in plants: a versatile molecule of the reactive oxygen species network. J lntegr. Plant Biol. 2008;50(1):2–18.

    CAS  Google Scholar 

  16. Ellis C, Turner JG. The Arabidopsis mutant cev1 has constitutively active jasmonate and ethylene signal pathways and enhanced resistance to pathogens. Plant Cell. 2001;13(5):1025–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Farmer EE, Ryan CA. Interplant communication: airborne methyl jasmonate induces synthesis of proteinase inhibitors in plant leaves. Proce Nati Acad Sci USA. 1990;87(19):713–6.

    Google Scholar 

  18. De Vos M, Van Oosten VR, Van Poecke RM, Van Pelt JA, Pozo MJ, Mueller MJ, et al. Signal signature and transcriptome changes of Arabidopsis during pathogen and insect attack. Mol Plant-Microbe Interact. 2005;18(9):923–37.

    Article  CAS  PubMed  Google Scholar 

  19. Brodbeck BV, Lii RFM, Andersen PC. Physiological and behavioral adaptations of three species of leafhoppers in response to the dilute nutrient content of xylem fluid. J Insect Physiol. 1993;39(1):73–81.

    Article  CAS  Google Scholar 

  20. Reymond P, Weber H, Damond M, Fanner EE. Differential gene expression in response to mechanical wounding and insect feeding in Arabidopsis. Plant Cell. 2000;12(5):707–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Mozoruk J, Hunnicutt LE, Cave RD, Hunter WB, Bausher MG. Profiling transcriptional changes Citrus sinensis (L.) Osbeck challenged by herbivory from the xylem-feeding leafhopper Homalodisca coagulate (say) by cDNA macroarray analysis. Plant Sci. 2006;170(6):1068–80.

    Article  CAS  Google Scholar 

  22. Schuler MA. The role of cytochrome P450 monooxygenases in plant-insect interactions. Plant Physiol. 1996;112(4):1411–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Treutter D. Significance of flavonoids in plant resistance: a review. Environ Chem Lett. 2006;4(3):147–57.

    Article  CAS  Google Scholar 

  24. Korkina LG. Phenylpropanoids as naturally occurring antioxidants: from plant defense to human health. Cell Mol Biol. 2007;53(1):15–25.

    CAS  PubMed  Google Scholar 

  25. Zhu JH, Verslues PE, Zheng XW, Lee BH, Zhan XQ, Manabe Y, et al. HOS10 encodes an R2R3-type MYB transcription factor essential for cold acclimation in plants. Proc Nati Acad Sci USA. 2010;102(28):9966–71.

    Article  CAS  Google Scholar 

  26. Rushton PJ, Somssich IE. Transcriptional control of plant genes to responsive to pathogens. Curr Opin Plant Biol. 1998;1(4):311–5.

    Article  CAS  PubMed  Google Scholar 

  27. Felton GW, Donato K, Del Vecchio RJ, Duffey SS. Activation of plant foliar oxidases by insect feeding reduces nutritive quality of foliage for herbivores. J Chem Ecol. 1989;15(12):2667–94.

    Article  CAS  PubMed  Google Scholar 

  28. Estruch JJ, Carozzi NB, Desai N, Duck NB, Warren GW, Koziel MG. Transgenic plants: an emerging approach to pest control. Nat Biotechnol. 1997;15(2):137–41.

    Article  CAS  PubMed  Google Scholar 

  29. Leszczynski B. Changes in phenols content and metabolism in leaves of susceptible and resistant wheat cultivars infested by Rhopalosiphum padi (L.) (Hom.:Aphididae). J Appl Entomol. 1985;100(4):343–8.

    CAS  Google Scholar 

  30. Huang YH. Phloem feeding regulates the plant defense pathways responding to both aphid infestation and pathogen infection, Biotechnology and sustainable agriculture 2006 and beyond. New York: Springer; 2007. p. 215–9.

    Google Scholar 

  31. Zhang SZ, Zhang F, Hua BZ. Enhancement of phenylalanine ammonia lyase, polyphenoloxidase, and peroxidase in cucumber seedlings by Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae) infestation. Agric Sci. 2008;7(1):82–7.

    CAS  Google Scholar 

  32. Wei Z, Hu W, Lin QS, Cheng XY, Tong MJ, Zhu LL, et al. Understanding rice plant resistance to the brown planthopper (Nilaparvata lugens): a proteomic approach. Proteomics. 2009;9(10):2798–808.

    Article  CAS  PubMed  Google Scholar 

  33. Du B, Wei Z, Wang ZQ, Wang XX, Peng XX, Du B, et al. Phloem-exudate proteome analysis of response to insect brown planthopper in rice. J Plant Physiol. 2015;183:13–22.

    Article  CAS  PubMed  Google Scholar 

  34. Royo J, Vancanneyt G, Pérez AG, Sanz C, Störmann K, Rosahl S, et al. Characterization of three potato lipoxygenases with distinct enzymatic activities and different organ-specific and wound-regulated expression patterns. J Biol Chem. 1996;271(35):21012–9.

    Article  CAS  PubMed  Google Scholar 

  35. Constabel CP, Ryan CA. A survey of wound- and methyl jasmonate-induced leaf polyphenol oxidase in crop plants. Phytochemistry. 1998;47(4):507–11.

    Article  CAS  Google Scholar 

  36. Constabel CP, Bergey DR, Ryan CA. Systemin activates synthesis of wound-inducible tomato leaf polyphenol oxidase via the octadecanoid defense signaling pathway. Proc Nati Acad Sci USA. 1995;92(2):407–11.

    Article  CAS  Google Scholar 

  37. Zebelo S, Piorkowski J, Disi J, Fadamiro H. Secretions from the ventral eversible gland of Spodoptera exigua caterpillars activate defense-related genes and induce emission of volatile organic compounds in tomato, Solanum lycopersicum. BMC Plant Biol. 2014;14(1):140.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Breusegem FV, Vranová E, Dat JF, Inzé D. The role of active oxygen species in plant signal transduction. Plant Sci. 2001;161(3):405–14.

    Article  Google Scholar 

  39. Yu T, Li YS, Chen XF, Hu J, Chang X, Zhu YG. Transgenic tobacco plants overexpressing cotton glutathione S-transferase (GST) show enhanced resistance to methyl viologen. J Plant Physiol. 2003;160(11):1305–11.

    Article  CAS  PubMed  Google Scholar 

  40. Gallé A, Csiszár J, Secenji M, Guóth A, Cseuz L, Tari I, et al. Glutathione transferase activity and expression patterns during grain filling in flag leaves of wheat genotypes differing in drought tolerance: response to water deficit. J Plant Physiol. 2009;166(17):1878–91.

    Article  PubMed  CAS  Google Scholar 

  41. Santos CVD, Rey P. Plant thioredoxins are key actors in the oxidative stress response. Trends Plant Sci. 2006;11(7):329–34.

    Article  CAS  Google Scholar 

  42. Reichheld JP, Mestres-Ortega D, Laloi C, Meyer Y. The multigenic family of thioredoxin h in Arabidopsis thaliana: specific expression and stress response. Plant Physiol Biochem. 2002;40(6–8):685–90.

    Article  CAS  Google Scholar 

  43. Laloi C, Mestresortega D, Marco Y, Meyer Y, Reichheld JP. The Arabidopsis cytosolic thioredoxin h5 gene induction by oxidative stress and its W-box-mediated response to pathogen elicitor. Plant Physiol. 2004;134(3):1006–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Zhou G, Wang X, Yan F, Wang X, Li R, Cheng J, et al. Genome-wide transcriptional changes and defence-related chemical profiling of rice in response to infestation by the rice striped stem borer Chilo Suppressalis. Physiol Plant. 2011;143(1):21–40.

    Article  CAS  PubMed  Google Scholar 

  45. Albeverio S, Khrennikov A. Biosynthesis and metabolism of jasmonates. J Plant Growth Regul. 2004;23(3):179–99.

    Article  CAS  Google Scholar 

  46. Gardner HW. Biological roles and biochemistry of the lipoxygenase pathway. Hortscience. 1995;30(2):197–205.

    CAS  Google Scholar 

  47. Grechkin A. Recent developments in biochemistry of the plant lipoxygenase pathway. Prog Lipid Res. 1998;37(5):317–52.

    Article  CAS  PubMed  Google Scholar 

  48. Fidantsef AL, Stout MJ, Thaler JS, Duffey SS, Bostock RM. Signal interactions in pathogen and insect attack: expression of lipoxygenase, proteinase inhibitor II, and pathogenesis-related protein P4 in the tomato, Lycopersicon esculentum. Physiol Mol Plant Pathol. 1999;54(3):97–114.

    Article  CAS  Google Scholar 

  49. Wang B, Hajano JUD, Ren Y, Lu C, Wang X. iTRAQ-based quantitative proteomics analysis of rice leaves infected by rice stripe virus reveals several proteins involved in symptom formation. Virol J. 2015;12(1):99–119.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Stintzi A, Weber H, Reymond P, Browse J, Farmer EE. Plant defense in the absence of jasmonic acid: the role of cyclopentenones. Proce Nati Acad Sci USA. 2001;98(22):12837–42.

    Article  CAS  Google Scholar 

  51. Johnson PR, Ecker JR. The ethylene gas signal transduction pathway in plants: a molecular perspective. Annu Rev Genet. 1998;32(32):227–54.

    Article  CAS  PubMed  Google Scholar 

  52. Yang SF, Hoffman NE. Ethylene biosynthesis and its regulation in higher plants. Annu Rev Plant Physiol. 2003;35(1):155–89.

    Article  Google Scholar 

  53. Niu YY, Chen M, Xu ZS, Li LC, Chen XP, Ma YZ. Characterization of ethylene receptors and their interactions with GmTPRA novel tetratricopeptide repeat protein (TPR) in soybean (Glycine max L.). JIA. 2013;12(4):571–81.

    Google Scholar 

  54. Nole-Wilson S, Krizek BA. DNA binding properties of the Arabidopsis floral development protein AINTEGMENGTA. Nucleic Acids Res. 2000;28(21):4076–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Zhang Z, Huang R. Enhanced tolerance to freezing in tobacco and tomato overexpressing transcription factor TERF2/LeERF2 is modulated by ethylene biosynthesis. Plant Mol Biol. 2010;73(3):241–9.

    Article  CAS  PubMed  Google Scholar 

  56. Monroy AF, Dhindsa RS. Low temperature signal transduction: induction of cold acclimation specific genes of alfalfa by calcium at 25 degrees C. Plant Cell. 1995;7(3):321–31.

    CAS  PubMed  PubMed Central  Google Scholar 

  57. Guilfoyle TJ. Aux/IAA proteins and auxin signal transduction. Trends Plant Sci. 1998;3(6):205–7.

    Article  Google Scholar 

  58. Guilfoyle TJ, Hagen G. Auxin response factors. Curr Opin Plant Biol. 2007;10(5):453–60.

    Article  CAS  PubMed  Google Scholar 

  59. Hagen G, Guilfoyle T. Auxin-responsive gene expression: genes, promoters and regulatory factors. Plant Mol Biol. 2002;49(3):373–85.

    Article  CAS  PubMed  Google Scholar 

  60. Park JE, Park JY, Kim YS, Staswick PE, Jeon J, Yun J, et al. GH3-mediated auxin homeostasis links growth regulation with stress adaptation response in Arabidopsis. J Biol Chem. 2007;282(13):10036–46.

    Article  CAS  PubMed  Google Scholar 

  61. Kant S, Bi YM, Zhu T, Rothstein SJ. SAUR39, a small auxin-up RNA gene, acts as a negative regulator of auxin synthesis and transport in rice. Plant Physiol. 2009;151(2):691–701.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Kant S, Rothstein S. Auxin-responsive SAUR39 gene modulates auxin level in rice. Plant Signal Behav. 2009;4(12):1174–5.

    Article  PubMed  PubMed Central  Google Scholar 

  63. León J, Rojo E, Sánchez-Serrano J. Wound signaling in plants. J Epx Bot. 2001;52(354):1–9.

    Article  Google Scholar 

  64. Afzal AJ, Wood AJ, Lightfoot DA. Plant receptor-like serine threonine kinases: roles in signaling and plant defense. Mol Plant-Microbe Interact. 2008;21(5):507–17.

    Article  CAS  PubMed  Google Scholar 

  65. Hasegawa PM, Bressan RA, Zhu JK, Bohnert HJ. Plant cellular and molecular responses to high salinity. Annu Rev Plant Physiol Plant Mol Biol. 2000;51:463–99.

    Article  CAS  PubMed  Google Scholar 

  66. Zipfel C. Pattern-recognition receptors in plant innate immunity. Curr Opin Immunol. 2008;20(1):10–6.

    Article  CAS  PubMed  Google Scholar 

  67. Chinnusamy V, Schumaker K, Zhu JK. Molecular genetic perspectives on cross-talk and specificity in abiotic stress signaling in plants. J Exp Bot. 2004;55(395):225–36.

    Article  CAS  PubMed  Google Scholar 

  68. Dey S, Ghose K, Basu D. Fusarium elicitor-dependent calcium influx and associated ros generation in tomato is independent of cell death. J Exp Bot. 2010;126(2):217–28.

    Google Scholar 

  69. Maffei ME, Mithöfer A, Arimura GI, Uchtenhagen H, Bossi S, Bertea CM, et al. Effects of feeding Spodoptera littoralis on lima bean leaves. III. Membrane depolarization and involvement of hydrogen peroxide. Plant Physiol. 2006;140(3):1022–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Mohanta TK, Occhipinti A, Zebelo SA, Foti M, Fliegmann J, Bossi S, et al. Ginkgo biloba responds to herbivory by activating early signaling and direct defenses. PLoS One. 2012;7(3):e32822.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Karban R, Kuć J. Induced resistance against pathogens and herbivores: an overview, Induced plant defenses against pathogens and herbivores; 1999. p. 1–15.

    Google Scholar 

  72. Alagar M, Suresh S, Saravanakumar D, Samiyappan R. Feeding-induced changes in defence enzymes and proteins and their implications in host resistance to Nilaparvata lugens. J Appl Entomol. 2010;134(2):123–31.

    Article  CAS  Google Scholar 

  73. Chen H, Wilkerson CG, Kuchar JA, Phinney BS, Howe GA. Jasmonate-inducible plant enzymes degrade essential amino acids in the herbivore midgut. Proc Nati Acad Sci USA. 2005;102(52):19237–42.

    Article  CAS  Google Scholar 

  74. Stensjo K, Pettersson J, Bryngelsson T, Jonsson L. Aphid infestation induces PR-proteins differently in barley susceptible or resistant to the birdcherry-oat aphid (Rhopalosiphum padi). Physiol Plant. 2008;110(4):496–502.

    Google Scholar 

  75. Van Loon LC, Van Strien EA. The families of pathogenesis related proteins, their activities, and comparative analysis of PR-1 type proteins. Physiol Mol Plant Pathol. 1999;55(2):85–97.

    Article  CAS  Google Scholar 

  76. Jwa NS, Agrawal GK, Rakwal R, Park CH, Agrawal VP. Molecular cloning and characterization of a novel jasmonate inducible pathogenesis-related class 10 protein gene, JIOsPR10, from rice (Oryza sativa L.) seedling leaves. Biochem Biophys Res Commun. 2001;286(5):973–83.

    Article  CAS  PubMed  Google Scholar 

  77. Uknes S, Winter AM, Delaney T, Vernooij B, Morse A, Friedrich L, et al. Biological induction of systemic acquired resistance in Arabidopsis. Mol Plant-Microbe Ineract. 1993;6(6):692–8.

    Article  Google Scholar 

  78. Green TR, Ryan CA. Wound induced proteinase inhibitors in plant leaves: a possible defense mechanism against insects. Science. 1972;175(4023):776–7.

    Article  CAS  PubMed  Google Scholar 

  79. Tamhane VA, Chougule NP, Giri AP, Dixit AR, Sainani MN, Gupta VS. In vivo and in vitro effect of Capsicum Annum proteinase inhibitors on Helicoverpa armigera gut proteinases. Biochimica Et Biophvsica Acta. 2005;1722(2):156–67.

    Article  CAS  Google Scholar 

  80. Gatehouse AMR, Hilder VA, Boulter D. Plant genetic manipulation for crop protection. Biotch Agric. 1992;22(6):1408–9(2).

    Google Scholar 

  81. Rao R, Fiandra L, Giordana B, De EM, Congiu T, Burlini N, et al. AcMNPV ChiA protein disrupts the peritrophic membrane and alters midgut physiology of Bombyx Mori larvae. Insect Biochem Mol Biol. 2004;34(11):1205–13.

    Article  CAS  PubMed  Google Scholar 

  82. Ding X, Gopalakrishnan B, Johnson LB, White FF, Wang XR, Morgan TD, et al. Insect resistance of transgenic tobacco expressing an insect chitinase gene. Transgenic Res. 1998;7(2):77–84.

    Article  CAS  PubMed  Google Scholar 

  83. Zhu-Salzman K, Salzman RA, Ahn JE, Koiwa H. Transcriptional regulation of sorghum defense determinants against a phloem-feeding aphid. Plant Physiol. 2004;134(1):420–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Kant MR, Ament K, Sabelis MW, Haring MA, Schuurink RC. Differential timing of spider mite-induced direct and indirect defenses in tomato plants. Plant Physiol. 2004;135(1):483–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Wang HB, Jin S, Tao Y. Insect resistance of the chitinase in Vicia faba leaves to Aphis craccivove. J Fudan University (Nat Sci). 1994;33(3):348–52.

    CAS  Google Scholar 

  86. Peumans WJ, Van Damme EJ. Lectins as plant defense proteins. Plant Physiol. 1995;109(2):347–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Komath SS, Kavitha M, Swamy MJ. Beyond carbohydrate binding: new direct ions in plant lectin research. Org Biomol Chem. 2006;4(6):973–88.

    Article  CAS  PubMed  Google Scholar 

  88. Boulter D, Edwards GA, Gatehouse AMR, Gatehouse JA, Hilder VA. Additive protective effects of incorporating two different higher plant insect resistence genes in transgenic tobacco plants. Corp Protect. 1990;9:351–4.

    Article  Google Scholar 

  89. Hilder VA, Powell KS, Gatehouse AMR, Gatehouse JA, Gatehouse LN, Shi Y, et al. Expression of snowdrop lectin in transgenic tobacco plants results in added protection against aphids. Transgenic Res. 1995;4(1):18–25.

    Article  CAS  Google Scholar 

  90. Van Damme EJM. Handbook of plant lectins: properties and biomedical applications. Chester: John Willey and Sons Ltd; 1998. p. 452.

    Google Scholar 

  91. Chrispeels MJ, Raikhel NV. Lectins, lectin genes, and their role in plant defense. Plant Cell. 1991;3(1):1–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. Lange BM, Rujan T, Martin W, Croteau R. Isoprenoid biosynthesis: the evolution of two ancient and distinct pathways across genomes. Proc Nati Acad Sci USA. 2000;97(24):13172–7.

    Article  CAS  Google Scholar 

  93. Herrmann KM. The shikimate pathway: early steps in the biosynthesis of aromatic compounds. Plant Cell. 1995;7(3):907–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Chaman ME, Copaja SV, Argandoña VH. Relationships between salicylic acid content, phenylalanine ammonia-lyase (PAL) activity, and resistance of barley to aphid infestation. J Agric Food Chem. 2003;51(8):2227–31.

    Article  CAS  PubMed  Google Scholar 

  95. Ritter H, Schulz GE. Structural basis for the entrance into the phenylpropanoid metabolism catalyzed by phenylalanine ammonia-lyase. Plant Cell. 2004;16(12):3426–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  96. Li JB, Fang LP, Zhang YN, Yang WJ, Guo Q, Li L, et al. The relationship between the resistance of ctton against cotton aphid, aphis gossypii, and the activity of phenylalanine ammonia-lyase. Chin Bull Entomol. 2008;45(3):422–5.

    CAS  Google Scholar 

  97. Hartley SE, Firn RD. Phenolic biosynthesis, leaf damage, and insect herbivory in birch (Betula Pendula). J Chem Ecol. 1989;15(1):275–83.

    Article  CAS  PubMed  Google Scholar 

  98. Kühnl T, Koch U, Heller RW, Wellmann E. Elicitor induced S-adenosyl-l-methionine: caffeoyl-CoA 3-O-methyltransferase from carrot cell suspension cultures. Plant Sci. 1989;60(1):21–5.

    Article  Google Scholar 

  99. Pakusch AE, Kneusel RE, Matern U. S-adenosyl-l-methionine: trans-caffeoyl-coenzyme a 3-O-methyltransferase from elicitor-treated parsley cell suspension culture. Arch Biochem Biophys. 1989;271(2):488–94.

    Article  CAS  PubMed  Google Scholar 

  100. Gowri G, Bugos RC, Campbell WH, Maxwell CA, Dixon RA. Stress responses in alfalfa (Medicago sativa L.): X. Molecular cloning and expression of S-adenosyl-l-methionine: caffeic acid 3-O-methyltransferase, a key enzyme of lignin biosynthesis. Plant Physiol. 1991;97(1):7–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Wyrambik D, Grisebach H. Purification and properties of isoenzymes of cinnamyl-alcohol dehydrogenase from soybean-cell-suspension cultures. Eur J Biochem. 1975;59(1):9–15.

    Article  CAS  PubMed  Google Scholar 

  102. Mansell RL, Gross GG, Stǒckigt J, Franke H, Zenk MH. Purification and properties of cinnamyl alcohol dehydrogenase from higher plants involved in lignin biosynthesis. Phytochemistry. 1974;13(11):2427–35.

    Article  CAS  Google Scholar 

  103. Henzell RF, Hall WT. Substituted phenols as repellents for male Costelytra zealandica (Goleoptera: Scarabaeidae) in mating flights. New Zealand J Zool. 1974;1(4):509–13.

    Article  CAS  Google Scholar 

  104. Schroeder FC, del Campo ML, Grant JB, Weibel DB, Smedley SR, Bolton KL, et al. Pinoresinol: a lignol of plant origin serving for defense in a caterpillar. Proc Nati Acad Sci USA. 2006;103(42):15497–501.

    Article  CAS  Google Scholar 

  105. Winkel-Shirley B. Biosynthesis of flavonoids and effects of stress. Curr Opin Plant Biol. 2002;5(3):218–23.

    Article  CAS  PubMed  Google Scholar 

  106. Koes RE, Quattrocchio F, Mol JNM. The flavonoid biosynthetic pathway in plants: function and evolution. BioEssays. 1994;16(2):123–32.

    Article  CAS  Google Scholar 

  107. Martin CR. Structure, function, and regulation of the chalcone synthase. Int Rev Cytol. 1993;147:233–84.

    Article  CAS  PubMed  Google Scholar 

  108. Muir SR, Collins GJ, Robinson S, Hughes S, Bovy A, Ric De Vos CH, et al. Overexpression of petunia chalcone isomerase in tomato results in fruit containing increased levels of flavonols. Nat Biotechnol. 2001;19(5):470–4.

    Article  CAS  PubMed  Google Scholar 

  109. Johnson ET, Ryu S, Yi H, Shin B, Cheong H, Choi G. Alteration of a single amino acid changes the substrate specificity of dihydroflavonol 4-reductase. Plant J. 2001;25(3):325–33.

    Article  CAS  PubMed  Google Scholar 

  110. Kumar V, Yadav SK. Overexpression of CsDFR and CsANR enhanced flavonoids accumulation and antioxidant potential of roots in tobacco. Plant Roots. 2013;7:65–76.

    Article  CAS  Google Scholar 

  111. Schuler MA, Daniele WR. Functional genomics of P450s. Annu Rev Plant Biol. 2003;54(1):627–9.

    Article  CAS  Google Scholar 

  112. Harvey PJ, Campanella BF, Castro PM, Harms H, Lichtfouse E, Schäffner AR, et al. Phytoremediation of polyaromatic hydrocarbons, anilines and phenols. Environ Sci Pollut Res Int. 2002;9(1):29–47.

    Article  CAS  PubMed  Google Scholar 

  113. Morant M, Bak S, Møller BL, Werck-Reichhart D. Plant cytochromes P450: tools for pharmacology, plant protection and phytoremediation. Curr Opin Biotechnol. 2003;14(2):151–62.

    Article  CAS  PubMed  Google Scholar 

  114. Bansal R, Mian M, Mittapalli O, Miche AP. RNA-Seq reveals a xenobiotic stress response in the soybean aphid, Aphis glycines, when fed aphid-resistant soybean. BMC Genomics. 2014;15:972.

    Article  PubMed  PubMed Central  Google Scholar 

  115. Kahn RA, Bak S, Svendsen J, Halkier BA, Møller BL. Isolation and reconstitution of cytochrome P450ox and in vitro reconstitution of the entire biosynthetic pathway of the cyanogenic glucoside dhurrin from sorghum. Plant Physiol. 1997;115(4):1661–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  116. Sanmartin M, Pateraki L, Chatzopoulou F, Kanellis AK. Differential expression of the ascorbate oxidase multigene family during fruit development and in response to stress. Planta. 2007;225(4):873–85.

    Article  CAS  PubMed  Google Scholar 

  117. Gaspard S, Monzani E, Casella L, Gullotti M, Maritano S, Marchesini A. Inhibition of aseorbate oxidase by phenolic compounds. Enzymatic and spectroscopic studies. Biochemistry. 1997;36(16):4852–9.

    Article  CAS  PubMed  Google Scholar 

  118. Sanmartin M, Drogoudi PD, Lyons T, Pateraki I, Barnes J, Kanellis AK. Over-expression of ascorbate oxidase in the apoplast of transgenic tobacco results in altered ascorbate and glutathione redox states and increased sensitivity to ozone. Planta. 2003;216(6):918–28.

    CAS  PubMed  Google Scholar 

  119. Barbehenn RV, Jaros A, Yip L, Tran L, Kanellis AK, Constabel CP. Evaluating ascorbate oxidase as a plant defense against leaf-chewing insects using transgenic poplar. J Chem Ecol. 2008;34(10):1331–40.

    Article  CAS  PubMed  Google Scholar 

  120. Fotopoulos V, Sanmartin M, Kanellis A. Effect of ascorbate oxidase over-expression on ascorbate recycling gene expression in response to agents imposing oxidative stress. J Exp Bot. 2006;57(14):3933–43.

    Article  CAS  PubMed  Google Scholar 

  121. Eulgem T, Rushton PJ, Schmelzer E, Hahlbrock K, Somssich IE. Early nuclear events in plant defence signalling: rapid gene activation by WRKY transcription factors. EMBO J. 1999;18(17):4689–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  122. Eulgem T, Rushton PJ, Robatzek S, Somssich IE. The WRKY superfamily of plant transcription factors. Trends Plant Sci. 2000;5(5):199–206.

    Article  CAS  PubMed  Google Scholar 

  123. Rushton PJ, Somssich IE, Ringler P, Shen QJ. WRKY transcription factors. Trends Plant Sci. 2010;15(5):247–58.

    Article  CAS  PubMed  Google Scholar 

  124. Zhou X, Jiang Y, Yu D. WRKY22 transcription factor mediates dark-induced leaf senescence in Arabidopsis. Mol Cell. 2011;31(4):303–13.

    Article  CAS  Google Scholar 

  125. Grunewald W, Karimi M, Wieczorek K, Van de Cappelle E, Chnitzki E, Grundler F, et al. A role for AtWRXY23 in feeding site establishment of plant-parasitic nematodes. Plant Physiol. 2008;148:358–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  126. Artico S, Ribeiro-Alves M, Oliveira-Neto OB, Leonardo Macedo LPD, Silveira S, et al. Transcriptome analysis of Gossypium hirsutum flower buds infested by cotton boll weevil (Anthonomus grandis) larvae. BMC Genomics. 2014;15(854):1–24.

    Google Scholar 

  127. Fang Y, You J, Xie K, Xie W, Xiong L. Systematic sequence analysis and identification of tissue-specific or stress-responsive genes of NAC transcription factor family in rice. Mol Gen Genomics. 2008;280(6):547–63.

    Article  CAS  Google Scholar 

  128. Nakashima K, Takasaki H, Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. NAC transcription factors in plant abiotic stress responses. Biochim Biophys Acta. 2012;1819(2):97–103.

    Article  CAS  PubMed  Google Scholar 

  129. Wang Z, Dane F. NAC (NAM/ATAF/CUC) transcription factors in different stresses and their signaling pathway. Acta Physiol Plant. 2013;35(5):1397–408.

    Article  CAS  Google Scholar 

  130. Hegedus D, Yu M, Baldwin D, Gruber M, Sharpe A, Parkin I, et al. Molecular characterization of brassicanapus NAC domain transcriptional activators induced in response to biotic and abiotic stresses. Plant Mol Biol. 2003;53(3):383–97.

    Article  CAS  PubMed  Google Scholar 

  131. Nogueira FTS, Schlogl PS, Camargo SR, Fernandez JH, de Jr Rosa VE, Pompermayer P, et al. SsNAC23, a member of the NAC domain protein family, is associated with cold, herbivory and water stress in sugarcane. Plant Sci. 2005;169(1):93–106.

    Article  CAS  Google Scholar 

  132. Kaur H, Heinzel N, Schöttner M, Baldwin IT, Gális I. R2R3-NaMYB8 regulates the accumulation of phenylpropanoid-polyamine conjugates, which are essential for local and systemic defense against insect herbivores in Nicotiana attenuata. Plant Physiol. 2010;152(3):1731–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  133. Toledo-Ortiz G, Huq E, Quail PH. The Arabidopsis basic/helix-loop-helix transcription factor family. Plant Cell. 2003;15(8):1749–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  134. Van Verk MC, Gatz C, Linthorst HJM. Transcriptional regulation of plant defense responses. J Regul Econ. 2009;20(2):191–211.

    Google Scholar 

  135. Fujimoto SY, Ohta M, Usui A, Shinshi H, Ohme-Takagi M. Arabidopsis ethylene-responsive element binding factors act as transcriptional activators or repressors of GCC box-mediated gene expression. Plant Cell. 2000;12(3):393–404.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  136. Leubner-Metzger G, Petruzzelli L, Waldvogel R, Vögeli-Lange R, Meins FJR. Ethylene responsive element binding protein (EREBP) expression and the transcriptional regulation of class I beta-1, 3-glucanase during tobacco seed germination. Plant Mol Biol. 1998;38(5):785–95.

    Article  CAS  PubMed  Google Scholar 

  137. Liu Q, Kasuga M, Sakuma Y, Miura S, Yamaguchi-Shinozaki K, Shinozaki K. Two transcription factors, DREB1 and DREB2, with an EREBP/AP2 DNA-binding domain separate two cellular signal transduction pathways in drought- and low-temperature-responsive gene expression in Arabidopsis. Plant Cell. 1998;10(8):1391–406.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  138. Park JM, Park CJ, Lee SB, Ham BK, Shin R, Paek KH. Overexpression of the tobacco Tsi1 gene encoding an EREBP/AP2-type transcription factor enhances resistance against pathogen attack and osmotic stress in tobacco. Plant Cell. 2001;13(5):1035–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  139. Jia W, Qiu K, He M, Song P, Zhou Q, Zhou F, et al. SOAPfuse: an algorithm for identifying fusion transcripts from paired-end RNA-Seq data. Genome Biol. 2013;14(2):R12.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  140. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-Seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7(3):562–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  141. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bionformatics. 2011;12(1):323.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

Not applicable

Funding

This work was supported by the Natural Science Fund of Guangxi (2013GXNSFDA019009; 2016GXNSFAA380238), and the Development Foundation of Guangxi Academy of Agricultural Sciences (2016JM03). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

All data were submitted to the National Center for Biotechnology Information (NCBI) under SRA number SRA549176.

Author information

Authors and Affiliations

Authors

Contributions

ZDS and WYZ conceived and designed the experiments. WYZ, ZYC, HZC, ZGL, SZY and XMT performed the experiments. WYZ analyzed the data. ZDS and WYZ contributed reagents/materials/analysis tools. WYZ and ZDS conceived the experiments and wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Zudong Sun.

Ethics declarations

Ethics approval and consent to participate

Not applicable for this research. Soybean seeds for this study were obtained from the Guangxi Academy of Agricultural Science.

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Table S1.

Total number of DEGs between HRK48 and HRK0. (XLS 163 kb)

Additional file 2: Table S2.

Total number of DEGs between HSK48 and HSK0. (XLS 111 kb)

Additional file 3: Table S3.

Total number of DEGs between HRK0 and HSK0. (XLS 119 kb)

Additional file 4: Table S4.

Total number of DEGs between HRK48 and HSK48. (XLS 101 kb)

Additional file 5: Table S5.

GO analysis the DEGs of HRK0-VS-HRK48. (XLS 164 kb)

Additional file 6: Table S6.

GO analysis the DEGs of HSK0-VS-HSK48. (XLS 141 kb)

Additional file 7: Table S7.

GO analysis the DEGs of HSK0-VS-HRK0. (XLS 137 kb)

Additional file 8: Table S8.

GO analysis the DEGs of HSK48-VS-HRK48. (XLS 119 kb)

Additional file 9: Table S9.

Pathway analysis the DEGs of HRK0-VS-HRK48. (XLS 108 kb)

Additional file 10: Table S10.

Pathway analysis the DEGs of HSK0-VS-HSK48. (XLS 92 kb)

Additional file 11: Table S11.

Pathway analysis the DEGs of HSK0-VS-HRK0. (XLS 90 kb)

Additional file 12: Table S12.

Pathway analysis the DEGs of HSK48-VS-HRK48. (XLS 90 kb)

Additional file 13: Table S13.

DEGs confirmed by the qRT-PCR. (XLS 119 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zeng, W., Sun, Z., Cai, Z. et al. Comparative transcriptome analysis of soybean response to bean pyralid larvae. BMC Genomics 18, 871 (2017). https://doi.org/10.1186/s12864-017-4256-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-017-4256-7

Keywords