Time-Course Transcriptome Analysis Reveals Global Gene Expression Profiling and Dynamic Developmental Signatures across Complete Life Cycle of Bombyx mori

Background: The silkworm (Bombyx mori) is an important lepidopteran model insect worldwide which undergoes a complete metamorphosis developmental process. Although genome sequencing has been long performed, no transcriptome data covering the complete life cycle are available. Methods: Herein, a total of 10 samples were collected consecutively at four different developmental stages, including eggs of 24 h after oviposition (Ed) and eggs of 24 h after artificial egg-hatching (E); larvae from fist to fifth instar (L1–L5); early and late pupa (P4 and P8); and adult moth (M), were subjected to Illumina RNA-Seq and time-course analysis. Results: The summations of the gene expression of the silkworm ten developmental stages show: at Ed stage, eggs develop towards diapause status, the total gene expression level is relatively low; at E stage, after artificial egg-hatching, the expression level improves rapidly; during larval stages from L1–L5, the expression level rises gradually and reaches a peak at L5 stage; during pupae and moth stages, the total gene expression decline stage by stage. The results revealed a dynamical gene expression profile exhibiting significant differential expressions throughout the silkworm life cycle. Moreover, stage-specific key genes were identified at different developmental stages, suggesting their functions mainly characterized in maintaining insect development and immunity homeostasis or driving metamorphosis. GO annotation and KEGG enrichment analysis further revealed the most significantly enriched and fundamentally biological processes during silkworm growth. Conclusion: Collectively, our omics data depicted the first comprehensive landscape of dynamic transcriptome throughout complete developmental processes of B. mori. Our findings also provide valuable references and novel insights into understanding the molecular developmental remodeling events for other Lepidoptera species.


Introduction
The silkworm, Bombyx mori (B. mori), is one of the most economically important insects in the sericulture industry, which has a long history of more than 5000 years of domestication [1]. With the development of life science and biotechnology, silkworm has currently gained attention for its potential applications in medicine and health research [2,3]. For example, the human macular carotenoid transporter (HMCT) was first discovered based on research achievements of B. mori carotenoid transporter system (BCTS) derived from silkworm yellow-cocoon strain [3]. Another case is that silkworm also has advantages as a model organism in toxicological research including health safety and environmental pollution assessment, mainly due to their high sensitivity to chemical compounds such as pesticides, drugs, and heavy metals [4].
Despite these promising practical values mentioned above, silkworms are not suitable for mass breeding due to their diversified genetic characteristics and sensitivity to Processes 2021, 9,1730 2 of 12 environmental condition factors (e.g., temperature and humidity) throughout the complete development life cycle [5]. Similar to most insects such as the mosquito or fruit fly, silkworm spans a complete life cycle characterized by four distinctive developmental stages, including embryo (egg), larva (caterpillar), pupa, and moth (adult) phenotypic metamorphosis within each generation [6][7][8]. Moreover, transitions from different developmental stages of life are usually controlled by tightly reined gene expression and regulation throughout transcriptional, epigenetic, and translational levels [9]. For example, diapause usually occurs in the silkworm egg stage [10]. After oviposition, diapause-destined eggs gradually enter diapause for 10 days, with a series of dramatic changes occurring during the onset of diapause [11]. In general, silkworm larvae are tetramolters that proceed through four instars, molting between each instar [12]. The last instar of silkworm larva is a period of high mulberry consumption and lasts approximately 7-10 days. When silkworms stop eating, they begin to spin cocoons to wrap themselves, which is named as the wandering stage. After 2 or 3 days, the silkworms finish the transformation from larvae to pupae in the cocoons. The newborn pupae have a soft and faint yellow cuticle. Several hours later, the cuticle becomes restrictive and rigid. After about 10 days, the silk moths break the chitinous cover and finish the transition from pupa to moth. However, due to the lack of studies on silkworm developmental stages such as, egg, larva, pupa, and moth causes limitations on understandings of silkworm growth and cultivation in sericulture industrial practice.
Therefore, in the present study, applying RNA-Seq based transcriptome analysis, we investigated gene expression dynamics throughout silkworm (B. mori) complete developmental processes with a coverage of 10 time-points sampling, to provide insights into the whole transcriptome remodeling. The results revealed the importance of manycore genes at specific developmental stages in driving complex metamorphosis in insects. Moreover, our study provides a comprehensive understanding of the physiology and biochemistry of Lepidoptera. It is also of important reference value in gene function research for maintaining homeostasis of development and immunity during insect feeding.

Silkworm Strain and Sample Preparation
The silkworm (B. mori) strain NB was obtained from the School of Life Sciences of Jiangsu University (Zhenjiang, China) and maintained in our laboratory as previously described [13]. Larvae were reared on fresh mulberry under a constant 12:12 day/night photoperiod at 25 ± 1 • C and 75% ± 3% relative humidity. Samples in this study were consecutively collected at four different developmental stages, including: (1) eggs at 24 h after oviposition and eggs at 24 h after artificial egg-hatching (namely Ed and E); (2) larvae from first to fifth instar (namely L1, L2, L3, L4, and L5, respectively); (3) early and late pupae (namely P4 and P8); (4) adult moth (namely M). During the larva stages, samples were collected at the time they had just awoken from dormancy. For each time point, several independent individuals were randomly sampled and homogenized as one integrated pool to minimize the individual genetic differences (Table S1). Samples were put into a 2-mL sterile EP tube after rinsing in RNase-free PBS buffer and dry with filter paper. Then the EP tubes were placed in liquid nitrogen and stored at −80 • C.

RNA Extraction, cDNA Library Construction, and RNA-Sequencing
RNA extraction, cDNA library construction, and RNA-Sequencing were all performed by Majorbio Co., Ltd. (Shanghai, China). Following the qualification of the sample RNAs, mRNAs were enriched using magnetic beads encapsulated Oligo (dT), then fragmentation buffer was used to fragment the mRNAs. After that, six-base random primers were utilized to synthesize cDNA strands, and then two-strand synthesis reactions were performed to construct double-stranded cDNAs. After qualifying the double-stranded cDNAs, the poly-A tail was ligated to the sequences. The connected sequences were then selected according to the Ampure XP Magnetic Beads (Beckman Coulter, Brea, CA, USA) kit instructions. Finally, PCR amplification of 350 bp cDNA fragments is performed to complete the cDNA library assembly. The Illumina Hiseq X Ten sequencer is used to perform paired-end sequencing after the library quality was determined by Agilent 2100. The obtained raw data, after removing low-quality reads, reads with adapters, and reads with the percent of N is greater than 10%, quality control (Q20, Q30, and GC content) was performed. Then the clean reads were mapped to the B. mori reference genome (version: GCA_000151625.1, Feb 2013) (http://metazoa.ensembl.org/Bombyx_mori/Info/Index, release 51 -May 2021) using hisat2 (Version 2.1.0) (http://ccb.jhu.edu/software/hisat2/index.shtml) and TopHat (Version v2.1.1) (http://ccb.jhu.edu/software/tophat/index.shtml). Default parameters were used in Hisat2 and Tophat2. Transcript abundance for each sample was measured by transcripts per kilobase of exon model per million mapped reads (TPM).

Bioinformatics Analysis
Omics data were further analyzed on the free online platform of Majorbio Cloud Platform (www.majorbio.com) [14]. Genes were annotated based on the BLASTX homology searching and comparison against public databases, including NCBI non-redundant (NR), gene ontology (GO), cluster of orthologous groups (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG), Swiss-Prot protein database (Swiss-Prot). The orthologous genes were further annotated to GO terms using the Blast2GO program. KEGG pathway enrichment was performed to identify interesting genes significantly enriched in metabolic or signal transduction pathways using KOBAS software. The significance of enrichment analysis was determined using the hypergeometric tests. An FDR of 0.05 was set as the threshold for screening the significantly enriched GO and KEGG terms.

Experimental Design and Overview of the Transcriptome Sequencing Data
The life cycle of the silkworm is characterized by a series of developmental events. The early larval stage, in which the larvae are fed on raw mulberry leaves, is followed by three distinct pupal stages and finally by the adult moth that emerges from its cocoon after about two weeks. During each of these stages, there are changes in gene expression patterns as well as alterations in cell numbers and morphology. As illustrated in Figure  1, the complete B. mori life expectancy was empirically spanned four discrete stages in their life cycle, including embryo (egg), larva (caterpillar), pupa, and moth (adult), where only the larval stage (from first to fifth larval instars) is a feeding period. To understand molecular mechanisms underlying B. mori developmental process, here we collected wholeanimal samples at a total of 10 representative time points and performed RNA sequencing (RNA-Seq) in this study. The embryonic sampling time points were set based on major stages of embryonic development: embryonic egging stage at 24 h after oviposition (namely Ed), and the late egging stage of embryogenesis at 24 h after artificial egg-hatching (namely E). For the larva stage, from the first to fifth instar larvae (namely L1, L2, L3, L4, and L5) were sampled and examined, respectively. For the pupa stage, pupae were collected separately from Days 4 to 8 after pupation (namely P4 and P8). After adult emergence, one-day-old virgin female moth adults were sampled (namely M). RNAs were extracted from each sample with three biological replicates and further subjected to RNA-sequencing on the Illumina HiSeq 2500 platform. After removing the adaptors and low-quality sequences, a total of 89.85 GB of high-quality clean reads were obtained. The cycle Q30% was >93.4% for each library, with GC contents ranging from 44.03% to 47.16%. After filtering, the clean reads were mapped to the B. mori reference genome (http://metazoa.ensembl.org/Bombyx_mori/Info/Index), with uniquely matching efficiencies ranging from 81.22~89.38%. The expression levels of all transcripts summation at the ten stages were shown by a heatmap (Figure 2A). The red color indicates a higher expression level, while the blue color indicates a lower expression level. The similarity between the ten stages with hierarchical clustering is shown in the above and left of the heatmap. Then we used principal component analysis (PCA) to examine the similarity among the ten groups of data ( Figure 2B). As expected, two egg stages, five larval stages, and three adult stages each clustered together, indicating these stages are closer at gene expression. Principal Component 1 (PC1) which accounted for 47.03% of the variance and PC2 which accounted for 17.06% of the total variance, separated two egg stages, five larval stages, and three adult stages from each other ( Figure 2B). All these quality control results indicated that our transcriptome sequencing data were reliable and acceptable, meeting the basic requirements of the follow-up bioinformatics analysis processes.
reads were mapped to the B. mori reference genome (http://metazoa.ensembl.org/Bombyx_mori/Info/Index), with uniquely matching efficiencies ranging from 81.22~89.38%. The expression levels of all transcripts summation at the ten stages were shown by a heatmap (Figure 2A). The red color indicates a higher expression level, while the blue color indicates a lower expression level. The similarity between the ten stages with hierarchical clustering is shown in the above and left of the heatmap. Then we used principal component analysis (PCA) to examine the similarity among the ten groups of data ( Figure 2B). As expected, two egg stages, five larval stages, and three adult stages each clustered together, indicating these stages are closer at gene expression. Principal Component 1 (PC1) which accounted for 47.03% of the variance and PC2 which accounted for 17.06% of the total variance, separated two egg stages, five larval stages, and three adult stages from each other ( Figure 2B). All these quality control results indicated that our transcriptome sequencing data were reliable and acceptable, meeting the basic requirements of the follow-up bioinformatics analysis processes. Figure 1. Schematic processes indicating a complete life cycle throughout four major developmental stages of B. mori, namely embryo (egg), larva (caterpillar), pupa, and moth (adult) respectively, where only the larva stage (from first to fifth larval instar) is a feeding period. Samples were collected at 10-time points in this study, including eggs at 24 h after oviposition and eggs at 24 h after artificial egg-hatching (Ed and E); larvae from first to fifth instar (L1-L5); early and late pupa (P4 and P8); and adult moth (M), were subjected to RNA-Seq followed by time-course transcriptome analysis. The life span of the B. mori goes through about 37 days in total. Figure 1. Schematic processes indicating a complete life cycle throughout four major developmental stages of B. mori, namely embryo (egg), larva (caterpillar), pupa, and moth (adult) respectively, where only the larva stage (from first to fifth larval instar) is a feeding period. Samples were collected at 10-time points in this study, including eggs at 24 h after oviposition and eggs at 24 h after artificial egg-hatching (Ed and E); larvae from first to fifth instar (L1-L5); early and late pupa (P4 and P8); and adult moth (M), were subjected to RNA-Seq followed by time-course transcriptome analysis. The life span of the B. mori goes through about 37 days in total.

Gene Expression Profile Dynamics throughout Silkworm Complete Life Cycle
Here we first observed global features describing gene expression dynamics during silkworm development life cycle and then compared them with data obtained from other life stages (embryo, larva, pupa, and moth) ( Figure 3). A total 15,745 expressed transcripts were identified in consecutive developmental stages among these samples by comparing with each other using DESeq package software. For further Venn analysis, we intersected four major developmental stages: egg (the union of sets Ed and E), larva (the union of sets L1 to L5), pupa (the union of sets P4 and P8), and moth ( Figure 3A). From the Venn analysis above, 207, 425, 93, and 107 genes were sequentially identified to be involved in four major developmental stages: egg, larva, pupa, and moth, respectively ( Figure 3A). Further GO enrichment analysis for the identified 10236 genes indicated that biological processes including "protein ubiquitination"and "cellular homeostasis" were most significant terms ( Figure 4A). KEGG pathway enrichment analysis showed "ribosome", "protein processing in endoplasmic reticulum" and "ubiquitin-mediated proteolysis" were most significantly enriched pathways ( Figure 4B). In addition, among the 10,236 genes, several genes were found to be constitutively expressed, the expression patterns of these transcripts are highly similar among four developmental stages, suggesting that they are conserved throughout the silkworm life cycle. Throughout the ten developmental stages, the summation of gene expression of all transcripts is increasing from the Ed to the E stage, as the larva stages. From Stages L1 to L5 the integral expression of all genes is upregulated and reaches a peak at the L5 stage, while from Stages P4 to M the expression becomes downregulated ( Figure 3B).

Gene Expression Profile Dynamics throughout Silkworm Complete Life Cycle
Here we first observed global features describing gene expression dynamics durin scripts are highly similar among four developmental stages, suggesting that they are conserved throughout the silkworm life cycle. Throughout the ten developmental stages, the summation of gene expression of all transcripts is increasing from the Ed to the E stage, as the larva stages. From Stages L1 to L5 the integral expression of all genes is upregulated and reaches a peak at the L5 stage, while from Stages P4 to M the expression becomes downregulated ( Figure 3B).   as the larva stages. From Stages L1 to L5 the integral expression of all genes is upregulated and reaches a peak at the L5 stage, while from Stages P4 to M the expression becomes downregulated ( Figure 3B).

Gene Expression Profiles during Embryogenesis
The silkworm egg development is a complex process mainly characterized by diapause, embryonic tropism, and larval maturation. In the present study, eggs from the second day of oviposition and the second day after artificial egg-hatching by mix acid liquid (Ed and E) were collected separately. Gene expression patterns were analyzed using next-generation sequencing technology to identify the expressed genes at Ed and E stages, respectively. The Venn diagrammatical analysis results showed that a total of 1328 specifically expressed genes were identified during the embryonic development process. Among them, 758 genes were found in the Ed Group, 570 genes were only detected in E Group, while 10,174 genes were commonly expressed throughout the silkworm egging stages ( Figure 5A). GO enrichment analysis for these identified commonly expressed genes further indicated that biological processes including "protein ubiquitination", "amino sugar metabolic process" and "glucosamine-containing compound metabolic process" were most involved in the egg development ( Figure 5B). KEGG pathway enrichment analysis showed Processes 2021, 9, 1730 7 of 12 "ribosome", "protein processing in endoplasmic reticulum", "amino sugar and nucleotide sugar metabolism" and "ubiquitin-mediated proteolysis" were significantly enriched pathways ( Figure 5C). In addition, we further annotated key genes related to these biological processes according to their predicted functions. The results showed that most highly expressed genes were associated with developmental processes such as morphogenesis. Some functional genes were found to be regulated by hormones such as ecdysone and juvenile hormone (JH). These results provided important information on specific gene expressions that participate in the silkworm embryonic development process.
( Figure 5A). GO enrichment analysis for these identified commonly expressed genes further indicated that biological processes including "protein ubiquitination", "amino sugar metabolic process" and "glucosamine-containing compound metabolic process" were most involved in the egg development ( Figure 5B). KEGG pathway enrichment analysis showed "ribosome", "protein processing in endoplasmic reticulum", "amino sugar and nucleotide sugar metabolism" and "ubiquitin-mediated proteolysis" were significantly enriched pathways ( Figure 5C). In addition, we further annotated key genes related to these biological processes according to their predicted functions. The results showed that most highly expressed genes were associated with developmental processes such as morphogenesis. Some functional genes were found to be regulated by hormones such as ecdysone and juvenile hormone (JH). These results provided important information on specific gene expressions that participate in the silkworm embryonic development process.

Gene Expression Differences throughout Silkworm Larval Development Stage
Previous studies have shown that the expression of several genes is altered during larval development in silkworms. However, little is known about how these changes occur during different developmental stages or whether they are conserved across instar transition. Here we investigated the transcriptional changes that occur throughout silkworm larval developmental stages (L1-L5) to identify genes whose expression change significantly between different developmental time points. Different genes exhibit different expression patterns at these five stages. A bunch of genes exhibited a wide range of expression patterns, including genes with peak expression at one of the stages from L1 to L5. In addition, there was no consistent pattern in which genes were upregulated or downregulated during each stage. For example, some key transcripts such as those encoding ribosomal proteins were consistently upregulated throughout all five larval instars, suggesting that their levels may be crucial for larval survival or growth rates ( Figure  S1). In addition, gene expression patterns of several key regulators were differentially expressed between L2 and L3 during mulberry leaf feeding. We also found that the number of expressed genes increased significantly during larval instar development, especially

Gene Expression Differences throughout Silkworm Larval Development Stage
Previous studies have shown that the expression of several genes is altered during larval development in silkworms. However, little is known about how these changes occur during different developmental stages or whether they are conserved across instar transition. Here we investigated the transcriptional changes that occur throughout silkworm larval developmental stages (L1-L5) to identify genes whose expression change significantly between different developmental time points. Different genes exhibit different expression patterns at these five stages. A bunch of genes exhibited a wide range of expression patterns, including genes with peak expression at one of the stages from L1 to L5. In addition, there was no consistent pattern in which genes were upregulated or downregulated during each stage. For example, some key transcripts such as those encoding ribosomal proteins were consistently upregulated throughout all five larval instars, suggesting that their levels may be crucial for larval survival or growth rates ( Figure S1). In addition, gene expression patterns of several key regulators were differentially expressed between L2 and L3 during mulberry leaf feeding. We also found that the number of expressed genes increased significantly during larval instar development, especially during the fifth stage (L5) when most silk proteins are synthesized. These candidate genes are most responsible for biological processes such as "cell communication/adhesion/interaction, cellular response to stimulus, macromolecule metabolic process". Functional enrichment analysis indicated gene ontology terms metabolic process (GO:0008152) and catalytic activity (GO:0003824) were more enriched GO terms during larval periods ( Figure 6A). Moreover, multiple developmental genes including Abcdc12a/abcdc12b were identified to be regulated through the Wnt signaling pathway during L5 stages ( Figure 6). Besides, genes such as BGIBMGA003326 and BGIBMGA014298 have a much higher expression at the L5 stage than L1 to L4 stages. Those two genes are annotated as hemolymph juvenile hormone binding protein (JHBP) at the database of clusters of orthologous genes (COGs) and Pfam database, according to previous research, BmJHBPd2 could inhibit the signal transduction of JH [15]. Collectively, these results demonstrate that regulation of these developmen-Processes 2021, 9, 1730 8 of 12 tal genes and biological pathways play important roles in regulating larval growth rate through timing their maturation process to optimize development and productivity.
sides, genes such as BGIBMGA003326 and BGIBMGA014298 have a much higher expression at the L5 stage than L1 to L4 stages. Those two genes are annotated as hemolymph juvenile hormone binding protein (JHBP) at the database of clusters of orthologous genes (COGs) and Pfam database, according to previous research, BmJHBPd2 could inhibit the signal transduction of JH [15]. Collectively, these results demonstrate that regulation of these developmental genes and biological pathways play important roles in regulating larval growth rate through timing their maturation process to optimize development and productivity.

Gene Expression Differences throughout Silkworm Adult Development Stages
Gene expression profiles were characterized throughout silkworm adult development stages from pupal development (P4, P8) to moth (M) stage associated with obvious morphological changes. Our results showed that the 184 genes were only expressed at the P4 stage, 257 genes at the P8 stage, and 700 genes at the M stage, while 9911 genes were commonly expressed throughout the adult development stage ( Figure 7A). Additionally, our GO enrichment analysis identified commonly expressed genes in biological processes including "protein ubiquitination", "chitin metabolic process" and "amino sugar metabolic process", which were involved in the adult development stage ( Figure 7B). KEGG pathway enrichment analysis showed "ribosome", "protein processing in endoplasmic reticulum", "ubiquitin-mediated proteolysis" and "amino sugar and nucleotide sugar metabolism" were the most significantly enriched pathways in the silkworm adult life stages ( Figure 7C). Moreover, the majority of highly expressed genes at pupa stages are those that translate putative cuticle proteins, such as BGIBMGA001862, BGIBMGA005259, and BGIBMGA012861. Among the highly expressed genes at M stage, we also discovered genes that encodes a chemosensory protein, IGF-like peptide, histone H4, and putative defense protein. These genes respectively are BGIBMGA004045, BGIBMGA012305, BGIBMGA009647 and BGIBMGA014360.

Gene Expression Differences throughout Silkworm Adult Development Stages
Gene expression profiles were characterized throughout silkworm adult development stages from pupal development (P4, P8) to moth (M) stage associated with obvious morphological changes. Our results showed that the 184 genes were only expressed at the P4 stage, 257 genes at the P8 stage, and 700 genes at the M stage, while 9911 genes were commonly expressed throughout the adult development stage ( Figure 7A). Additionally, our GO enrichment analysis identified commonly expressed genes in biological processes including "protein ubiquitination", "chitin metabolic process" and "amino sugar metabolic process", which were involved in the adult development stage ( Figure 7B). KEGG pathway enrichment analysis showed "ribosome", "protein processing in endoplasmic reticulum", "ubiquitin-mediated proteolysis" and "amino sugar and nucleotide sugar metabolism" were the most significantly enriched pathways in the silkworm adult life stages ( Figure 7C). Moreover, the majority of highly expressed genes at pupa stages are those that translate putative cuticle proteins, such as BGIBMGA001862, BGIBMGA005259, and BGIBMGA012861. Among the highly expressed genes at M stage, we also discovered genes that encodes a chemosensory protein, IGF-like peptide, histone H4, and putative defense protein. These genes respectively are BGIBMGA004045, BGIBMGA012305, BGIB-MGA009647 and BGIBMGA014360.

Discussion
The silkworm (B. mori) has numerous advantages such as a suitable life cycle with a distinctive boundary of the development stages, moderate body size as well as high reproductive rate, making it a potential model organism for lifespan and developmental research of Lepidoptera insects [16]. Interestingly, the life cycle of silkworm is character-

Discussion
The silkworm (B. mori) has numerous advantages such as a suitable life cycle with a distinctive boundary of the development stages, moderate body size as well as high reproductive rate, making it a potential model organism for lifespan and developmental research of Lepidoptera insects [16]. Interestingly, the life cycle of silkworm is characterized by a series of developmental events, and the transcriptome of each developmental stage has different functions and biological meanings. As illustrated in Figure 1, the complete B. mori life expectancy was empirically spanned four discrete stages in the life cycle, including embryo (egg), larva (caterpillar), pupa, and moth (adult), where only the larval stage (from first to fifth larval instars) is a feeding period. For example, (1) egg is the initial stage from embryogenesis to the larval formation and the quality of silkworm eggs directly affects subsequent development of larvae, pupae, and moths. Diapause occurs in the silkworm egging stage [17]. Genes for basic energy metabolic pathways were expressed more actively during non-feeding, ovum development, and embryogenesis. (2) The larval stage is another next important period for insect growth and development [18]. Therefore, the mulberry feeding efficiency of silkworm larvae is essential as it accounts for nutrients absorption and accumulation [19]. During this process, silkworm larvae show significant weight gain and body-color variations at the first to third larval stages. Typically, silkworms undergo four molts (i.e., tetramolters) and spin cocoons during the late fifth larval instar. The duration of the silkworm larval stage usually depends on B. mori strain and rearing conditions (i.e., temperature and humidity) [12]. (3) The pupa is the transition stage from larva to moth. Similar to many other insects such as mosquito and fruit fly, the most marked biological processes during the silkworm pupa stage were mainly morphogenesis-associated, such as multicellular organism (i.e., eye, antennal sensilla, and wing) development, tube development, and imaginal disc morphogenesis [20]. (4) The moth is the final reproductive stage usually coordinated with three key events in offspring reproduction: to mate, lay eggs, and die.
To fully dissect molecular mechanisms underlying B. mori developmental process, here we presented RNA-Seq based and time-course transcriptome analysis to reveal global gene expression profiling and dynamic developmental signatures across silkworm complete life cycle. We firstly identified the expression patterns of core genes during silkworm complete developmental stages and then investigated how these genes regulating the silkworm life cycle. The results showed that more than 65% of total genes were expressed in all four developmental stages. Functional ontology results revealed that ribosome biogenesis, gene transcription, and translation, mitochondrion, protein binding, metabolic processes, cell cycle were the most enriched terms. For example, most ribosomal genes BGIBMGA010139, BGIBMGA010571, and BGIBMGA002572 were identified to exhibit similar expression patterns during both larval and adult development, suggesting that ribosome biogenesis is an essential factor throughout insect metamorphosis. Apart from ribosomal genes, we also found core genes such as BGIBMGA003186, BGIBMGA007469, which are related to eukaryotic translation initiation function, are also highly and consecutively expressed during all the ten developmental stages. BGIBMGA003186 is eukaryotic translation initiation factor 4A, while BGIBMGA007469 is eukaryotic initiation factor 5A. BGIBMGA004074, whose gene name is BMORCPG19, translates putative cuticle protein.
Silkworm eggs are generally divided into diapause and non-diapause types. The former is characterized by the formation of a thick shell, which has different levels of resistance to desiccation and low sensitivity to light. However, it is not clear how diapaused eggs differ from non-diapaused eggs at molecular levels. To explore the underlying gene expression mechanisms of the silkworm embryonic development, we performed RNA-Seq analysis and compared the transcriptional profiles of eggs from the second day of oviposition and the second day after artificial egg-hatching by mix acid liquid (Ed and E). The results showed that most genes such as BGIBMGA006239, BGIBMGA013413, and BGIBMGA002381 mainly related to metabolism were expressed highly at both Ed and E stages, indicating that metabolic pathways play important roles in silkworm egg devel-opment, regardless of whether they have entered diapause or not. Moreover, the results also showed that there were differences in mRNA levels of some key genes involved in the embryonic process between Ed and E, including BGIBMGA011468, BGIBMGA003444, and BGIBMGA004633. BGIBMGA011468 belongs to the insect hormone biosynthesis pathway. BGIBMGA006146, whose gene name is BmWnt1, expresses higher at the Ed stage than at the E stage. Wnt gene family members associated with the Wnt signaling pathway encoding conserved glycoproteins participate in a variety of biological processes including embryo development, cell proliferation, and differentiation, and tissue regeneration [15]. Hoxa13/Hoxd10 are related to head morphogenesis. BGIBMGA003444 was previously cloned at diapause-destined embryo before [21]. These results indicated that these key genes play important roles in silkworm embryonic development and egging formation.
The silkworm feeds and grows quickly during the larval stages. Therefore, feeding efficiency is very important as it directly affects silkworm larval growth rate and development [19]. In addition, silkworm larvae have become an ideal research model for depicting gene expression dynamics underlying insect major developmental processes such as cell division and morphogenesis. In this study, among the genes continuously expressed throughout the silkworm larval development process, BGIBMGA013945 and BGIBMGA002259 were found to be highly expressed during the first and second larval stages and gradually lower expressed at third, fourth, and fifth larval stages. In silkworms, the third of the fifth instar is considered the boundary for the whole larval development stages [12]. During the fifth instar, the silk gland of the silkworm produces sericin and silk proteins, synthesizing large amounts of silk over several days [12,22]. Our transcriptome data also identified genes in L5 silkworm larvae which were most highly expressed and mainly associated with secretion of proteins, tissue development, and metabolism. These results provided important information about the molecular basis of key events in silkworm larval development of Lepidoptera insect.
After completion of silk spinning, silkworms proceed with larval-pupal metamorphosis [12]. During the pupal stage, insects require large amounts of materials and energy for moth adulating development including sexual maturity and the formation of embryos, as of midgut degeneration. It has been reported that the establishment of energy storage during larval stages are usually allocated for fueling pupal development and early adult stages [23,24]. We further analyzed RNA-Seq results and expression profiles to filter identified genes that were actively transcribed throughout the P4-P8-M developmental stage. For example, BGIBMGA009714 and BGIBMGA012861 were highly expressed in P4 and P8 pupal stages which might translate to putative cuticle proteins. Further, BGIBMGA005781 was expressed in the P4-P8-pupal stage, which translates heat shock protein. The pupal stage is also a critical period of insect growth and development and involves cell migration and division [25]. Our results also indicated that ribosome, protein processing in the endoplasmic reticulum, ubiquitin-mediated proteolysis, amino sugar, and nucleotide sugar metabolism, and calcium signaling pathways were significantly enriched during those periods. Interestingly, during the final stage of silkworm metamorphosis (M), bulk of enzymes, cocoons, for example, one of the serine proteases capable of hydrolyzing silk protein which enables the adult moth to emerge, was largely produced. There are several other similar enzymes such as lysozyme, peptidyl-prolyl cis-trans isomerase, and superoxide dismutase [Cu-Zn].
In conclusion, our omics data have depicted a comprehensive landscape of dynamic gene expression profiling throughout the complete life cycle in B. mori. It is also of important reference value in gene function research for maintaining homeostasis of development and immunity during insect feeding.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/pr9101730/s1, Figure S1: Gene expression of several key genes mentioned in the article, Table  S1: The amount of individuals homogenized in one integrated pool, Table S2: Gene expression of several key genes mentioned in the article.
Author Contributions: X.W. performed the sample collection, library preparation, gene annotation and wrote manuscript. Y.F. and Q.G. and J.X. participated in the bioinformatics and statistical analyses. R.H.T. reviewed this manuscript for language improvement. Y.Y. and K.C. conceived the research, wrote and revised the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This work was supported by the National Natural Science Foundation of China (31861143051, 31900359, and 31872425).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.