International Journal of
Molecular Sciences Article
Transcriptional Response of Silkworm (Bombyx mori) Eggs to O2 or HCl Treatment Jing Gong, Sha Tian, Xia Zhou, Huan Yang, Yong Zhu * and Yong Hou * State Key Laboratory of Silkworm Genome Biology, College of Biotechnology, Southwest University, Chongqing 400715, China; [email protected]
(J.G.); [email protected]
(S.T.); [email protected]
(X.Z.); [email protected]
(H.Y.) * Correspondence: [email protected]
(Y.Z.); [email protected]
(Y.H.); Tel.: +86-23-6825-0191 (Y.Z.); +86-23-6825-0099 (Y.H.); Fax: +86-23-6825-0191 (Y.Z.); +86-23-6825-1128 (Y.H.) Academic Editor: Kun Yan Zhu Received: 9 September 2016; Accepted: 1 November 2016; Published: 7 December 2016
Abstract: Diapause is a common biological phenomenon that occurs in many organisms, including fish, insects, and nematodes. In the silkworm (Bombyx mori), diapause generally occurs in the egg stage. Treatment with O2 , HCl, or other compounds can prevent egg diapause. Here, we characterized the transcriptomic responses of newly laid eggs treated with O2 or HCl. Digital gene expression analysis showed that 610 genes in O2 -treated eggs and 656 in HCl-treated eggs were differentially expressed. Of these, 343 genes were differentially expressed in both treatments. In addition to trehalases, sorbic acid dehydrogenases, and some enzymes involved in the carbohydrate metabolism, we also identified heat shock proteins, cytochrome P450, and GADD45, which are related to stress tolerance. Gene ontology enrichment analysis showed differentially expressed genes in O2 -treated eggs were involved in oxidoreductase activity as well as in binding, catalytic, and metabolic processes. The Kyoto Encyclopedia of Genes and Genomes analysis showed that the pathways for ribosome biogenesis, spliceosome, and circadian rhythm were significantly enriched in HCl-treated eggs. The reliability of the data was confirmed by qRT-PCR analysis. Our results improved the understanding of the mechanism of diapause blocking in silkworm eggs treated with O2 or HCl and identified novel molecular targets for future studies. Keywords: silkworm; diapause; trehalase; RNA biogenesis; heat shock protein
1. Introduction Diapause, defined as a period of arrested development and reduced energy consumption, is a strategy that helps animals survive extreme environmental conditions, such as cold, heat, and drought. Diapause occurs in both vertebrates and invertebrates [1,2]. In insects, diapause occurs at any stage of the development, including egg, larva, pupa, and adult, and the underlying biological mechanism has been studied extensively . The silkworm, Bombyx mori, is an economically important insect with a breeding history of more than 5000 years in China. Diapause in silkworm usually occurs during the late gastrula stage of embryogenesis and is called embryonic diapause [4,5]. Most silkworm strains in northern China produce eggs that undergo diapause and hatching occurs once or twice per year, whereas in southern China, most silkworm strains do not produce diapause eggs under natural conditions and reproduce several times per year . In addition to geography, the environmental conditions of egg development affect whether silkworm eggs enter diapause. For instance, Dazao, a bivoltine silkworm strain, lays eggs that experience diapause if they are incubated under long-day conditions (18 h light:6 h dark) at 25 ◦ C during the egg-development stage, whereas when the eggs are incubated in continuous darkness at 15 ◦ C, they do not experience diapause . Int. J. Mol. Sci. 2016, 17, 1838; doi:10.3390/ijms17121838
Int. J. Mol. Sci. 2016, 17, 1838 Int. J. Mol. Sci. 2016, 17, 1838
2 of 13 2 of 13
Although the mechanism of diapause has not been fully elucidated, previous molecular Although the mechanism of diapause has not been fully elucidated, previous molecular studies studies showed that the phenomenon is caused by a neuropeptide diapause hormone (DH) released showed that the phenomenon is caused by a neuropeptide diapause hormone (DH) released from the from the neurosecretory cells of the subesophageal ganglion . In the diapause stage, glucose is neurosecretory cells of the subesophageal ganglion . In the diapause stage, glucose is transformed transformed to trehalose and glycerol, which help the eggs to overcome the adverse environmental to trehalose and glycerol, which help the eggs to overcome the adverse environmental conditions, conditions, whereas in the diapause termination stage, trehalose and glycerol are degraded . whereas in the diapause termination stage, trehalose and glycerol are degraded . Some conditions can prevent egg diapause. HCl treatment has been widely used in sericulture Some conditions can prevent egg diapause. HCl treatment has been widely used in sericulture for for preventing diapause-destined eggs from entering diapause; however, the underlying preventing diapause-destined eggs from entering diapause; however, the underlying mechanism of mechanism of this effect remains unclear. Previous studies have suggested that the termination of this effect remains unclear. Previous studies have suggested that the termination of egg diapause is egg diapause is mainly caused by the presence of H+ or O−, as H2SO4 and other acids, O2, and mainly caused by the presence of H+ or O− , as H2 SO4 and other acids, O2 , and hydrogen peroxide can hydrogen peroxide can interfere with diapause [10,11]. In addition, dimethyl sulfoxide and other interfere with diapause [10,11]. In addition, dimethyl sulfoxide and other substances with high alkyl substances with high alkyl chain structures can prevent diapause . chain structures can prevent diapause . In order to investigate the mechanism of embryonic diapause in silkworm, large-scale In order to investigate the mechanism of embryonic diapause in silkworm, large-scale screenings screenings were performed, including genome-wide microarray and shotgun liquid were performed, including genome-wide microarray and shotgun liquid chromatography–tandem chromatography–tandem mass spectrometry, combined with bioinformatics, to identify mass spectrometry, combined with bioinformatics, to identify differentially expressed genes (DEGs) differentially expressed genes (DEGs) in diapause and non-diapause eggs [13–15]. Our goal was to in diapause and non-diapause eggs [13–15]. Our goal was to study the underlying mechanism of study the underlying mechanism of diapause and diapause blocking by treating diapause-destined diapause and diapause blocking by treating diapause-destined eggs with O2 or HCl. The DEGs of eggs with O2 or HCl. The DEGs of treated and non-treated eggs were analyzed by RNA-sequencing treated and non-treated eggs were analyzed by RNA-sequencing followed by bioinformatics analysis, followed by bioinformatics analysis, and the molecular impact of the two treatments is discussed. and the molecular impact of the two treatments is discussed. 2. Results 2. Results Transcription Patterns Patterns 2.1. RNA-Sequencing and General Transcription Non-treated eggs entered diapause, whereas the O O22-treated and the HCl-treated eggs avoided avoided diapause and anddeveloped developedinto into larva (Figure In order to explore the mechanism of diapause diapause larva (Figure 1). In1).order to explore the mechanism of diapause blocking blocking by O 2 or HCl, we performed RNA-sequencing and obtained 12,920,255 and 11,131,963 by O2 or HCl, we performed RNA-sequencing and obtained 12,920,255 and 11,131,963 clean reads from clean from the13,598,605 O2-treated eggs; 13,598,605 and 13,514,169 reads from the O2reads -treated eggs; and 13,514,169 clean reads from theclean HCl-treated eggs;the andHCl-treated 12,313,326 eggs;12,506,248 and 12,313,326 and 12,506,248 clean readseggs from the 1). non-treated eggs (Table 1). The read and clean reads from the non-treated (Table The read counts were converted to counts were converted to Reads Per Kilobase of exon model per Million mapped reads (RPKM) Reads Per Kilobase of exon model per Million mapped reads (RPKM) (Table S1) and showed that (Table S1) and showed the distributions of similar gene expression levels wereand similar betweeneggs the the distributions of genethat expression levels were between the treated non-treated treated S1). and The non-treated eggs S1).that Themapped proportion of reference total reads that mapped the 95.67% reference (Figure proportion of (Figure total reads to the genome rangedto from to genome ranged from 95.67% to 96.49%, and correlation values were significantly higher between 96.49%, and correlation values were significantly higher between the duplicated samples than among duplicated (Figure samplesS2). than among the treatments (Figure S2). the treatments
1. Egg incubation after different treatments. Half Figure 1. Half the the eggs eggs (left side) from each moth were the other other half half (right (right side) side) remained remained untreated untreated (control). (control). (a) (a) O O22-treated treated with O O22 or HCl, and the eggs and and (b) (b) HCl-treated HCl-treated eggs eggs and and non-treated non-treated eggs. eggs. eggs and non-treated eggs
Int. J. Mol. Sci. 2016, 17, 1838
3 of 13
Int. J. Mol. Sci. 2016, 17, 1838
3 of 13
Table 1. Statistical analysis of DEGseq data obtained from O2 -treated eggs and HCl-treated eggs. Table 1. Statistical analysis of DEGseq data obtained from O2-treated eggs and HCl-treated eggs.
Sample Sample Name Name O2-treated_1 O2 -treated_1 O2-treated_2 O2 -treated_2 HCl-treated_1 HCl-treated_1 HCl-treated_2 HCl-treated_2 Non-treated_1 Non-treated_1 Non-treated_2 Non-treated_2
Raw Reads 12,984,364 12,984,364 11,192,761 11,192,761 13,670,896 13,670,896 13,595,227 13,595,227 12,361,440 12,361,440 12,559,734 12,559,734
Clean Clean Reads Reads
Clean Clean Bases Bases
TotalMapped Mapped Total
12,920,255 12,920,255 11,131,963 11,131,963 13,598,605 13,598,605 13,514,169 13,514,169 12,313,326 12,313,326 12,506,248 12,506,248
0.65 G G 0.65 0.56 G G 0.56 0.68 G G 0.68 0.68 G G 0.68 0.62 G G 0.62 0.63 0.63 G G
12,360,505 12,360,505(95.67%) (95.67%) 10,740,908 10,740,908(96.49%) (96.49%) 13,068,581 13,068,581(96.1%) (96.1%) 12,993,548 12,993,548(96.15%) (96.15%) 11,827,559 11,827,559(96.05%) (96.05%) 11,993,516(95.9%) (95.9%) 11,993,516
Error Q20Q20 Q30 Q30 GC Content GC Error Rate Content (%) Rate (%)(%) (%)(%) (%)(%) (%) 0.01 97.9 94.07 43.13 0.01 97.9 94.07 43.13 0.01 43.38 0.01 98.84 98.84 96.07 96.07 43.38 0.01 43.44 0.01 97.92 97.92 94.13 94.13 43.44 0.01 45.03 0.01 98.64 98.64 95.87 95.87 45.03 0.01 43.15 0.01 97.91 97.91 94.12 94.12 43.15 0.01 98.16 98.16 94.66 94.66 0.01 43.543.5
Differentially ExpressedGenes Genes(DEGs) (DEGs) in in O O2-Treated and HCl-Treated Eggs 18.104.22.168. Differentially Expressed 2 -Treated and HCl-Treated Eggs To observe the gene expression patterns among O2-treated eggs, HCl-treated eggs, and To observe the gene expression patterns among O2 -treated eggs, HCl-treated eggs, and non-treated non-treated eggs, hierarchical clustering was performed based on the log10 RPKMs. The gene eggs, hierarchical clustering was performed based on the log10 RPKMs. The gene expression patterns expression patterns were similar between the HCl-treated eggs and O2-treated eggs, but differed in were similar between the HCl-treated eggs and O2 -treated eggs, but differed in the treatment groups the treatment groups from those of the non-treated eggs (Figure 2). In the O2-treated eggs, we from those of 610 the non-treated 2). In the we identified 610 DEGsin(239 2 -treated identified DEGs (239 eggs with (Figure upregulated andO371 with eggs, downregulated expression); thewith upregulated andeggs, 371 with downregulated expression); the HCl-treated we identified 656 DEGs HCl-treated we identified 656 DEGs (298 within upregulated and eggs, 358 with downregulated (298expression). with upregulated and 358 with downregulated expression). Only 194 genes showed differential Only 194 genes showed differential expression between the HCl-treated eggs and expression between the HCl-treated and O2 -treated eggs (Figure 2, Tables S2 and S3). O2-treated eggs (Figure 2, Tables S2 eggs and S3).
Figure 2. Hierarchical cluster analysis, Venn diagram, and volcano plot of differentially expressed
Figure 2. Hierarchical cluster analysis, Venn diagram, and volcano plot of differentially expressed genes genes (DEGs) in O2-treated eggs and HCl-treated eggs. (a) Hierarchical clustering of DEGs. Blue (DEGs) in O2 -treated eggs and HCl-treated eggs. (a) Hierarchical clustering of DEGs. Blue indicates indicates low expression, white indicates moderate expression, and red indicates high expression. low expression, white indicates moderate expression, and red indicates high expression. The log10 The log10 Reads Per Kilobase of exon model per Million mapped reads (RPKM) values were used to Reads Per Kilobase of exonlevel; model Million mapped (RPKM) values wereinused to estimate estimate the expression (b) per Venn diagram of genereads number and distribution O2-treated eggs, the expression level; (b) Venn diagram of gene number and distribution in O -treated eggs, HCl-treated 2 and non-treated eggs; HCl-treated eggs, and non-treated eggs; (c) DEGs between O2-treated eggs eggs, eggs; (c) DEGs between non-treated eggs;upregulated and (d) DEGs 2 -treated eggs andand (d)non-treated DEGs between HCl-treated eggs andOnon-treated eggs.and Red dots represent between eggs and non-treated eggs. Red dots represent upregulated genes and green dots genes HCl-treated and green dots downregulated genes. downregulated genes.
Int. Int.J.J.Mol. Mol.Sci. Sci.2016, 2016,17, 17,1838 1838
44of of13 13
2.3. Common DEGs between O2-Treated and HCl-Treated Eggs 2.3. Common DEGs between O2 -Treated and HCl-Treated Eggs A total of 343 DEGs (109 upregulated and 232 downregulated) were common to both A total accounting of 343 DEGsfor (109 upregulated downregulated) common to both treatments, treatments, 56.2% of DEGsand in 232 O2-treated eggs andwere 52.3% of DEGs in HCl-treated accounting 56.2% of Ontology DEGs in O(GO) -treated eggs and 52.3% of DEGs in HCl-treated eggs (Table S4). eggs (Table for S4). Gene analysis showed that the common DEGs were mainly 2 Gene Ontology (GO) analysis that the common DEGs were mainly category distributed in the the cellular binding distributed in the binding andshowed catalytic terms of the molecular function and and metabolic catalytic terms of the of molecular function category and the cellular3).and processes of the and processes the biological process category (Figure Of metabolic the 180 genes annotated biological process category (Figure Of the in 180nucleic genes annotated successfully by GO, genes successfully by GO, 25 genes were 3). involved acid binding (GO:0003676), 22 25 genes in were zinc involved in nucleic acid binding (GO:0003676), 22 genes in zinc ion binding 16 genes in ion binding (GO:0008270), 16 genes in ATP binding (GO:0005524), and 11 (GO:0008270), genes in oxidoreductase ATP binding (GO:0005524), activity (GO:0016491) (Tableand S4).11 genes in oxidoreductase activity (GO:0016491) (Table S4).
Figure 3. 3. Gene GeneOntology Ontology(GO) (GO)categories categoriesof ofcommon commondifferentially differentially expressed expressed genes genes (DEGs) (DEGs) between between Figure O22-treated -treatedeggs eggs and and HCl-treated HCl-treated eggs. used to eggs. The The Web WebGene GeneOntology OntologyAnnotation AnnotationPlot Plot(WEGO) (WEGO)was was used O plot GO annotation. to plot GO annotation.
2.4. GO GOand andKEGG KEGGEnrichment Enrichment Analysis Analysis of of DEGs DEGs 2.4. A total total of of 734 734 genes genes (79.52% (79.52% of of all DEGs) were annotated successfully successfully by by GO analysis. In In the the A O22-treated -treated eggs, enriched genes genes were associated with cofactor binding and oxidoreductase O eggs,significantly significantly enriched were associated with cofactor binding and activity in the molecular function category. In the HCl-treated eggs, no significantly enriched terms oxidoreductase activity in the molecular function category. In the HCl-treated eggs, no significantly were identified, although 63 genes were involved in were oxidoreductase activity in the molecular function enriched terms were identified, although 63 genes involved in oxidoreductase activity in the category (Figure 4, category Table S5).(Figure 4, Table S5). molecular function
Int. J. Mol. Sci. 2016, 17, 1838
Int. J. Mol. Sci. 2016, 17, 1838
5 of 13
5 of 13
Figure 4. Gene Ontology (GO) enrichment in (A) -treatedeggs eggsand and(B) (B)HCl-treated HCl-treatedeggs. eggs.The The statistical significance Figure 4. Gene Ontology (GO) enrichmentofofdifferentially differentiallyexpressed expressed genes genes (DEGs) (DEGs) in (A) O O22-treated statistical significance of of functional GOGO Slim enrichment was analyzed process,oxidoreductase oxidoreductaseactivity, activity,and and cofactor binding were significantly functional Slim enrichment was analyzedusing usingGOseq. GOseq.The Theterms terms of of oxidation-reduction oxidation-reduction process, cofactor binding were significantly enriched in O -treated eggs. of **indicate indicatethe thesignificant significant enrichment p < 0.05 during Ontology analysis. 2-treated eggs.The The symbol symbol of enrichment at pat < 0.05 during GeneGene Ontology analysis. enriched in2O
Int. J. Mol. Sci. 2016, 17, 1838
6 of 13
Int. J. Mol. Sci. 2016, 17, 1838
2.5. Validation of 2.5. Validation of RNA-seq RNA-seq Data Data by by qRT-PCR qRT-PCR The The 15 15 DEGs DEGs that that were were selected selected randomly randomly for for qRT-PCR qRT-PCR validation. validation. In In general, general, the the patterns patterns of of upregulation and downregulation were consistent with those obtained from digital gene expression upregulation and downregulation were consistent with those obtained from digital gene expression (DGE) (DGE) analysis, analysis, confirming confirming that that both both methods methods were were reliable reliable(Figure (Figure5,5,Table TableS6). S6).
Figure 5. 5. Quantitative Quantitative real-time real-time polymerase polymerase chain of differentially differentially expressed expressed genes genes (DEGs) (DEGs) Figure chain reaction reaction of in O22-treated eggs and HCl-treated eggs. The relative gene expression was normalized against the significant at at p