Formation Mechanisms of the Ellipsoid Egg in Silkworm (Bombyx mori): Insights from Transcriptomic Profiling

TL;DR

This study explores the genetic basis of elongated egg formation in silkworms using transcriptomic data.

Contribution

The study identifies specific genes and biological processes linked to elongated egg formation in silkworms.

Findings

374 differentially expressed genes were identified in elongated versus normal silkworm eggs.

DEGs were enriched in protein hydrolysis, DNA metabolism, and chromatin organization.

Transcriptional activity was higher in elongated eggs, especially on day 3 of the pupal stage.

Abstract

Background/Objectives: The elongated egg is a morphological mutant of silkworm (Bombyx mori) eggs, yet the biochemical processes and molecular mechanisms underlying this trait remain unclear. Methods: In this study, we performed transcriptome sequencing on the ovaries of female pupae from the Nistari silkworm strain (comparing normal and elongated eggs) during the first three days post-pupation using high-throughput sequencing. Results: A total of 153.56 Gb of filtered data was obtained, identifying 23,366 genes and 35,798 mRNAs. Comparative analysis across three control groups revealed 374 differentially expressed genes (DEGs), with 131 upregulated and 243 downregulated genes in the elongated egg group. Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses indicated that these DEGs were primarily associated with protein hydrolysis, DNA metabolic processes, and euchromatin/heterochromatin organization. Trend expression analysis revealed that transcriptional activity in elongated eggs was significantly higher than in normal eggs, particularly on day 3 of the pupal stage. Conclusions: Weighted gene co-expression network analysis (WGCNA) classified gene expression patterns into twelve modules, with two modules showing specificity. Thirteen hub genes were identified, which are functionally linked to translation initiation, protein density regulation, post-translational modification, and protein turnover. These findings provide foundational insights into the molecular mechanisms driving the formation of the elongated egg in silkworms.

Figures

Full text

1. Introduction

The domestic silkworm (Bombyx mori) is an economically significant lepidopteran insect, playing a vital role in China’s silk economy [1]. The hatching rate is a crucial indicator of silkworm egg quality. It is influenced by several factors, including the season of egg production, nutritional content, and the duration of cold storage [2]. It directly impacts cocoon yield and quality, as well as the economic benefits for farmers [3]. Several factors influence the hatching rate, including genetic, physiological, pathological, and environmental factors [4]. Among these, the quality of the silkworm eggs, such as egg shape and egg weight, significantly impacts the hatching rate [2,5]. 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 [5]. 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 [6]. 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 [7]. This process leads to the coalescence and solidification of the chitin, ultimately forming the eggshell [7]. The normal silkworm egg is typically shortoval in shape, slightly flattened, with a pointed anterior edge at the micropyle end [8]. However, there are various ovoid mutation types in nature, one of which is the elongated egg [9]. The elongated egg is characterized by an extended oval shape, with a polygonal cell structure forming a reticulated pattern on the central surface, and the egg also has a long axis [9,10].

The elongated egg trait exhibits mimetic matrilineal inheritance, with a one-generation delay in phenotypic dominant–recessive segregation compared to normal inheritance [9,10]. This trait is controlled by a recessive mutant gene, elp, located at locus 16.1 on linkage group 18 [11]. 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 [11]. So far, the molecular mechanism underlying the formation of elongated eggs remains unclear. Investigating the molecular processes involved in the formation of elongated silkworm eggs can provide valuable insights into the genetic and molecular regulation of egg shape, offering a theoretical foundation for optimizing silkworm breeding and improving the hatching rate.

In this study, we used normal silkworm eggs from the Nistari breed and its elon-gated 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 they hatched. After hatching, the silkworms were reared on fresh mulberry leaves at 25 °C and 70–85% humidity in darkness. 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. The normal egg group served as the control, and the elongated egg group served as the experimental group.

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 (Invi-trogen, 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, Bremen, 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 DEGs with distinct expression patterns were selected for qRT-PCR verification. Silkworm β-actin 3 (BmActin 3) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) were used as internal reference controls. Primers were designed using Oligo 7.6 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) [12], and clean reads were obtained after filtering. The genome database of the University of Tokyo (https://kaikocdna.dna.affrc.go.jp/) was used as the reference genome sequence. Genome alignment and transcript reconstruction were performed using HISAT2 (2.2.1) [13] and StringTie (2.2.3) [14].

2.4. Differential Expression Gene Screening

Gene expression abundance in each sample was calculated using Cuffdiff [15] and StringTie [14], 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 (2.10) software package [16]. 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

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 KEGG, https://www.genome.jp/kegg/ (accessed on 10 March 2025). 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(GENE DEENOVO, Guangzhou, China) [17] for data analysis.

2.6. Trend Analysis and Gene Co-Expression Network Analysis (WGCNA)

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 (1.3.13) [18]. 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 (2.2.10) [19], and the resulting gene regulatory networks were visualized with Cytoscape (3.9.1) [20]. Module analysis was also conducted to explore the functional relationships between the genes.

3. Results and Analysis

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 (2.2.1) 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) [21]. These results indicate that the sequencing data provide 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. The sequences of the specific primers used are listed in Table 2. 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 three days before pupation are crucial for follicle development [22]. 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. Using the normal egg groups (C1, C2, C3) as a reference, a total of 374 genes were found to be significantly differentially expressed in the elongated egg groups (E1, E2, E3) during the first three days of the pupal stage, including 131 up-regulated genes and 243 down-regulated genes in the elongated egg group. Specifically, in the comparison of C1 vs. E1, there were 60 up-regulated and 35 down-regulated genes in the elongated egg group. In C2 vs. E2, there were 38 up-regulated and 22 down-regulated genes in the elongated egg group. And in C3 vs. E3, there were 49 up-regulated and 192 down-regulated genes in the elongated egg group. A complete list of all differentially expressed genes is provided in Table S1 of the Supplementary Materials. The specific pathways associated with these genes are presented in Table 3. 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

Using the C1 group as a reference, the DEGs of the E1 group were annotated based on the whole transcriptome background, along with 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 all the enriched GO terms, there are 156 up-regulated genes and 70 down-regulated genes. 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).

Using the C2 group as a reference, the DEGs of the E2 group were annotated in the same way, based on the whole transcriptome background along with 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 all the enriched GO terms, there are 72 up-regulated genes and 37 down-regulated genes. 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).

Using the C3 group as a reference, the DEGs of the E3 group were annotated in the same manner. GO enrichment analysis revealed that 311 DEGs were associated with biological processes, 143 with cellular components, and 143 with molecular functions. Among all the enriched GO terms, there are 83 up-regulated genes and 514 down-regulated genes. 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 (Table 2 and Figure 7).

3.5. Trend Analysis

To investigate the expression patterns of DEGs between the control group (C1, C2, C3) and the experimental group (E1, E2, E3), trend analysis was performed using the STEM (1.3.13) [18]. 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 1800 genes, cluster 4 contained 1793 genes, cluster 1 contained 1696 genes, and cluster 10 contained 1002 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 1971 genes, profile 1 containing 1838 genes, profile 10 containing 1646 genes, and profile 5 containing 1333 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) [23] 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.98, p < 0.01), while the pink module showed a strong positive correlation (correlation coefficient = −0.89, p < 0.01) (Figure 9B). Additionally, the scatter plots in Figure 9C,D show the correlations between Module Membership (MM) and Gene Significance (GS) for these two modules, which were 0.6 and 0.91, respectively, indicating strong associations. All the genes contained within the pink module and the yellow modules are included in the Supplementary Materials (Table S2).

To better understand the gene co-expression relationships, we used Cytoscape (3.9.1) software [20] 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 [24]. It is down-regulated only in the C3 vs. E3 comparison. BGIBMGA002542 encodes α-tubulin, a key component of the cytoskeleton [25]. This gene is up-regulated in both the C1 vs. E1 and C3 vs. E3 comparisons. BGIBMGA000462 encodes 27 kDa hemolymph proteins with an unknown function DUF1397 [26,27]. 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 [28]. It is up-regulated only in the C1 vs. E1 comparison. The protein encoded by BGIBMGA000562 (up-regulation) is associated with the sex-linked chocolate recessive color mutation in silkworm larvae [29], but its specific function remains unknown. Chen et al. reported that the expression level of the BGIBMGA004157 (up-regulation) gene is more significant in oocytes during the steroid biosynthesis process, suggesting its involvement in the initiation of diapause [30]. BGIBMGA002058 (down-regulation) is a gene linked to the Z sex chromosome [31]. BGIBMGA011120 (up-regulation) encodes a protein involved in plasmalogen biosynthesis [30], while BGIBMGA001201 (up-regulation) encodes muscular protein 20, which contains a putative actin-binding surface [32]. The functions of the remaining genes, BGIBMGA011121, BGIBMGA012351, BGIBMGA002576, and BGIBMGA006755, have not been reported.

4. Discussion

The morphology of silkworm eggs is determined by the eggshell structure secreted by follicular cells [9,33]. The proliferation, differentiation, and secretory functions of these cells directly influence the final morphology of the eggs [33]. Matsuzaki [34] used electron microscopy to observe that the formation of the silkworm eggshell involves a process where follicular cells gradually secrete substances such as chitin and yolk membrane proteins. The arrangement of follicular cells is closely linked to the elliptical shape of the egg [6,34]. Our GO enrichment analysis revealed that the DEGs during the first 1–3 days of the pupal stage were significantly enriched in processes related to protein hydrolysis, DNA metabolic process, and euchromatin/heterochromatin organization. These biological processes are directly associated with the division and proliferation of follicular cells, as well as the regulation of cell nucleus functions. The activity in the DNA metabolic process may reflect the rapid proliferation of follicular cells during the pupal stage to meet the demand for eggshell synthesis. Additionally, differences in chromatin organization may influence the expression timing of proteins secreted by follicular cells by regulating gene transcription efficiency.

Lu et al. [9] were the first to report that the elongated egg phenotype (elp mutation) in silkworms exhibited maternal mimicry inheritance. The phenotypic segregation ratio was delayed by one generation compared to normal inheritance, and the elp gene was mapped to the 18th linkage group at position 16.1 [9], though its function was not clarified. Our results show that the transcriptional activity of ovarian tissue from elongated eggs was significantly higher than that of normal eggs during the first three days of the pupal stage. Among the DEGs, 243 genes showed downregulated expression, including several candidate genes related to eggshell structural proteins. Wang [35] found that an excessively high ovality index in silkworm eggs leads to a decrease in the fertilization rate, with the fertilization rate of elongated eggs being approximately 30% lower than that of normal eggs. Based on the above information, we hypothesize that the formation of elongated eggs may result from the elp gene inhibiting the expression of specific structural proteins in follicle cells. This prevents the normal short-elliptical shape of the eggshell from forming during deposition, causing it to elongate longitudinally instead. Additionally, the abnormal morphology of the follicle cells, such as longitudinal stretching, may further affect the sperm entry channel (e.g., changes in the structures around the egg pore), ultimately leading to a decrease in the hatching rate.

Furthermore, Telfer [7] proposed that oocyte formation in Lepidoptera insects occurs in three stages: yolk accumulation, eggshell secretion, and follicular cell degeneration, with each stage being regulated by hormones such as ecdysone. KEGG enrichment analysis in this study revealed that DEGs of elongated oocytes were significantly enriched in the MAPK signaling pathway (ko04013) and the Hippo signaling pathway (ko04391) during the first three days of the pupal stage. Both pathways are involved in the regulation of cell proliferation and apoptosis. Previous studies have shown that the MAPK pathway can activate the secretion function of follicular cells via a phosphorylation cascade [30], while the Hippo pathway regulates cell size and tissue morphology by inhibiting YAP/TAZ protein activity [36]. In this study, the differential expression of these pathways (e.g., the key gene BGIBMGA004157 of the MAPK pathway being up-regulated by 2.3 times in elongated oocytes) may disrupt the normal balance of follicular cell proliferation and apoptosis, leading to uncontrolled longitudinal growth and ultimately the formation of elongated oocytes.

Meanwhile, WGCNA identified two specific co-expression modules, revealing that 13 central genes are involved in steroid biosynthesis, phosphatidylcholine synthesis, cytoskeleton composition, and other processes. Notably, the core hub gene BGIBM-GA000562 encodes a protein related to the sex-linked chocolate recessive spot mutation in silkworm larvae [29]. Through searches in other model animal databases, we found that the orthologs of BGIBMGA000562 encode a protein that plays a crucial role in cell signaling and embryonic development [37]. KCTD1 can bind to the important transcription factor AP-2α and recruit corepressors, such as histone deacetylases, to the promoter regions of target genes, thereby inhibiting gene expression activated by AP-2α. This indirectly regulates many genes involved in cell proliferation, differentiation, and apoptosis [38]. Additionally, KCTD1 can promote the degradation of β-catenin in the Wnt signaling pathway through ubiquitination, thus suppressing excessive Wnt signaling, which is essential for normal embryonic development [39]. However, this gene is significantly down-regulated in elongated eggs, particularly on the third day of pupation, which may affect the normal early development of the eggs. Another core gene, BGIBMGA002058, is associated with the Z chromosome [31], which is particularly noteworthy given that Chen et al. [40] reported that the gene responsible for giant egg mutants in the Yun7 silkworm variety is also located on the Z chromosome. This raises the possibility that BGIBMGA002058 contributes to the formation of elongated egg mutants in the Nistari variety, through a mechanism distinct from the elp and sp genes in terms of chromosomal location [41]. Additionally, four core genes (BGIBMGA011121, BGIBMGA012351, BGIBMGA002576, and BGIBMGA006755) remain functionally uncharacterized. Interestingly, these genes, belonging to the pink module, are significantly associated with the GO term eggshell structure formation (e.g., GO:0005576, extracellular matrix). Their expression patterns closely mirror those of known eggshell protein genes, peaking on the third day of the pupal stage. These findings suggest that they may represent novel candidate genes involved in the formation of elongated eggs. Their functions should be further investigated through RNA interference (RNAi) or gene knockout experiments to verify their roles.