Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Differentially expressed genes in the head of the 2nd instar pre-molting larvae of the nm2 mutant of the silkworm, Bombyx mori

  • Pingyang Wang ,

    Contributed equally to this work with: Pingyang Wang, Simin Bi

    ‡ PYW and SMB are co-first authors on this work. FW, PZX and XJS also contributed equally to this work.

    Affiliations School of Biotechnology, Jiangsu University of Science and Technology, Zhenjiang Jiangsu, China, The Sericulture Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang Jiangsu, China

  • Simin Bi ,

    Contributed equally to this work with: Pingyang Wang, Simin Bi

    ‡ PYW and SMB are co-first authors on this work. FW, PZX and XJS also contributed equally to this work.

    Affiliations School of Biotechnology, Jiangsu University of Science and Technology, Zhenjiang Jiangsu, China, The Sericulture Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang Jiangsu, China

  • Fan Wu ,

    ‡ PYW and SMB are co-first authors on this work. FW, PZX and XJS also contributed equally to this work.

    Affiliation Industrial Crops Institute, Hubei Academy of Agricultural Sciences, Wuhan, China

  • Pingzhen Xu ,

    ‡ PYW and SMB are co-first authors on this work. FW, PZX and XJS also contributed equally to this work.

    Affiliations School of Biotechnology, Jiangsu University of Science and Technology, Zhenjiang Jiangsu, China, The Sericulture Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang Jiangsu, China

  • Xingjia Shen ,

    ‡ PYW and SMB are co-first authors on this work. FW, PZX and XJS also contributed equally to this work.

    Affiliations School of Biotechnology, Jiangsu University of Science and Technology, Zhenjiang Jiangsu, China, The Sericulture Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang Jiangsu, China

  • Qiaoling Zhao

    qlzhao302@126.com

    Affiliations School of Biotechnology, Jiangsu University of Science and Technology, Zhenjiang Jiangsu, China, The Sericulture Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang Jiangsu, China

Abstract

Molting is an important physiological process in the larval stage of Bombyx mori and is controlled by various hormones and peptides. The silkworm mutant that exhibits the phenotype of non-molting in the 2nd instar (nm2) is incapable of molting in the 2nd instar and dies after seven or more days. The ecdysone titer in the nm2 mutant is lower than that in the wildtype, and the mutant can be rescued by feeding with 20E and cholesterol. The results of positional cloning indicated that structural alteration of BmCPG10 is responsible for the phenotype of the nm2 mutant. To explore the possible relationship between BmCPG10 and the ecdysone titer as well as the genes affected by BmCPG10, digital gene expression (DGE) profile analysis was conducted in the nm2 mutant, with the wildtype strain C603 serving as the control. The results revealed 1727 differentially expressed genes, among which 651 genes were upregulated and 1076 were downregulated in nm2. BLASTGO analysis showed that these differentially expressed genes were involved in various biological processes, cellular components and molecular functions. KEGG analysis indicated an enrichment of these differentially expressed genes in 240 pathways, including metabolic pathways, pancreatic secretion, protein digestion and absorption, fat digestion and absorption and glycerolipid metabolism. To verify the accuracy of the DGE results, quantitative reverse transcription PCR (qRT-PCR) was performed, focusing on key genes in several related pathways, and the results were highly consistent with the DGE results. Our findings indicated significant differences in cuticular protein genes, ecdysone biosynthesis genes and ecdysone-related nuclear receptors genes, but no significant difference in juvenile hormone and chitin biosynthesis genes was detected. Our research findings lay the foundation for further research on the formation mechanism of the nm2 mutant.

Introduction

The silkworm, Bombyx mori, a holometabolous lepidopteran, is an experimental model for molecular entomology in the fields of genetics, physiology and biochemistry [1]. The insect epidermis lacks elasticity; as a result, insects must molt several times during the larval period to allow insect growth and development. Molting is an important physiological process in the silkworm and is controlled by prothoracicotropic hormone (PTTH), 20-hydroxyecdysone (20E), juvenile hormone (JH) and neuropeptides. When it becomes necessary to molt, ecdysone is synthesized in the prothoracic gland promoted by PTTH, after which activated 20E is combined with the heterodimeric nuclear receptor EcR-USP to promote the expression of early genes. Early-late genes are then activated by the early genes, followed by a series of late genes activated by the early-late genes, such as eclosion hormone biosynthesis genes, neuropeptide genes, cuticular protein genes and hydrolase genes, and the activating signals are amplified by the key nuclear receptor βFTZ-F1 [2]. When the 20E titer decreases in the late phase of molting, numerous pigmentation-related enzymes are expressed and transported to the epidermis. The resulting pigments participate in the formation of a new epidermis [3, 4]. Molting is mainly controlled by molting hormone and juvenile hormone. However, the pupating process requires the disappearance of juvenile hormone, and ecdysone plays a dominant role [5, 6].

During molting, cuticular proteins play key roles. A large number of cuticular proteins are synthesized during the formation of a new epidermis, and the same proteins are degraded in the apolysis of the old epidermis. It has been demonstrated that more than 1% of the genes are cuticular protein genes (CPGs) in almost all insect genomes [7]. In the Bombyx mori genome, more than 200 CPGs are predicted [7, 8] and can be classified into 12 families based on the characteristics of conservative motifs [9]. Cuticular proteins act mainly as structural proteins of the epidermis [10]. Additionally, they play roles in other processes, such as locomotion, feeding, stress resistance and immunity [1114]. Thus, CPGs play irreplaceable roles in insect development and reproduction.

Studying the functions of molting-related genes using mutants is a good way to analyze the mechanisms leading to molting. A 2nd instar non-molting (nm2) mutant that we discovered in the silkworm variety C603 provides a good model for studying the mechanisms underlying molting. The nm2 mutant appears normal in the first instar, but at the beginning of the pre-molting stage of the second instar, it is unable to molt and becomes lustrous, which lasts for 6 to 8 days, followed by death of the larva. Genetic and clone mapping analyses revealed that a deletion of 217 bp in the open reading frame of the cuticular protein gene BmCPG10 is responsible for the phenotype of the nm2 mutant [15]. In this study, the head of the nm2 mutant in the pre-molting stage was subjected to digital gene expression (DGE) profile analysis, and the wildtype C603 served as the control. The differentially expressed genes according to DGE analysis and the key genes in ecdysone, juvenile hormone and chitin biosynthesis, cuticular protein genes and ecdysone-induced nuclear receptor genes were analyzed and verified using quantitative reverse transcription polymerase chain reaction (qRT-PCR). The possible mechanism underlying the development of the nm2 mutant phenotype is explained based on the results of our DGE analysis and previous studies.

Materials and methods

Silkworm rearing and tissue isolation

The nm2 mutant and wildtype C603 silkworm strains were supplied by the Sericulture Research Institute, Chinese Academy of Agricultural Sciences (Zhenjiang, China). The larvae were fed fresh mulberry leaves under standard conditions: 25 ±2°C temperature, 12-h light/12-h dark photoperiod and 65 ±5% relative humidity. The 2nd instar pre-molting larvae of the wildtype and the nm2 mutant were dissected and subsequently stored at −80°C.

RNA extraction

Total RNA was extracted from the head of the 2nd instar pre-molting larvae using RNAiso Plus (TaKaRa, China), according to the instructions of the manufacturer. After treatment with DNase I (TaKaRa, China) to degrade DNA contamination, the concentration of total RNA was determined using a NANODROP1000 microspectrophotometer (Thermo, USA). Total RNA quality was determined based on the 260/280 absorbance ratio and electrophoresis and then stored at −80°C.

Digital gene expression profile

The treated mRNA was enriched using oligo (dT) magnetic beads. By mixing with fragmentation buffer, the mRNA was fragmented into short fragments. Then, the first strand of cDNA was synthesized using random hexamer primers. Buffer, dNTPs, RNase H and DNA polymerase I were subsequently added to synthesize the second strand. Next, the double-stranded cDNA was purified with magnetic beads. End repair and 3’-end single A (adenine) nucleotide addition were then performed. Finally, sequencing adaptors were ligated to the fragments, and the fragments were enriched via PCR amplification. During the quality control (QC) step, an Agilent 2100 Bioanaylzer and an ABI StepOnePlus Real-Time PCR System were used to qualify and quantify the sample library. The library products were then ready for sequencing via the Illumina HiSeqTM 2000 platform at the Beijing Genomics Institute (BGI).

DGE data analysis

The primary sequencing data produced by the Illumina HiSeqTM 2000 platform, referred to as raw reads, were subjected to quality control (QC) procedures to determine whether a resequencing step was needed. After QC, the raw reads were filtered into clean reads, which were aligned to reference sequences from the silkworm genome database, SilkDB (http://silkworm.swu.edu.cn/silkdb/). QC of the alignment was performed to determine whether resequencing was needed. The alignment data were utilized to calculate the read distribution and mapping ratio based on the reference genes. If the alignment results passed QC, downstream analyses, including gene expression and deep analyses, proceeded based on the observed gene expression (e.g., principal component analysis/correlation/screening of differentially expressed genes). Deep analyses based on DGE were also performed, including Gene Ontology (GO) enrichment analysis and KEGG pathway enrichment analysis.

RSEM (RNA-Seq by Expectation Maximization) [16] is a quantification tool that computes maximum likelihood abundance estimates using the expectation maximization (EM) algorithm in 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. Expression levels were calculated according to the FPKM (Fragments Per Kb per Million fragments) method. The differentially expressed genes were filtered by referring to the digital gene expression profiles reported by Audic et al. [17]. P-values corresponding to differential gene expression were calibrated via multiple hypothesis testing, and the P-value threshold was determined by controlling the FDR (false discovery rate) [18]. The significance of gene expression differences was judged using |Log2 (nm2/wildtype)| ≥ 1 and an FDR ≤ 0.001 as the default thresholds.

Quantitative reverse transcriptase PCR

The heads of nm2 mutant and wildtype C603 larvae in the pre-molting stage were collected. Total RNA was extracted using RNAiso Plus (TaKaRa, China), according to the instructions of the manufacturer. After treatment with DNase I (TaKaRa, China), cDNA was synthesized using Reverse Transcriptase M-MLV (RNase H-) kit (TaKaRa, China), according to the instructions of the manufacturer, and was then diluted to 50 ng/μL to serve as the template for qRT-PCR. The 20-μL reaction system included 1 μL of primers (10 μmol/L, S1 Table), 1 μL of cDNA, 10 μL of 2×SYBR® Premix Ex Taq™ (Tli RNaseH Plus) (TaKaRa, China) and 8 μL of ddH2O. After transient centrifugation, qRT-PCR was performed using a LightCycle 96 Real-time PCR System (Roche, Switzerland) with the following reaction program: three-step reaction protocol consisting of 45 cycles of 95°C for 10 s, 58°C for 10 s and 72°C for 10 s, after a 10-min step of pre-degeneration, followed by melting. The quality of the qRT-PCR product was tested through melting curve analysis, and relative expression was calculated using the 2–ΔΔCt method [19], with the average of three house-keeping genes, RPL3 (BGIBMGA013567), GAPDH (BGIBMGA007490) and BmActin3 (BGIBMGA005576), serving as the reference. The qRT-PCR results were compared with the DGE results to verify the accuracy of the DGE data.

Results

Analysis of DGE libraries

From the DGE libraries, we acquired 12,912,619 and 12,900,347 clean reads from the 13,129,213 and 13,129,324 raw reads after sifting of 216,594 and 216,594 adapter reads, 21,750 and 23,019 reads in which unknown bases were more than 10% and some low quality reads that the percentage of low quality bases was over 50% in a read in the wildtype C603 and the nm2 mutant, respectively. The ratio of clean reads to raw reads was 98.35% in the wildtype and 98.26% in the nm2 mutant (Table 1). The GC% for nm2 (50.88%) was slightly higher than that for the wildtype (49.65%). The clean reads from the two libraries contained unknown N bases but at proportions as low as 0.0034% and 0.0036%, respectively (Table 1).

thumbnail
Table 1. General information on reads from the two libraries.

https://doi.org/10.1371/journal.pone.0180160.t001

Evaluation of the reads showed that the base distribution of the reads in the two libraries was decentralized in the first 10 bases, which was normal because of the sequencing adaptors. Thereafter, the four type of bases and unknown N bases became uniform, thus meeting the requirements for follow-up analysis (Fig 1A). The analysis of the quality distribution showed a high quality for bases from the 6th to the last base and even for the first 6 bases (Fig 1B). In addition, the results of Q20 and Q30 analysis all met the requirements for follow-up analysis (Table 1).

thumbnail
Fig 1. Analysis of read data.

(A) Percentage composition of bases in the reads. The four types of base and unknown bases (N) were uniform from the 10th base. (B) Distribution of quality. High quality was observed from the 6th base to the last and even for the first 6 bases.

https://doi.org/10.1371/journal.pone.0180160.g001

Read mapping analysis

By aligning the clean reads to the reference genome database SilkDB (http://silkworm.genomics.org.cn/), we found that 81.61% and 81.13% of the clean reads could be aligned to the genome in the wildtype and the nm2 mutant, respectively, among which 57.11% and 57.64% were completely aligned to the genome, and 24.5% and 23.49% were partially aligned to genome. In addition, 65.46% and 64.56% of the clean reads could be mapped to unique sequences, while 16.15% and 16.56% mapped to sequences at multiple positions in the wildtype and the nm2 mutant, respectively. The clean reads were also aligned to 14623 known genes [20], and the percentages of total mapped reads, perfectly mapped reads, partially mapped reads, uniquely mapped reads and multiply mapped reads in the wildtype were 55.48%, 38.66%, 16.82%, 50.81% and 4.67%, respectively, but in the nm2 mutant, the corresponding values were 60.12%, 42.58%, 17.54%, 53.28% and 6.83% (Table 1). Quality assessments, including the unique gene mapping ratio and the genome mapping ratio, also met the requirements for follow-up analysis (Table 1).

Reads random analysis showed that most of the clean reads could be mapped to the middle of genes and fewer to the 5 'and 3' ends (Fig 2A). Sequence saturation analysis showed that in the range of 0 to 70%, saturation increased in nearly linearly with an increase in the number of reads. Only approximately 1 Mega reads reached a saturation of 75%, while the clean reads obtained from DGE analysis reached 13 Mega, and the sequencing depth met the requirements for follow-up analysis (Fig 2B).

thumbnail
Fig 2. Analysis of mapping data.

(A) Reads random from the wildtype and the nm2 mutant. Most of the clean reads could be mapped to the middle of genes. (B) Sequence saturation of the wildtype and the nm2 mutant. The saturation increased almost linearly with the increase in the number of reads in the range of 0 to 70%.

https://doi.org/10.1371/journal.pone.0180160.g002

Analysis of differentially expressed genes

In the two libraries, 11,856 genes were found to be expressed in the wildtype (81.08% of all 14,623 genes) and 11,751 in the nm2 mutant (80.36%), among which 11,231 genes were expressed in both the wildtype and the nm2 mutant. In addition, 625 genes were expressed specifically in the wildtype, whereas 520 were expressed specifically in the nm2 mutant (S2 Table).

The scatter plots presented the distribution of differentially expressed genes in screening threshold dimensions (Fig 3). We observed 651 upregulated genes in the nm2 mutant and 1076 downregulated genes, according to the default thresholds of an FDR≤0.001 and an FPKM ratio = |log2 (nm2/wildtype)| ≥1 (S3 Table).

thumbnail
Fig 3. Scatter plots of all expressed genes in each pairwise.

Blue blots mean down-regulation genes, orange blots mean up-regulation genes and brown blots mean non-regulation genes. The screening threshold was on top legend.

https://doi.org/10.1371/journal.pone.0180160.g003

The GO enrichment analysis showed that these 1727 differentially expressed genes (651 upregulated genes and 1076 downregulated genes) were mainly involved in biological processes such as metabolic processes, cellular processes, single-organism processes, localization, the response to stimulus and biological regulation; cellular components such as the cell, cell parts, membranes and organelles; and molecular functions such as catalytic activity, binding and transporter activity (Fig 4, S4 Table).

thumbnail
Fig 4. GO enrichment of biological processes, cellular components and molecular functions.

The differentially expressed genes were mainly related to biological processes such as metabolic processes, cellular processes, single-organism processes, localization, the response to stimulus and biological regulation; cellular components such as the cell, cell parts, membranes and organelles; and molecular functions such as catalytic activity, binding and transporter activity.

https://doi.org/10.1371/journal.pone.0180160.g004

KEGG cluster analysis was performed for these 1727 differentially expressed genes using the KEGG database (http://www.genome.jp/kegg/kegg2). Approximately 1130 genes were mapped to 240 pathways, among which 273 genes (159 upregulated and 114 downregulated) were mapped to the first of the top 5 pathway categories: metabolic pathways. The other categories in the top 5 were pancreatic secretion (42 upregulated genes among 68 total genes), protein digestion and absorption (28 upregulated genes among 49 total genes), fat digestion and absorption (26 upregulated genes among 46 total genes) and glycerolipid metabolism (27 upregulated genes among 45 total genes) (Fig 5, S4 Table).

thumbnail
Fig 5. Statistics of top 20 pathway enrichment of differentially expressed genes in each pairwise.

Rich factor is the ratio of differentially expressed gene numbers annotated in this pathway term to all gene numbers annotated in this pathway term. Greater rich factor means greater intensiveness. Q-value is corrected P-value ranging from 0~1, and less Q-value means greater intensiveness.

https://doi.org/10.1371/journal.pone.0180160.g005

Verifying the accuracy of the DGE data

To verify the accuracy of the DGE data, we analyzed the expression of several housekeeping genes as well as several target genes that had been investigated in previous studies. Additionally, some key genes, including partial differentially expressed genes based on DGE, ecdysone, juvenile hormone and chitin biosynthesis and metabolism genes, and several nuclear receptor genes, were verified using qRT-PCR (Table 2, S1 Table).

The fold-changes of three housekeeping genes, BmActin3 (Actin), GAPDH (glyceraldehyde-3-phosphate dehydrogenase) and RPL3 (ribosomal protein L3), between the wildtype and the nm2 mutant were approximately equal to one, which preliminarily verified the accuracy of the DGE data. BmCPG10 [15] and BmCP-like [21] were upregulated in the nm2 mutant, with the fold-changes of 2.37 and 86.16, respectively. These results were consistent with previous studies, which indicated that the DGE data were accurate and reliable.

The results of qRT-PCR for 29 genes among 39 total genes were completely consistent with the DGE results. The results for 6 of the remaining 10 genes were completely inconsistent with the DGE results. The remaining 4 genes did not show a significant difference between the DGE and qRT-PCR results, which are indicated with # in Table 2. Thus, 33 genes were consistent between the DGE and qRT-PCR analyses, and the accuracy rate reached 84.62%, demonstrating that the DGE results were reliable.

Discussion

Differentially expressed genes based on DGE analysis

A large number of biosynthesis and decomposition processes occur during molting, such as carbohydrate and lipid metabolism, protein metabolism and nucleic acid metabolism [22, 23], and many substances are decomposed in apolysis, while many substances are synthesized in the formation of a new epidermis [24, 25]. Because of the non-molting phenotype of the nm2 mutant, the synthesis and metabolism of proteins, lipids, glucose and nucleic acids are greatly altered. We verified 15 genes selected from the DGE results using qRT-PCR, 12 of which were consistent with the DGE analysis. Four of these genes encode glycoside hydrolase; 2 encode cuticular proteins; 2 encode nuclear receptors; 3 encode glucose-methanol-choline oxidoreductase; and one encodes peptidase.

Hydrolase is a gene superfamily including protein hydrolase, lipid hydrolase, glycoside hydrolase and nuclear acid hydrolase [2628]. Because nm2 cannot molt in the pre-molting stage of the 2nd instar, the expression of hydrolases differed between the nm2 mutant and the wildtype. The epidermis is one of the most important and largest tissues in insects, with the main function of protecting the insect from danger. The CPGs are a large family, with more than 1% of the genes in the insect genome encoding cuticular proteins [7]. Because the main phenotype of the nm2 mutant is that it cannot molt in the pre-molting stage of the 2nd instar, many CPGs might be differentially expressed in this mutant. Glucose is a molecule of prime importance, functioning as an energy substance, as a substrate for synthesizing other substances, as a component of glycoproteins and in recognition processes [2931]. Glucose-methanol-choline oxidoreductase (GMC oxidoreductase) is homologous to D. melanogaster glucose dehydrogenase, E. coli choline dehydrogenase, A. niger glucose oxidase, H. polymorpha methanol oxidase and B. sterolicum cholesterol oxidase, which all contain a canonical ADP-binding beta alpha beta-fold close to their amino termini [32, 33]. Additionally, Chitin can bind cuticular proteins to form the cuticle, and chitin is synthesized from glucose [29]. In the nm2 mutant, the ecdysone titer was lower than in the wildtype, and the mutant could be rescued by feeding with cholesterol. The lack of cholesterol might be related to cholesterol oxidase, which belongs to the GMC oxidoreductase family [33]. Chitin can combine with cuticular proteins to participate in the formation of the epidermis, and the nm2 mutant accordingly cannot molt. In addition, the nm2 mutant ingests smaller amounts of mulberry leaves and shows almost no development [15], which might be related to abnormalities in glucose metabolism.

Differentially expressed genes in the juvenile hormone synthesis pathway

Juvenile hormone, a physiologically active substance produced in the corpora allata, has a variety of natural active forms. It plays roles in many physiological process, such as the maintenance of larval morphological traits, promotion of gonad maturity, adult diapause and pheromone production [34, 35]. Juvenile hormone always works together with ecdysone to regulate insect growth and development. The results of the DGE and qRT-PCR analyses showed no obviously differentially expressed genes in the juvenile hormone biosynthesis pathway between the wildtype and the nm2 mutant. The key genes JHDH and JHDK were only slightly downregulated, without any significant difference (Fig 6A). Based on the above results, we speculated that the biosynthesis of juvenile hormone was not altered in the nm2 mutant. The content of juvenile hormone also did not present a significant difference, which could be preliminarily verified by the finding that the characteristics of the nm2 mutant exhibited no improvement after feeding with juvenile hormone (data not shown).

thumbnail
Fig 6.

Synthesis pathways of juvenile hormone A and ecdysone B. Genes with green background represent their downregulation and red for upregulation in nm2 mutant, and yellow background represent no significant differences between C603 and nm2. Three key genes in the juvenile hormone biosynthesis pathway, JHDH, JHDK and JHAMT, were without any significant difference in the nm2 mutant. In the ecdysone synthesis pathway, neverland, spook and sad were upregulated in the nm2 mutant, and CYP314A1 and CYP18A1 were downregulated. There were no significant difference for Phm and Dib.

https://doi.org/10.1371/journal.pone.0180160.g006

Differentially expressed genes in the ecdysone synthesis pathway

Ecdysone is one of the most important hormones in metamorphosis during insect development. During the larval stages, ecdysone promotes the molting process in every stage, together with juvenile hormone. In the metamorphosis stage, ecdysone individually promotes metamorphosis [36]. The synthesis of ecdysone from the substrate cholesterol involves a series of halloween genes, and ecdysone then functions as the activated form 20-hydroxyecdysone (20E) [37]. It has been demonstrated that the 20E titer in the nm2 mutant is significantly lower than in the wildtype. Additionally, the mutant characteristics can be rescued by feeding with 20E, cholesterol or 7-dehydrocholesterol (7dC) [15]. These results confirmed the lack of 20E in the nm2 mutant, which might be caused by a lack of cholesterol.

In the ecdysone synthesis pathway, neverland, spook and sad were found to be upregulated in the nm2 mutant and CYP314A1 and CYP18A1 to be downregulated (Fig 6B). Because of the lack of 20E in the nm2 mutant, early genes such as neverland (fold = 11.64) and spook (fold = 2.84) were observed to be upregulated in an attempt to increase the production of 20E, which is controlled by feedback regulation. Other genes in the ecdysone biosynthesis pathway also exhibited a slight increase expression, with the sad gene being particularly upregulated (5.47-fold). However, the upregulation of these genes could not rescue the absence of 20E because of the lack of cholesterol. Thus, very little ecdysone was transformed to 20E, and interestingly, we found that the CYP314A1 gene (fold = 0.05), which plays a role in the conversion from ecdysone to 20E [37, 38], was markedly downregulated in the nm2 mutant. The CYP18A1 gene was also downregulated (0.43-fold) because of the lack of 20E (shown in Table 2 and Fig 6B).

Differentially expressed genes in the chitin synthesis pathway

Chitin, whose main component is glycosaminoglycan [39], forms the cuticle of the insect epidermis, trachea and peritrophic membrane of digestive tube, by binding proteins to protect the insect from damage [40]. The synthesis of chitin begins with glucose, which is ultimately converted to UDP-N-acetylglucosamine [29]. Chitin synthesis is a complex process, and many enzymes function in this pathway [41]. We chose three key enzymes for analysis: glucose-6-phosphate isomerase (G6PI), glutamine: fructose-6-phosphate-aminotransferase (GFPA) and UDP-N-acetylglucosamine pyrophosphorylase (UNAP), of which UNAP is the rate-limiting enzyme [41]. In the chitin metabolism pathway, after chitin is degraded to oligomers via chitinase, β-N-acetyl-glucosaminidase (NAG) converts the oligomers to monomers [42]. The results of DGE and qRT-PCR showed no obviously differences between nm2 and the wildtype, which indicated that the synthesis of chitin occurred normally and that the content of chitin might remain at a normal level in the nm2 mutant. Chitin can bind cuticular proteins to participate in epidermis construction following the appropriate signal, which is consistent with the finding that the nm2 mutant can be rescued by feeding 20E.

Differentially expressed cuticle protein genes

The epidermis is an important insect organ that can protect the insect from harm. The epidermis is prerequisite for growth, reproduction and adaptation to the complex and changeable living environment experienced by insects [7, 43]. In the insect genome, there are many CPGs, encoding a variety of cuticular proteins, and there are at least 200 CPGs in the silkworm [7, 8]. The analysis of 208 CPGs between nm2 and the wildtype showed that 59 CPGs were not expressed in the wildtype and nm2; 81 CPGs were upregulated in nm2; and 68 CPGs downregulated (S5 Table). After removing 61 genes exhibiting low expression, with an FPKM<5, in both the wildtype and nm2, 56 genes were upregulated, and 32 genes were downregulated in nm2. Considering genes showing a fold-change ≥2, ≥5 or ≥10 led to the identification of 48, 33 and 21 upregulated genes and 30, 26 and 23 downregulated genes in nm2, respectively (see Fig 7). Based on the CPG statistics, we identified many CPGs that were expressed differentially between the wildtype and nm2, including many genes exhibiting very large differences. Many CPGs have ecdysone binding sites, and their expression is regulated by the ecdysone titer [4447]. The ecdysone titer was lower in the nm2 mutant than in the wildtype, and CPGs controlled by ecdysone were therefore expressed differentially.

thumbnail
Fig 7. Expression of cuticle protein genes.

Many CPGs were expressed differentially between nm2 and C603, with more than 20 differentially expressed genes showing a fold-change of 10 or greater.

https://doi.org/10.1371/journal.pone.0180160.g007

Differentially expressed ecdysone-induced nuclear receptor genes

The nuclear receptor gene superfamily consists of transcription factors that participate in many biological processes, including molting metamorphosis, embryonic development, cell differentiation, reproduction and other physiological processes [9, 48], and many nuclear receptor genes are involved in the ecdysone-induced signaling pathway. Ashburner [49] first established this signaling pathway in 1973; after verification by a number of scholars, Bonneton et al. [50] amend the pathway in 2008. In this pathway, early genes such as E75, Broad, E74, E93 and HR39 are first activated after the binding of 20E and EcR-Usp heterologous dimers, followed by the activation of early-late genes such as HR3, HR4, HR38 and E78 and, finally, the activation of late genes, by amplifying the activation signal through βFTZ-F1. By analyzing silkworm nuclear receptor genes [51], we found that the key gene βFTZ-F1 was significantly upregulated (5.04-fold); the early-late genes HR3 and HR4 were significantly downregulated; and no significant difference was observed for other nuclear receptor genes, with usp being slightly upregulated, whereas the others were slightly downregulated (Table 2).

βFTZ-F1 is expressed in late embryonic development and in the larval and pupal stages at each developmental transition in insects and is controlled by 20E. βFTZ-F1 expression significantly increases when the 20E titer is low [52]. In addition, βFTZ-F1 promotes ecdysone biosynthesis during molting [53]. In the nm2 mutant, the 20E titer is much lower than in the wildtype, and the mutant can be rescued by feeding with 20E. A low 20E titer causes an increase in βFTZ-F1, which is consistent with the significant upregulation of βFTZ-F1 observed in the nm2 mutant. The functions of HR3 and HR4, which are early-late genes in the ecdysone-induced signaling pathway, are similar in the activation of relatively late genes. HR3 is closely related to insect molting [54], and HR4 is involved in ovum formation [55]. The nm2 mutant is non-molting in the 2nd instar, which might be related to downstream genes regulated by HR3 and HR4. The DGE results showed many CPGs were differentially expressed between the wildtype and nm2, and these CPGs might be controlled by HR3 and HR4 and be differentially expressed along with the downregulation of HR3 and HR4.

BmCPG10 is upregulated in the nm2 mutant

An approximately 217 bp deletion in the open reading frame (ORF) of the BmCPG10 gene was found to be the key cause of the nm2 mutant phenotype. The results of a semi-quantitative reverse-transcription polymerase chain reaction (RT-PCR) analysis showed BmCPG10 to be highly expressed in the wildtype C603 and not expressed in the nm2 mutant [15]. However, our DGE results indicated that BmCPG10 was upregulated in the nm2 mutant (see Table 2). After analyzing BmCPG10 in the nm2 mutant, we found that the RT-PCR primers targeted the mutant region where the 217 bp deleted sequence was located, which resulted in the false observation that BmCPG10 was not expressed in the nm2 mutant because the 217 bp deleted sequence did not exist in the nm2 mutant, leading to the absence of an RT-PCR product. Hence, we used a pair of primers targeting the normal region of the gene for qRT-PCR, and the results showed that BmCPG10 was indeed upregulated in the nm2 mutant (Table 2).

Possible formation mechanism of the nm2 mutant

According to previous research, cuticular proteins are thought to be regulated by the 20E titer, as late genes in the ecdysone-induced signaling pathway [56]. When the 20E titer reaches a threshold value, 20E combines with the ecdysone receptor and promotes a series of genes, including a series of CPGs, to complete the molting process [57, 58]. A large number of genes in the ecdysone-induced signaling pathway were expressed abnormally because of the low titer of 20E in the nm2 mutant. The nobo gene is essential for ecdysteroid biosynthesis via regulating the behavior of cholesterol [59], and this gene was found to be slightly upregulated. Most genes in the ecdysone synthesis pathway were upregulated, while the CYP314A1 gene, catalyzing the transformation from ecdysone to 20E, was markedly downregulated. The nuclear receptor genes in the ecdysone-induced signaling pathway were also expressed differentially, with most of these receptor genes being downregulated, except for FTZ-F1 and USP.

The previous research confirmed that the cuticular protein gene BmCPG10 is responsible for the phenotype of the nm2 silkworm mutant. The 20E titer in nm2 is lower than that in the wildtype, and the nm2 mutant can be rescued by feeding with 20E [15], indicating that the 20E titer is regulated by BmCPG10 to some extent. However, CPGs always act as effectors in the ecdysone-induced signaling pathway [2, 56]. Based on these findings, we conjecture that BmCPG10 might be endowed with new functions. A possible mechanism underlying the development of the nm2 mutant phenotype is suggested in Fig 8.

thumbnail
Fig 8. The ecdysone-induced signaling pathway.

The 20E titer in the nm2 mutant was lower than in the wildtype, leading to upregulation of ecdysone biosynthesis genes and downregulation of CYP314A1, which takes part in the conversion from ecdysone to 20E. Two nuclear receptor genes, HR3 and HR4, were markedly downregulated, whereas the key nuclear receptor gene βFTZ-F1 was significantly upregulated. BmCPG10 might act as an EPDFP monitor to regulate the molting process by controlling the biosynthesis of ecdysone.

https://doi.org/10.1371/journal.pone.0180160.g008

The cuticle of insects does not exhibit elasticity and cannot grow with the growth and development of larvae. Therefore, molting must occur every time the cuticle is no longer sufficiently large to accommodate the growth and development of larvae [60]. In the molting process, molting hormone and juvenile hormone play key roles. Additionally, illumination, temperature, nutrition and the expansion pressure and distraction forces of polypides (EPDFP) can promote the process of molting [61, 62]. However, the possible mechanism promoting molting is still unclear. We speculate that BmCPG10 is one type of molting promoter that may monitor EPDFP. When EPDFP occurs at a low level, BmCPG10 is highly expressed, but when EPDFP reaches a threshold value, the expression of BmCPG10 declines to a low level; during this process, the molting process is promoted. Furthermore, BmCPG10 is necessary for cholesterol generation via ecdysteroid synthesis or ecdysteroid transportation. In the nm2 mutant, the function of BmCPG10 is lost, and it cannot play the role monitoring EPDFP. Hence, when the silkworm larvae meet the requirements for molting in the 2nd instar, BmCPG10 cannot respond to EPDFP, and the expression of BmCPG10 remains at a high level (3.56-fold); thus, the molting process is not promoted. In the 1st instar, another CPG might act as the EPDFP monitor to promote the first molt, explaining why the phenotype of non-molting in the 2nd instar occurs in the nm2 mutant.

Conclusion

Several ecdysone synthesis genes were found to be significantly upregulated in nm2 as a result of a lack of 20E, while CYP314A1 is significantly downregulated in nm2 due to a lack of ecdysone. As the effectors of the ecdysone-induced signaling pathway, many CPGs, including the mutated BmCPG10 gene of the nm2 mutant, are expressed differentially between the wildtype and nm2. Three ecdysone-induced nuclear receptor genes (FTZ-F1, HR3 and HR4) are greatly affected by the lack of 20E.

Supporting information

S1 Table. Primers for qRT-PCR, functional description and qRT-PCR and FPKM results.

https://doi.org/10.1371/journal.pone.0180160.s001

(XLSX)

S2 Table. FPKMs of all genes in nm2 and C603.

https://doi.org/10.1371/journal.pone.0180160.s002

(XLSX)

S3 Table. Differentially expressed genes between nm2 and C603.

https://doi.org/10.1371/journal.pone.0180160.s003

(XLSX)

S4 Table. GO enrichment and KEGG enrichment analyses of differentially expressed genes.

https://doi.org/10.1371/journal.pone.0180160.s004

(XLSX)

S5 Table. Expression of cuticular protein genes.

https://doi.org/10.1371/journal.pone.0180160.s005

(XLSX)

Acknowledgments

This work was supported by grants from the National Natural Science Foundation of China (No. 31372378).

This work was also supported by the Graduate Student Innovation Program of Jiangsu University of Science and Technology (No. YCX16B-02).

We are very grateful to Minglei Wang, Fen He, Weiyang Wei, Yuan Zhang, and Chenjie Yang for their help and assistance in the experiments.

Author Contributions

  1. Conceptualization: QLZ PYW.
  2. Data curation: PYW FW.
  3. Formal analysis: PYW SMB.
  4. Funding acquisition: QLZ.
  5. Investigation: PYW SMB.
  6. Methodology: PYW QLZ.
  7. Project administration: PYW SMB FW.
  8. Resources: QLZ PZX XJS.
  9. Supervision: QLZ PYW SMB FW.
  10. Validation: PYW SMB FW.
  11. Visualization: PYW QLZ.
  12. Writing – original draft: PYW QLZ.
  13. Writing – review & editing: PYW.

References

  1. 1. Goldsmith MR, Shimada T, Abe H. The genetics and genomics of the silkworm, Bombyx mori. Annu Rev Entomol. 2005; 50: 71–100. pmid:15355234.
  2. 2. Christiaens O. Key role and diversity of EcR/USP and other nuclear receptors in selected Arthropoda species. Ghent University. 2013.
  3. 3. Hiruma K, Carter MS, Riddiford LM. Characterization of the dopa decarboxylase gene of Manduca sexta and its suppression by 20-hydroxyecdysone. Dev Biol. 1995; 169(1): 195–209. pmid:7750638.
  4. 4. Cerenius L, Soderhall K. The prophenoloxidase-activating system in invertebrates. Immunol Rev. 2004; 198: 116–126. pmid:15199959.
  5. 5. Kinjoh T, Kaneko Y, Itoyama K, Mita K, Hiruma K, Shinoda T. Control of juvenile hormone biosynthesis in Bombyx mori: cloning of the enzymes in the mevalonate pathway and assessment of their developmental expression in the corpora allata. Insect Biochem Mol Biol. 2007; 37(8): 808–818. pmid:17628279.
  6. 6. Muramatsu D, Kinjoh T, Shinoda T, Hiruma K. The role of 20-hydroxyecdysone and juvenile hormone in pupal commitment of the epidermis of the silkworm, Bombyx mori. Mech Dev. 2008; 125(5–6): 411–420. pmid:18331786.
  7. 7. Futahashi R, Okamoto S, Kawasaki H, Zhong YS, Iwanaga M, Mita K, et al. Genome-wide identification of cuticular protein genes in the silkworm, Bombyx mori. Insect Biochem Mol Biol. 2008;38(12):1138–1146. pmid:19280704.
  8. 8. Liang J, Zhang L, Xiang Z, He N. Expression profile of cuticular genes of silkworm, Bombyx mori. BMC Genomics. 2010; 11(1): 1702–1716. pmid:20226095.
  9. 9. Kingjones K, Thummel CS. Nuclear receptors—a perspective from Drosophila. Nature Reviews Genetics. 2005; 6(4): 311–323. pmid:15803199.
  10. 10. Soares MPM, Silvatorres FA, Eliasneto M, Nunes FMF, Simões ZLP, Bitondi MMG. Ecdysteroid-dependent expression of the tweedle and peroxidase genes during adult cuticle formation in the honey bee, Apis mellifera. Plos One. 2011; 6(5): e20513. pmid:21655217.
  11. 11. Jasrapuria S, Specht CA, Kramer KJ, Beeman RW, Muthukrishnan S. Gene families of cuticular proteins analogous to peritrophins (CPAPs) in Tribolium castaneum have diverse functions. Plos One. 2012; 7(11): e49844. pmid:23185457.
  12. 12. Asano T, Taoka M, Shinkawa T. Identification of a cuticle protein with unique repeated motifs in the silkworm, Bombyx mori. Insect Biochemistry & Molecular Biology. 2013; 43(4): 344–351. pmid:23376333.
  13. 13. Qiao L, Xiong G, Wang RX, He SZ, Chen J, Tong XL, et al. Mutation of a cuticular protein, BmorCPR2, alters larval body shape and adaptability in silkworm, Bombyx mori. Genetics. 2014; 196(4): 1103–1115. pmid:24514903.
  14. 14. Teets NM, Denlinger DL. Surviving in a frozen desert: environmental stress physiology of terrestrial Antarctic arthropods. Journal of Experimental Biology. 2014; 217(217): 84–93. pmid:24353207.
  15. 15. Wu F, Wang P, Zhao Q, Kang L, Xia D, Qiu Z, et al. Mutation of a cuticle protein gene, BmCPG10, is responsible for silkworm non-moulting in the 2nd instar mutant. PLoS One. 2016; 11(4): e0153549. pmid:27096617.
  16. 16. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011; 12(1): 93–9. pmid:21816040.
  17. 17. Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Research. 1997; 7(10): 986–995. pmid:9331369.
  18. 18. Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Annals of Statistics. 2001; 29(4): 1165–1188. PMID:18298808.
  19. 19. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001; 25(4): 402–408. pmid:11846609.
  20. 20. Genome IS. The genome of a lepidopteran model insect, the silkworm Bombyx mori. Insect biochemistry and molecular biology. 2008; 38(12): 1036–1045. pmid:19121390.
  21. 21. Wu F, Kang L, Wang P, Zhao Q. The expression analysis of cysteine proteinase-like protein in wild-type and nm2 mutant silkworm (Lepidoptera: Bombyx mori). Gene. 2016; 586(1): 170–175. pmid:27080953.
  22. 22. Ashburner M. Chromosomal action of ecdysone. Nature. 1980; 285(5765): 435–436. pmid:6250033.
  23. 23. Catalan RE, Martinez AM, Aragones MD. In vitro effect of ecdysterone on protein kinase activity. Comparative Biochemistry & Physiology Part B Comparative Biochemistry. 1982; 71(2): 301–303.
  24. 24. Ote M, Mita K, Kawasaki H, Daimon T, Kobayashi M, Shimada T. Identification of molting fluid carboxypeptidase A (MF-CPA) in Bombyx mori. Comparative Biochemistry & Physiology Part B Biochemistry & Molecular Biology. 2005; 141(3): 314–322. pmid:15936966.
  25. 25. Sui YP, Liu XB, Chai LQ, Wang JX, Zhao XF. Characterization and influences of classical insect hormones on the expression profiles of a molting carboxypeptidase A from the cotton bollworm (Helicoverpa armigera). Insect Molecular Biology. 2009; 18(3): 353–363. pmid:19523067.
  26. 26. Lammens W, Le Roy K, Schroeven L, Van Laere A, Rabijns A, Van den Ende W. Structural insights into glycoside hydrolase family 32 and 68 enzymes: functional implications. J Exp Bot. 2009; 60(3): 727–740. pmid:19129163. 19129163.
  27. 27. Bacik JP, Whitworth GE, Stubbs KA, Vocadlo DJ, Mark BL. Active site plasticity within the glycoside hydrolase NagZ underlies a dynamic mechanism of substrate distortion. Chem Biol. 2012; 19(11): 1471–1482. pmid:23177201.
  28. 28. Graebin NG, Schoffer Jda N, Andrades D, Hertz PF, Ayub MA, Rodrigues RC. Immobilization of glycoside hydrolase families GH1, GH13, and GH70: state of the art and perspectives. Molecules. 2016; 21(8): 1074. pmid:27548117.
  29. 29. Candy DJ, Kilby BA. Studies on chitin synthesis in the desert locust. Journal of Experimental Biology. 1962; 39(1): 129–140.
  30. 30. Sharon N, Lis H. Carbohydrates in cell recognition. Scientific American. 1993; 268(1): 82–89. pmid:7678182
  31. 31. Lairson LL, Henrissat B, Davies GJ, Withers SG. Glycosyltransferases: structures, functions, and mechanisms. Annual Review of Biochemistry. 2008; 77(77): 521–555. pmid:18518825.
  32. 32. Cavener DR. GMC oxidoreductases. A newly defined family of homologous proteins with diverse catalytic activities. Journal of Molecular Biology. 1992; 223(3): 811–814. pmid:1542121.
  33. 33. Li J, Vrielink A, Brick P, Blow DM. Crystal structure of cholesterol oxidase complexed with a steroid substrate: implications for flavin adenine dinucleotide dependent alcohol oxidases. Biochemistry. 1993; 32(43): 11507–11515. pmid:8218217.
  34. 34. Mauchamp B, Darrouzet E, Malosse C, Couillaud F. 4′-OH-JH-III: an additional hydroxylated juvenile hormone produced by locust corpora allata in vitro. Insect Biochemistry & Molecular Biology. 1999; 29(6): 475–480.
  35. 35. Steiner B, Pfisterwilhelm R, Grossniklausbürgin C, Rembold H, Treiblmayr K, Lanzrein B. Titres of juvenile hormone I, II and III in Spodoptera littoralis (Noctuidae) from the egg to the pupal moult and their modification by the egg-larval parasitoid Chelonus inanitus (Braconidae). Journal of Insect Physiology. 1999; 45(4): 401–413. pmid:12770366.
  36. 36. Riddiford LM, Hiruma K, Zhou X, Nelson CA. Insights into the molecular basis of the hormonal control of molting and metamorphosis from Manduca sexta and Drosophila melanogaster. Insect Biochemistry & Molecular Biology. 2003; 33(12): 1327–1338. pmid:14599504.
  37. 37. Marchal E, Vandersmissen HP, Badisco L, Van dVS, Verlinden H, Iga M, et al. Control of ecdysteroidogenesis in prothoracic glands of insects: a review. Peptides. 2010; 31(3): 506–519. pmid:19723550.
  38. 38. Ito Y, Yasuda A, Sonobe H. Synthesis and phosphorylation of ecdysteroids during ovarian development in the silkworm, Bombyx mori. Zoolog Sci. 2008; 25(7): 721–77. pmid:18828659.
  39. 39. Merzendorfer H. Insect chitin synthases: a review. Journal of Comparative Physiology B. 2006; 176(1): 1–15. pmid:16075270.
  40. 40. Merzendorfer H, Zimoch L. Chitin metabolism in insects: structure, function and regulation of chitin synthases and chitinases. Journal of Experimental Biology. 2003; 206(Pt 24): 4393–4412. pmid:14610026.
  41. 41. Cohen E. Chitin synthesis and inhibition: a revisit. Pest Management Science. 2001; 57(10): 946–950. pmid:11695188.
  42. 42. Fukamizo T, Kramer KJ. Mechanism of chitin hydrolysis by the binary chitinase system in insect moulting fluid. Insect Biochemistry. 1985; 15(2): 141–145.
  43. 43. Busse MS, Arnold CP, Towb P, Katrivesis J, Wasserman SA. A κB sequence code for pathway-specific innate immune responses. Embo Journal. 2007; 26(16): 3826–3835. pmid:17660749.
  44. 44. Zhong YS, Mita K, Shimada T, Kawasaki H. Glycine-rich protein genes, which encode a major component of the cuticle, have different developmental profiles from other cuticle protein genes in Bombyx mori. Insect Biochemistry & Molecular Biology. 2006; 36(2): 99–110. pmid:16431278.
  45. 45. Okamoto S, Futahashi R, Kojima T, Mita K, Fujiwara H. Catalogue of epidermal genes: Genes expressed in the epidermis during larval molt of the silkworm Bombyx mori. BMC Genomics. 2008; 9(1): 396–410. pmid:18721459.
  46. 46. Nita M, Wang HB, Zhong YS, Mita K, Iwanaga M, Kawasaki H. Analysis of ecdysone-pulse responsive region of BMWCP2 in wing disc of Bombyx mori. Comparative Biochemistry & Physiology Part B Biochemistry & Molecular Biology. 2009; 153(1): 101–108. pmid:19233308.
  47. 47. Charles JP. The regulation of expression of insect cuticle protein genes. Insect Biochemistry & Molecular Biology. 2010; 40(3): 205–213. pmid:20060042.
  48. 48. Ashburner M. Sequential gene activation by ecdysone in polytene chromosomes of Drosophila melanogaster: II. The effects of inhibitors of protein synthesis. Developmental Biology. 1974; 39(39): 141–157. pmid:4209831.
  49. 49. Ashburner M. Sequential gene activation by ecdysone in polytene chromosomes of Drosophila melanogaster: I. Dependence upon ecdysone concentration. Developmental Biology. 1973; 35(1): 47–61. pmid:4207111.
  50. 50. Bonneton F, Chaumot A, Laudet V. Annotation of Tribolium nuclear receptors reveals an increase in evolutionary rate of a network controlling the ecdysone cascade. Insect Biochemistry & Molecular Biology. 2008; 38(4): 416–429. pmid:18342247.
  51. 51. Cheng D, Xia Q, Duan J, Wei L, Huang C, Li Z, et al. Nuclear receptors in Bombyx mori: insights into genomic structure and developmental expression. Insect Biochem Mol Biol. 2008; 38(12): 1130–1137. pmid:18992339.
  52. 52. Cruz J, Nieva C, Manépadrós D, Martín D, Bellés X. Nuclear receptor BgFTZ-F1 regulates molting and the timing of ecdysteroid production during nymphal development in the hemimetabolous insect Blattella germanica. Developmental Dynamics. 2008; 237(11): 3179–3191. pmid:18816439.
  53. 53. Yamada M, Murata T, Hirose S, Lavorgna G, Suzuki E, Ueda H. Temporally restricted expression of transcription factor betaFTZ-F1: significance for embryogenesis, molting and metamorphosis in Drosophila melanogaster. Development. 2000; 127(23): 5083–5092. pmid:11060234.
  54. 54. Kostrouchova M, Krause M, Kostrouch Z, Rall JE. Nuclear hormone receptor CHR3 is a critical regulator of all four larval molts of the nematode Caenorhabditis elegans. Proceedings of the National Academy of Sciences of the United States of America. 2001; 98(13): 7360–7365. pmid:11416209.
  55. 55. Xu J, Tan A, Palli SR. The function of nuclear receptors in regulation of female reproduction and embryogenesis in the red flour beetle, Tribolium castaneum. Journal of Insect Physiology. 2010; 56(56): 1471–1480. pmid:20416316.
  56. 56. Ali MS, Rahman RF, Swapon AH. Transcriptional regulation of cuticular protein glycine-rich13 gene expression in wing disc of Bombyx mori, Lepidoptera. Journal of Insect Science. 2015; 15(1): 71–80. pmid:25843580.
  57. 57. Dhadialla TS, Carlson GR, Le DP. New insecticides with ecdysteroidal and juvenile hormone activity. Entomology. 1998; 43(43): 545–569. pmid:9444757.
  58. 58. Hiruma K, Riddiford LM. Regulation of transcription factors MHR4 and betaFTZ-F1 by 20-hydroxyecdysone during a larval molt in the tobacco hornworm, Manduca sexta. Developmental Biology. 2001; 232(1): 265–274. pmid:11254363.
  59. 59. Enya S, Ameku T, Igarashi F, Iga M, Kataoka H, Shinoda T, et al. A Halloween gene noppera-bo encodes a glutathione S-transferase essential for ecdysteroid biosynthesis via regulating the behaviour of cholesterol in Drosophila. Sci Rep. 2014; 4: 6586–6595. pmid:25300303.
  60. 60. Mesce KA, Fahrbach SE. Integration of endocrine signals that regulate insect ecdysis. Front Neuroendocrinol. 2002; 23(2): 179–199. pmid:11950244.
  61. 61. Zhao XF. Molecular mechanism of insect molting and application. Chinese Bulletin of Entomology. 2007; 44(3): 323–326.
  62. 62. Guo EE, Li S, Cao Y. The Hormonal Regulation Network of Insect Molting. Science of Sericulture. 2008; 34(2): 370–374.