- Research article
- Open Access
Comparative transcriptome analysis of soybean response to bean pyralid larvae
BMC Genomicsvolume 18, Article number: 871 (2017)
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.
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.
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.
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 . 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 . So bean pyralid is different from other leaf-feeding insects with chewing mouthparts which cause holes or incisions by means of encroachment . 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 . 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 . 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 . 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 . 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.
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.
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.
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.
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.
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).
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).
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).
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).
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 . 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 . Research demonstrated that plant transcription factors, including NAC, MYB, WRKY, were involved in the plant defense responses . 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).
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.
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 . Meanwhile, POD has been found to have significant insecticidal effects on Lepidoptera and Coleoptera . Previous studies have shown that POD could be induced by insects in wheat , sorghum , cucumber  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 . 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 . 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 , and plays an important signal factor in plant induction defense pathways . 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 . 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 . 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 . Ethylene receptor is an upstream component of ethylene signal transduction pathway, that plays an important role in plant growth and responses to adversity stresses . 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 . 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 . 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 . 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 . 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  and Ginkgo biloba . 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 .
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 . 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 . 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 . When insect eats plants with sustained chitinase expression, its digestive tract becomes damaged and its epidermis cannot form normally . Chitinase could be induced by aphid and spider mites in sorghum  and tomato , 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 .
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 , Myzus persicae  and Lacanobia oleracea  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 .
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 . 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 . 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 . 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 . 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 . 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 . 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 . Previous studies have shown that over-expressions of DFR could improve the content and antioxidant capacities of flavonoid . 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 . CYP genes can be induced by aphid in soybean . 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 . 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 . For example, one WRKY23 gene was induced by Heterodera schachtii in Arabidopsis thaliana  and 20 WRKY transcription factors were induced by cotton boll weevil in Gossypium hirsutum . 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 , sugarcane  and Gossypium hirsutum . MYB is one of the central regulators of development, metabolism, and response to abiotic and biotic stresses . 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 . MYB could be induced by small cabbage white caterpillars  and cotton boll weevil . 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 . 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 . 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.
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.
The tested materials Gantai-2-2 (highly resistant line)  and Wan82–178 (highly susceptible line)  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 . 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/) .
Genes and isoforms expression level are quantified by a software package: RSEMV1.2.12 (RNASeq by Expectation Maximization) . 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:
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.
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.
- ACC oxidase:
Aminocyclopropane carboxylate oxidase
Caffeoyl-coA 3-O-methyl transferase
Caffeic acid 3-O-methyltransferase
Differentially expressed genes
Ethylene responsive factor 1
Kyoto Encyclopedia of Genes and Genomes
12-oxophytodienoic acid reductase
Quantitative real time-PCR
Reactive oxygen species
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.
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.
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.
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.
Sun ZD, Chen HZ, Wei DW. A study on leaf-feeding insect species on soybeans in Nanning. Guangxi Agric Sci. 2001;2:104–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.
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.
Xing GN, Zhao TJ, Gai JY. Inheritance of resistance to Lamprosema indicata Fabricius in soybean. Acta Agron Sin. 2008;34(1):8–16.
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
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.
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.
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.
Foyer CH, Shigeoka S. Understanding oxidative stress and antioxidant functions to enhance photosynthesis. Plant Physiol. 2011;155(1):93–100.
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.
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.
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.
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.
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.
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.
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.
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.
Schuler MA. The role of cytochrome P450 monooxygenases in plant-insect interactions. Plant Physiol. 1996;112(4):1411–9.
Treutter D. Significance of flavonoids in plant resistance: a review. Environ Chem Lett. 2006;4(3):147–57.
Korkina LG. Phenylpropanoids as naturally occurring antioxidants: from plant defense to human health. Cell Mol Biol. 2007;53(1):15–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.
Rushton PJ, Somssich IE. Transcriptional control of plant genes to responsive to pathogens. Curr Opin Plant Biol. 1998;1(4):311–5.
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.
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.
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.
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.
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.
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.
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.
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.
Constabel CP, Ryan CA. A survey of wound- and methyl jasmonate-induced leaf polyphenol oxidase in crop plants. Phytochemistry. 1998;47(4):507–11.
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.
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.
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.
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.
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.
Santos CVD, Rey P. Plant thioredoxins are key actors in the oxidative stress response. Trends Plant Sci. 2006;11(7):329–34.
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.
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.
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.
Albeverio S, Khrennikov A. Biosynthesis and metabolism of jasmonates. J Plant Growth Regul. 2004;23(3):179–99.
Gardner HW. Biological roles and biochemistry of the lipoxygenase pathway. Hortscience. 1995;30(2):197–205.
Grechkin A. Recent developments in biochemistry of the plant lipoxygenase pathway. Prog Lipid Res. 1998;37(5):317–52.
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.
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.
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.
Johnson PR, Ecker JR. The ethylene gas signal transduction pathway in plants: a molecular perspective. Annu Rev Genet. 1998;32(32):227–54.
Yang SF, Hoffman NE. Ethylene biosynthesis and its regulation in higher plants. Annu Rev Plant Physiol. 2003;35(1):155–89.
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.
Nole-Wilson S, Krizek BA. DNA binding properties of the Arabidopsis floral development protein AINTEGMENGTA. Nucleic Acids Res. 2000;28(21):4076–82.
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.
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.
Guilfoyle TJ. Aux/IAA proteins and auxin signal transduction. Trends Plant Sci. 1998;3(6):205–7.
Guilfoyle TJ, Hagen G. Auxin response factors. Curr Opin Plant Biol. 2007;10(5):453–60.
Hagen G, Guilfoyle T. Auxin-responsive gene expression: genes, promoters and regulatory factors. Plant Mol Biol. 2002;49(3):373–85.
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.
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.
Kant S, Rothstein S. Auxin-responsive SAUR39 gene modulates auxin level in rice. Plant Signal Behav. 2009;4(12):1174–5.
León J, Rojo E, Sánchez-Serrano J. Wound signaling in plants. J Epx Bot. 2001;52(354):1–9.
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.
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.
Zipfel C. Pattern-recognition receptors in plant innate immunity. Curr Opin Immunol. 2008;20(1):10–6.
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.
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.
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.
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.
Karban R, Kuć J. Induced resistance against pathogens and herbivores: an overview, Induced plant defenses against pathogens and herbivores; 1999. p. 1–15.
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.
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.
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.
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.
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.
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.
Green TR, Ryan CA. Wound induced proteinase inhibitors in plant leaves: a possible defense mechanism against insects. Science. 1972;175(4023):776–7.
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.
Gatehouse AMR, Hilder VA, Boulter D. Plant genetic manipulation for crop protection. Biotch Agric. 1992;22(6):1408–9(2).
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.
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.
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.
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.
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.
Peumans WJ, Van Damme EJ. Lectins as plant defense proteins. Plant Physiol. 1995;109(2):347–52.
Komath SS, Kavitha M, Swamy MJ. Beyond carbohydrate binding: new direct ions in plant lectin research. Org Biomol Chem. 2006;4(6):973–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.
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.
Van Damme EJM. Handbook of plant lectins: properties and biomedical applications. Chester: John Willey and Sons Ltd; 1998. p. 452.
Chrispeels MJ, Raikhel NV. Lectins, lectin genes, and their role in plant defense. Plant Cell. 1991;3(1):1–9.
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.
Herrmann KM. The shikimate pathway: early steps in the biosynthesis of aromatic compounds. Plant Cell. 1995;7(3):907–19.
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.
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.
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.
Hartley SE, Firn RD. Phenolic biosynthesis, leaf damage, and insect herbivory in birch (Betula Pendula). J Chem Ecol. 1989;15(1):275–83.
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.
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.
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.
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.
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.
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.
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.
Winkel-Shirley B. Biosynthesis of flavonoids and effects of stress. Curr Opin Plant Biol. 2002;5(3):218–23.
Koes RE, Quattrocchio F, Mol JNM. The flavonoid biosynthetic pathway in plants: function and evolution. BioEssays. 1994;16(2):123–32.
Martin CR. Structure, function, and regulation of the chalcone synthase. Int Rev Cytol. 1993;147:233–84.
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.
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.
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.
Schuler MA, Daniele WR. Functional genomics of P450s. Annu Rev Plant Biol. 2003;54(1):627–9.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Eulgem T, Rushton PJ, Robatzek S, Somssich IE. The WRKY superfamily of plant transcription factors. Trends Plant Sci. 2000;5(5):199–206.
Rushton PJ, Somssich IE, Ringler P, Shen QJ. WRKY transcription factors. Trends Plant Sci. 2010;15(5):247–58.
Zhou X, Jiang Y, Yu D. WRKY22 transcription factor mediates dark-induced leaf senescence in Arabidopsis. Mol Cell. 2011;31(4):303–13.
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.
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.
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.
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.
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.
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.
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.
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.
Toledo-Ortiz G, Huq E, Quail PH. The Arabidopsis basic/helix-loop-helix transcription factor family. Plant Cell. 2003;15(8):1749–70.
Van Verk MC, Gatz C, Linthorst HJM. Transcriptional regulation of plant defense responses. J Regul Econ. 2009;20(2):191–211.
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.
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.
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.
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.
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.
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.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bionformatics. 2011;12(1):323.
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.
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
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Total number of DEGs between HRK48 and HRK0. (XLS 163 kb)
Total number of DEGs between HSK48 and HSK0. (XLS 111 kb)
Total number of DEGs between HRK0 and HSK0. (XLS 119 kb)
Total number of DEGs between HRK48 and HSK48. (XLS 101 kb)
GO analysis the DEGs of HRK0-VS-HRK48. (XLS 164 kb)
GO analysis the DEGs of HSK0-VS-HSK48. (XLS 141 kb)
GO analysis the DEGs of HSK0-VS-HRK0. (XLS 137 kb)
GO analysis the DEGs of HSK48-VS-HRK48. (XLS 119 kb)
Pathway analysis the DEGs of HRK0-VS-HRK48. (XLS 108 kb)
Pathway analysis the DEGs of HSK0-VS-HSK48. (XLS 92 kb)
Pathway analysis the DEGs of HSK0-VS-HRK0. (XLS 90 kb)
Pathway analysis the DEGs of HSK48-VS-HRK48. (XLS 90 kb)
DEGs confirmed by the qRT-PCR. (XLS 119 kb)