Mechanisms Behind the Formation of Ellipsoid Eggs in Silk-Worms (Bombyx mori): Insights from Transcriptomic Profiling

1. Introduction

The domestic silkworm (Bombyx mori) is an economically significant lepidopteran insect, playing a vital role in China's agricultural economy [1]. The hatching rate is a key indicator of silkworm egg quality, influenced by several factors such as egg quality, the season of egg production, nutritional content, and the duration of cold storage [2]. Among these, the shape index and weight of the silkworm eggs are particularly impactful on the hatching rate [2]. Therefore, studying the effect of the shape index on the hatching rate and understanding its molecular mechanisms are of great importance for improving the hatching rate, as well as for the selection and preservation of high-quality silkworm eggs.

The shape of the silkworm egg is determined by the eggshell, which is secreted by the follicular cells [3]. During ovarian development in the female moth, the follicular cells secrete and deposit a thick layer of yolk membrane, followed by the programmed secretion of chitin, which is sequentially deposited outside the oocyte [4]. This process leads to the coalescence and solidification of the chitin, ultimately forming the eggshell [4]. The normal silkworm egg is typically short-oval in shape, slightly flattened, with a pointed anterior edge at the micropyle end [5]. However, there are various oval mutations in nature, one of which is the elongated egg [6]. The elongated egg is characterized by an extended oval shape, with a polygonal cell structure forming a reticulated pattern on the central surface. The long axis of the egg aligns with the egg's long axis [6,7].

The elongated egg trait exhibits mimetic matrilineal inheritance, with a one-generation delay in phenotypic dominant-recessive segregation compared to normal inheritance [6,7]. This trait is controlled by a recessive mutant gene, elp, located at locus 16.1 on linkage group 18 [8]. However, this gene has yet to be finely localized and cloned. Since embryos of the elongated eggs are capable of developing, hatching, and surviving normally, it is hypothesized that the elp gene likely elongates the shape of the follicular cells but does not affect their secretory function [8]. So far, the molecular mechanism underlying the formation of elongated eggs remains unclear.

In this study, we used normal silkworm eggs from the Nistari breed and its elongated egg mutant as research subjects. Ovarian tissues from the first three days after pupation were subjected to high-throughput transcriptome sequencing to analyze differences in transcript levels between normal and elongated eggs during days 1 to 3 of the pupal stage. By identifying differentially expressed genes during this period, we aimed to uncover the key genes associated with the formation of the elongated egg mutation. This research will provide insights into the molecular mechanisms underlying the development of elongated eggs.

2. Materials and Methods

2.1. Animals and Rearing Methods

The normal egg strain of the Nistari variety and the elongated egg strain of silkworms were selected as experimental materials. Both strains were provided by the Key Laboratory of Sericulture, Research Institute of the Chinese Academy of Agricultural Sciences, and Jiangsu University of Science and Technology. The eggs were incubated at 25°C until hatching. After hatching, the silkworms were reared on fresh mulberry leaves at 25°C and 70-85% humidity. At the pupal stage, ovarian tissues were dissected and collected for analysis. Three replicates were taken for each day across a three-day period, with transcriptome sequencing conducted on each of the collected samples. Simultaneously, normal ovarian tissues extracted from the pupal stage during the first three days were labeled as C1, C2, and C3, while elongated ovarian tissues from the same period were labeled as E1, E2, and E3, respectively.

2.2. Extraction of Total RNA and qRT-PCR Analysis

Ovaries were collected from female pupae on days 1 to 3 of pupation, with three biological replicates for each day. Total RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The RNA quality (OD_{260 nm}/OD_{280 nm}) and concentration (ng/μL) were assessed using a Thermo Scientific NanoDropTM 1000 spectrophotometer (Thermo Scientific, Germany). Qualified samples were selected for the construction of sequencing libraries and subsequent Illumina transcriptome high-throughput sequencing.

To validate the reliability of the RNA-Seq data, 20 differentially expressed genes (DEGs) with distinct expression patterns were selected for qRT-PCR verification. Silkworm β-actin 3 and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) were used as internal reference controls. Primers were designed using Oligo 7 software, and each sample included three biological replicates and three technical replicates. The relative expression levels of target genes were calculated using the 2^−ΔΔCT method, based on the Cq values of the target and reference genes.

2.3. Transcriptome Sequencing and Assembly

Transcriptome sequencing of normal and elongated eggs from ovaries during the first three days of the pupal stage was conducted by Nanjing Jisi Huiyuan Biotechnology Co., Ltd. in Nanjing, China. Sequencing was performed with three biological replicates, and 18 libraries were constructed. After sequencing, the raw reads were subjected to quality control (QC) using the NGS QC Toolkit (v2.3.3) [9], and clean reads were obtained after filtering. The genome database of the University of Tokyo (http://kaikocdna.dna.affrc.go.jp/) was used as the reference genome sequence. Genome alignment and transcript reconstruction were performed using HISAT2 [10] and StringTie [11].

2.4. Differential Expression Gene Screening

Gene expression abundance in each sample was calculated using Cuffdiff [12] and StringTie [11], based on sequence similarity comparison. The expression levels of transcripts or genes were measured using the FPKM method. Differential expression analysis was performed using the Negative Binomial Distribution Test (NB) in the DESeq software package [13]. The NB test assessed the significance of differences in read counts, with the base mean value used to estimate gene expression. Differentially expressed genes were selected using the following criteria: |log2(Fold Change)| ≥1 and FDR < 0.05. To correct for multiple testing, the Benjamini-Hochberg method was applied to adjust the p-values, and FDR was used as the key index for identifying differentially expressed genes.

2.5. GO and KEGG Functional Enrichment

Gene Ontology (GO, http://www.geneontology.org/) is an international standard classification system for gene function, divided into three categories: molecular function, biological process, and cellular component. After screening for differentially expressed genes, GO enrichment analysis was conducted to investigate their distribution across these categories and to clarify the functional differences in the samples. The enrichment analysis is based on the hypergeometric distribution, which calculates the probability of observing a given number of differential genes within a specific GO term compared to the whole genome. A small p-value indicates significant enrichment of differential genes in that GO category.

KEGG pathway analysis was conducted on the differentially expressed genes using annotations from the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/). The significance of gene enrichment in each pathway was assessed using the hypergeometric distribution test, where a small p-value indicates significant enrichment in that pathway. Both GO and KEGG enrichment analyses were performed and visualized using the OmicShare online platform [14] for data analysis.

2.6. Trend Analysis and Gene Co-Expression Network Analysis

Expression trend analysis of differentially expressed genes during the first three days of the pupal stage (days 1, 2, and 3) was performed using the STEM software [15]. The analysis used default parameters, with FDR < 5% and a fold change greater than 2 between any two time points. Gene co-expression network analysis was carried out using the WGCNA module in TBtools software [16], and the resulting gene regulatory networks were visualized with Cytoscape [17]. Module analysis was also conducted to explore the functional relationships between the genes.

3. Results

3.1. Quality Control and Sequence Splicing of Transcriptome Data

Using sequencing-while-synthesizing technology, each sample generated approximately 32-35 million raw reads, with a total read length of about 3 Gb, approximately seven times the size of the silkworm genome. Most of the bases achieved or exceeded the Q30 quality score. After removing low-quality reads (where more than 50% of the bases had a quality score of Q ≤ 5), a total of 153.56 Gb of clean data was obtained. The base quality evaluation showed that 89.77% of the bases had a Q30 score, as detailed in Table 1.

Reference genome alignment using HISAT2 revealed that about 81%-84% of the clean reads from each sample aligned with the silkworm genome. The coverage of all detected genes was approximately 80%, with 45%-60% of the reads mapping to known exonic regions, 11%-31% to intronic regions, and 7%-12% to intergenic regions. After merging our sequencing data with publicly available mRNA sequences, genome splicing and gene reconstruction were performed using StringTie. This analysis predicted a total of 23,366 genes and 35,798 mRNAs, significantly surpassing the 14,624 genes predicted by the Silkworm Genome Database (SilkDB) [18]. These results indicate that the sequencing data provides a more comprehensive representation of gene expression in silkworm tissues.

3.2. qPCR Assay and DEGs Screening

The expression levels of 20 genes were assessed using qRT-PCR according to the manufacturer's instructions. A strong correlation (r = 0.939, Figure 1) was observed between the RNA-Seq data and the qRT-PCR results. The expression profiles obtained from qRT-PCR closely matched those derived from RNA-Seq, validating the reliability of the RNA-Seq data in this study.

The elongated eggs of the silkworm are a characteristic of the eggshell trait, secreted by the egg follicle cells. The three days before pupation are crucial for follicle development [19]. To investigate the gene expression patterns associated with this trait, a comparative analysis of normal and elongated egg samples was performed. The data were normalized, and differentially expressed genes (DEGs) were identified using FDR ≤ 0.05 and log2 fold change > 1 as screening criteria. A total of 374 DEGs were detected between the two groups during the first three days of the pupal stage, including 131 up-regulated and 243 down-regulated genes. Specifically, the C1 vs E1 group showed 60 up-regulated and 35 down-regulated genes. The C2 vs E2 group showed 38 up-regulated and 22 down-regulated genes, while the C3 vs E3 group had 49 up-regulated and 192 down-regulated genes. The specific pathways associated with these genes are presented in Table 2. Differential expression was further analyzed using Venn and MA plots (Figure 2 and Figure 3), which revealed changes in gene expression patterns between the two groups during the early pupal stage.

3.3. GO Enrichment Analysis of DEGs

In the C1 vs E1 group, DEGs were annotated using the whole transcriptome background and GO/KO annotations. GO enrichment analysis revealed that 138 DEGs were associated with biological processes, 52 with molecular functions, and 36 with cellular components. Among the biological process-related genes, most were involved in multicellular biological processes, protein hydrolysis, DNA metabolism, and GPI anchor biosynthesis. Cellular component-related genes were mainly involved in signal transduction, heterochromatin, euchromatin, integral components of the membrane, and multilineage chromosome bands. For molecular functions, the majority of genes were related to RNA-directed DNA polymerase activity, peptidyl chain endonuclease activity, DNA binding, and hydrolase activity (Figure 4).

In the C2 vs E2 group, DEGs were annotated similarly using the whole transcriptome background and GO/KO annotations. GO enrichment analysis showed that 63 DEGs were assigned to biological processes, 18 to cellular components, and 45 to molecular functions. Among the biological process-related genes, most were involved in RNA-directed DNA polymerase activity, somatic muscle development, protein hydrolysis, and DNA integration. The cellular component-related genes were primarily involved in signal transduction, membrane integration, cytoplasm, euchromatin, and protein complex. Genes related to molecular functions were mainly associated with RNA-directed DNA polymerase activity, zinc-finger structure, DNA binding, and neuropeptide Y receptor activity (Figure 5).

For the C3 vs E3 group, DEGs were annotated in the same way, and GO enrichment analysis revealed that 311 DEGs were linked to biological processes, 143 to cellular components, and 143 to molecular functions. Biological process-related genes were predominantly involved in translation initiation, RNA phosphodiester bond hydrolysis, somatic muscle development, and nucleus development. Cellular component-related genes mainly contributed to signal transduction, heterochromatin, negative regulation of smooth signaling pathways, dorsal closure, and ion transport across membranes. Molecular function-related genes were primarily linked to monooxygenase activity, deconjugating enzyme activity, cysteine-type peptidyl endonuclease activity, and calcium-dependent phospholipid binding (Figure 6).

3.4. KEGG Pathway Enrichment Analysis of DEGs

KEGG enrichment analysis identified a total of 123 DEGs annotated to 53 KEGG pathways across all three sample pairs. After filtering for pathways with a p-value ≤ 0.05, several pathways were found to be significantly enriched in DEGs. In the C1 vs E1 group, the pathways mismatch repair, fanconi anemia pathway, metabolic pathway, and GPI-anchor biosynthesis were significantly enriched. The C2 vs E2 group enriched only one KEGG pathway, namely, endocytosis. In the C3 vs E3 group, seven pathways were significantly enriched, including those related to cysteine and methionine metabolism, protein processing in the endoplasmic reticulum, spliceosome, endocytosis, ribosome biosynthesis, and GPI-anchor biosynthesis (Figure 7). These enriched pathways provide valuable insights into the specific processes, functions, and pathways associated with stunting in the silkworm B. mori.

3.5. Trend Analysis

To investigate the expression patterns of differentially expressed genes (DEGs) in the normal egg group (C1, C2, C3) and the elongated egg group (E1, E2, E3), trend analysis was performed using the Short Time-series Expression Miner (STEM) software [20]. The 13,678 DEGs in the normal egg group were clustered into 10 expression patterns, with clusters 8, 4, 1, and 10 showing significantly enriched trends (p < 0.05). Specifically, cluster 8 contained 1,800 genes, cluster 4 contained 1,793 genes, cluster 1 contained 1,696 genes, and cluster 10 contained 1,002 genes (Figure 8). In the elongated egg group, 14,005 DEGs were clustered into 10 expression patterns. Profiles 4, 1, 10, and 5 exhibited significant enrichment trends (p < 0.05), with profile 4 containing 1,971 genes, profile 1 containing 1,838 genes, profile 10 containing 1,646 genes, and profile 5 containing 1,333 genes (Figure 8).

For each expression pattern, GO functional enrichment analysis was performed. In the normal egg group, genes in profiles 4 and 10 contained the majority of GO function entries, while genes in profile 1 were mainly enriched in cellular components and molecular functions. In the elongated egg group, profiles 4 and 1 contained the most GO function entries, while genes in profile 5 were enriched in cellular components and molecular functions. Additionally, genes in profile 8 were concentrated in the biological process and molecular function categories. In terms of expression trends, the first three days after pupation showed higher transcriptional activity in the elongated eggs compared to the normal eggs, particularly on the third day of pupation.

3.6. Weighted Gene Co-Expression Network Analysis (WGCNA)

To further investigate the relationship between gene expression in the ovaries during the first three days of the pupal stage and the formation of elongated eggs, we employed Weighted Gene Co-expression Network Analysis (WGCNA) [21] to examine the expression patterns of 23,366 genes filtered from 18 samples. This analysis identified 12 distinct gene co-expression modules based on their expression patterns (Figure 9A). Among them, two modules—yellow and pink—exhibited specific gene expression characteristics. The yellow module showed a negative correlation with the number of days until pupation (correlation coefficient = −0.89, p < 0.01), while the pink module showed a strong positive correlation (correlation coefficient = 0.59, p < 0.01). Additionally, the correlations between Module Membership (MM) and Gene Significance (GS) for these two modules were 0.6 and 0.91, respectively, indicating robust associations (Figure 9B and 9C).

To better understand the gene co-expression relationships, we used Cytoscape software [17] to visualize the co-expression network of differentially expressed genes with a weight greater than 0.2 in the yellow and pink modules. In the yellow module, four hub genes were identified: BGIBMGA005781, BGIBMGA002542, BGIBMGA000462, and BGIBMGA013932 (Figure 10). In contrast, the pink module contained nine core genes: BGIBMGA000562, BGIBMGA004157, BGIBMGA002058, BGIBMGA011121, BGIBMGA011120, BGIBMGA012351, BGIBMGA002576, BGIBMGA001201, and BGIBMGA006755 (Figure 11).

The core gene BGIBMGA005781 encodes a heat shock protein-like protein, which may play an important role in the formation of lipid bodies in silkworms [22]. BGIBMGA002542 encodes α-tubulin, a key component of the cytoskeleton [23]. BGIBMGA000462 encodes 27-kDa hemolymph proteins with an unknown function DUF1397 [24,25]. However, the specific function of this protein has not yet been reported. BGIBMGA013932 encodes a Serpin7 protein, which is involved in inhibiting the degradation of silk proteins and the generation of oxidative radicals, thereby maintaining homeostasis during silk formation [26]. The protein encoded by BGIBMGA000562 is associated with the sex-linked chocolate recessive color mutation in silkworm larvae [27], but its specific function remains unknown. Chen et al. reported that the expression level of the BGIBMGA004157 gene is more significant in oocytes during the steroid biosynthesis process, suggesting its involvement in the initiation of diapause [28]. BGIBMGA002058 is a gene linked to the Z sex chromosome [29]. BGIBMGA011120 encodes a protein involved in plasmalogen biosynthesis [28], while BGIBMGA001201 encodes muscular protein 20, which contains a putative actin-binding surface [30]. The functions of the remaining genes, BGIBMGA011121, BGIBMGA012351, BGIBMGA002576, and BGIBMGA006755, have not been reported.

4. Discussion

The hatching rate of silkworm eggs is a crucial indicator for evaluating the quality of silkworm breeding stock, directly affecting the yield, quality of cocoons, and the economic benefits for farmers [31]. Several factors influence the hatching rate, including genetic, physiological, pathological, and environmental factors [32]. Among these, the quality of the silkworm eggs, such as egg shape and egg weight, significantly impacts the hatching rate [33]. An excessively large egg shape index or a relatively long egg shape can hinder fertilization, leading to a decrease in the fertilization rate. Conversely, an excessively small egg shape index may cause difficulties in embryo development, leading to mortality in the later stages, thus reducing the hatching rate [33]. Therefore, investigating the molecular mechanisms underlying the formation of elongated silkworm eggs can provide valuable insights into the genetic and molecular basis of egg shape regulation, offering a theoretical foundation for optimizing silkworm breeding and enhancing the hatching rate.

In previous studies, we identified a mutant silkworm from the Nistari variety, which exhibited oval-shaped eggs that elongated into a more cylindrical form [33]. The ratio of the long axis to the short axis of these mutant eggs was 1.7, significantly higher than the 1.26 ratio of normal eggs, and the hatching rate was 59.49%, much lower than the 86.43% observed in normal eggs [34]. However, the molecular mechanism underlying the formation of elongated eggs remains unclear. In this study, we used normal Nistari eggs as the control group and performed transcriptome sequencing on the ovaries of both elongated and normal eggs from the first to the third day of the pupal stage. We then identified and analyzed differentially expressed genes (DEGs) between the two groups through GO and KEGG enrichment analysis. The results indicated that the DEGs were primarily involved in DNA replication, protein synthesis and degradation, mRNA splicing, snRNA regulatory complex activity, and the composition of the cytoskeleton and chromosomes. From a cellular perspective, the formation of elongated eggs involves the stretching of normal follicles into an elongated shape, a process controlled by the elp gene [8]. This process likely requires cytoskeletal rearrangement and energy conversion, which may, in turn, influence the transcriptional regulation of related genes.

Through WGCNA analysis, two specific co-expression modules were identified, revealing 13 hub genes that play roles in processes such as steroid biosynthesis, phosphatidylcholine synthesis, and cytoskeleton composition. Notably, the core hub gene BGIBMGA000562 encodes a protein related to the sex-linked chocolate recessive spot mutation in silkworm larvae [27], although its specific function remains to be further investigated. The core gene BGIBMGA002058 is associated with the Z chromosome [29], which is of particular interest, as Chen et al. reported that the target gene responsible for the formation of giant egg mutants in the Yun7 silkworm variety is also located on the Z chromosome [35]. This suggests that BGIBMGA002058 may play a role in the formation of elongated egg mutants in the Nistari silkworm variety, a mechanism that differs significantly from the elp and sp genes in terms of chromosomal location [36]. Additionally, the functions of the remaining four core genes (BGIBMGA011121, BGIBMGA012351, BGIBMGA002576, and BGIBMGA006755) have yet to be reported. Further research into the functions of these genes may help uncover the molecular mechanisms behind the formation of elongated eggs.