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.

Keywords: Bean pyralid; Differentially expressed genes (DEGs); Soybean; Transcriptome sequencing.

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

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

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

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

Fig. 5

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

Fig. 6

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

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

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