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

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 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 2 nd 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 2 nd 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 |Log 2 (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 1 Premix Ex Taq™ (Tli RNaseH Plus) (TaKaRa, China) and 8 μL of ddH 2 O. After transient centrifugation, qRT-PCR was performed using a LightCycle 96 Realtime 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.

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).
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 6 th 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).

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

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 = |log 2 (nm2/wildtype)| !1 (S3 Table).
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).
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). Differentially expressed genes in nm2 mutant of the silkworm, Bombyx mori

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, Differentially expressed genes in nm2 mutant of the silkworm, Bombyx mori 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.

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], 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 Differentially expressed genes in nm2 mutant of the silkworm, Bombyx mori 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. Differentially expressed genes in nm2 mutant of the silkworm, Bombyx mori Hydrolase is a gene superfamily including protein hydrolase, lipid hydrolase, glycoside hydrolase and nuclear acid hydrolase [26][27][28]. Because nm2 cannot molt in the pre-molting stage of the 2 nd 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 2 nd 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 [29][30][31]. 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).  Presenting one gene, BmCPG10 [15]. Two results from qRT-PCR based on two primers are located at different sites.
**: The results for these genes were contradictory between DGE and qRT-PCR.

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 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 nm2 mutant of the silkworm, Bombyx mori 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 [44][45][46][47]. The ecdysone titer was lower in the nm2 mutant than in the wildtype, and CPGs controlled by ecdysone were therefore expressed differentially.

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 2 nd 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. 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 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 Differentially expressed genes in nm2 mutant of the silkworm, Bombyx mori 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 2 nd 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 1 st instar, another CPG might act as the EPDFP monitor to promote the first molt, explaining why the phenotype of non-molting in the 2 nd 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