Selection and evaluation of dual reference genes in silkworm (Bombyx mori ) based on transcriptome data

Abstract

Objective

The stability of reference genes is critical for qPCR data accuracy, especially in insect developmental studies with stage-specific transcriptional variation. This study aimed to identify stable reference genes for Bombyx mori’s entire developmental cycle and validate their utility in target gene normalization.

Methods

Transcriptomic data covering the entire developmental cycle of B. mori were screened for stable reference genes using principal component analysis (PCA) and Venn analysis. A total of 125 highly expressed and stable candidate genes were identified and nine candidate reference genes were selected based on the coefficient of variation (CV), with actin included as a traditional control. The stability of these nine genes was systematically evaluated using three algorithms: Delta CT, BestKeeper, and NormFinder, while the geNorm algorithm determined the optimal reference gene pair. The performance of the identified dual-reference system was validated by normalizing seven target genes and comparing results with normalization using actin alone.

Results

The optimal reference genes varied across the developmental stages of B. mori. Although no single gene was universally stable during all stages, the geNorm algorithm identified two genes (006219 and 000526) as the most stable pair for normalization throughout the entire developmental cycle. Normalization of seven target genes using this dual-reference system yielded consistent and reliable expression trends, whereas normalization with actin alone resulted in significant expression biases across developmental stages.

Conclusions

The dual-reference gene system effectively minimizes technical errors, enhancing the accuracy and reproducibility of qPCR data. Actin is not a universally stable reference gene for B. mori, and this study provides a reliable set of reference genes for qPCRbased gene expression analysis in B. mori and that serves as a valuable reference for similar studies on other insects. Further validation of this dual-reference system under diverse experimental conditions will contribute to advanced molecular research in insects.

1. Introduction

Quantitative PCR (qPCR) is a crucial tool for measuring gene expression levels in biological systems, with its accuracy depending on high-quality RNA and reference gene stabilitys. ¹ In holometabolous insects, development progresses through four distinct stages, namely egg, larva, pupa, and adult, and is accompanied by organ remodelling that results in stage-specific gene expression patterns,^2,3where, genes essential in lifecycle regulation are stably expressed across all these stages.^4,5 The silkworm (Bombyx mori), is a model holometabolous insect and widely used in studies on insect growth, development, and metabolic regulation.^6,7

Transcriptomic technology (RNA-Seq) enables the analysis of gene expression across different developmental stages and tissues in insects, ⁸ and stably expressed genes that serve as reference genes for RT-qPCR have been studied in various species.^9–11 Current transcriptomic studies on B. mori mainly focus either on specific organs or selected developmental stages,^12,13 while the identification of stable reference genes across the entire lifecycle remains unresolved. Consequently, researchers often rely on actin or Tubulin as reference genes, but actin expression fluctuates due to organ remodelling, potentially compromising the accuracy and reliability of qPCR results. ¹⁴ Therefore, identifying reference genes that maintain stable expression throughout all the developmental stages of B. mori, independent of metamorphic transitions-is essential for reliable data normalization.

Previous studies using oligonucleotide microarray technology have identified transcription initiation factors as potential qPCR reference genes in B. mori.^15,16 The stability of reference genes can be assessed using the RefFinder online tool, which integrates the Delta CT, NormFinder, BestKeeper, and geNorm algorithms-to assign weighted scores and calculate geometric means for ranking candidate genes. ¹⁷ This approach has been successfully applied for the selection of reference genes in Abelmoschus esculentus (okra), ¹⁸ and Pseudomonas aeruginosa. ¹⁹

In this study, integrated transcriptomic data from the eggs, all larval instars, pupae, and adults of B. mori were used to identify candidate reference genes with fragments per kilobase per million (FPKM) ≥ 190 and low coefficient of variation (CV < 0.3). The stability of these candidates was further evaluated with the RefFinder tool to establish a dual-reference gene normalization system suitable for all the developmental stages of B. mori. The universal applicability of the 006219 and 000526 gene combination was validated in cross-stage experiments to provide a methodological framework for the precision analysis of developmental regulatory networks in B. mori. This serves as a reference for the selection of stable reference genes in other holometabolous insects.

2. Materials and methods

2.1. Diapause release and silkworm rearing

Diapause eggs of the NB strain Bombyx mori were preserved at 4 °C for three months and maintained at the School of Life Sciences, Jiangsu University. A total of 30 silkworm eggs were treated per batch for diapause termination, with three independent biological replicates set up in this experiment. The diapause release followed emersion in a formaldehyde solution for five minutes at room temperature, transfer to a hydrochloric acid solution (specific gravity 1.075) containing 2% formaldehyde and incubation at 46°C for 15 minutes and rinsing under slowrunning water for 40 minutes at room temperature. Hatched larvae were reared at 25°C and 70% humidity, with a 12-hour light/12-hour dark photoperiod, and fed fresh mulberry leaves until pupation. Mulberry leaves were sampled from the mulberry plantation managed by our research group at Jiangsu University. With long-term cultivation, the mulberry trees maintain stable leaf quality, and only mature mulberry leaves were provided for larval feeding of Bombyx mori. A total of 30 Bombyx mori larvae were maintained under group rearing conditions in an 80 cm × 50 cm breeding container.

2.2. Sample collection and RNA sequencing

Samples were collected from each developmental stage as follows: diapause eggs were sampled at 24 hours post-oviposition; activated eggs were collected at 24 hours post-diapause release; and larvae at each instar (1st to 5th instar) were sampled at 24 hours after moulting. Additional samples were taken on the fourth day (Pupa D4) and eighth day (Pupa D8) after pupation, as well as 24 hours post-eclosion (Moth). The pupation process of Bombyx mori was monitored and recorded by visual observation. Each sample consisted of 30 eggs pooled for each of the Diapause Egg and Activated Egg stages, 30 larvae for Instar 1, and five individuals for Instar 2, Instar 3, Instar 4, Instar 5, Pupa D4, and Pupa D8. Three biological replicates were set for each developmental stage. All samples were pooled together, and only one set of representative data was used for subsequent analysis. Each sample was supplemented with 1 mL of Trizol reagent, placed in a 1.5 mL EP tube, and stored at −80 °C in a refrigerator.

Total RNA was extracted from pooled samples at each time point using the phenolchloroform method, and sequencing libraries were constructed using 1 μg of total RNA and sequenced on the Illumina HiSeq X Ten platform of Majorbio (Shanghai, China). Sequencing data were analysed using the Majorbio Cloud platform, with reference to the B. mori genome (https://metazoa.ensembl.org/Bombyx_mori/Info/Index, Version: GCA_000151625.1). Pairedend sequencing was performed on the Illumina NovaSeq 6000 platform with an insert size of approximately 350 bp. Each sample obtained more than 6 Gb of clean data, which provided sufficient sequencing depth to support gene expression quantification and transcript structure analysis of Bombyx mori. Raw data were filtered following standard procedures to generate clean data. The bioinformatics tools used included HISAT2 (v2.1.0) and TopHat (v2.1.1) for genome alignment, RSeQC-2.3.6 for sequencing quality evaluation, StringTie for transcript assembly, and RSEM for gene expression quantification.

2.3. RNA extraction and RT-qPCR

Validation of candidate reference gene stability used, samples collected from diapause eggs (Egg), each larval instars (sampled 24 hours after moulting), pupae on the fourth day after pupation (Pupa), and moths one day after eclosion (Moth). The RNA was extracted from each sample using the phenol-chloroform method, and first-strand cDNA was synthesized from 1 μg of total RNA using the HiScript II 1st Strand cDNA Synthesis Kit (+gDNA wiper) (Vazyme, China). A qPCR was performed using the 2×SYBR Green qPCR Mix, and the cycle threshold (CT) values of the candidate reference genes were determined across different developmental stages. Primer sequences for the candidate reference genes are listed in Table S1. The qPCR cycling conditions were performed as follows: initial denaturation at 94 °C for 3 min (1 cycle), followed by 40 amplification cycles consisting of denaturation at 94 °C for 10 s and annealing/extension at 60 °C for 40 s, and ending with a melting curve analysis program including 95 °C for 15 s (1 cycle), 60 °C for 1 min (1 cycle), and 95 °C for 1 s (1 cycle). Three technical replicates were performed for each qPCR sample to reduce experimental and instrumental variation and ensure the reliability of quantitative results. The relative expression levels of the target genes were determined using the 2^{-^ΔΔCT} comparative threshold cycle method. All experiments were conducted at 4°C.

2.4. CT stability evaluation

The CT values obtained from the qPCR experiments were uploaded to RefFinder (https://blooge.cn/RefFinder/) to compute comprehensive rankings and evaluate the expression stability of candidate reference genes across the examined developmental stages.

2.5. Statistical analyses

A one-way ANOVA method was used to analyse the data, and statistical significance is indicated by *, **, and *** for p < 0.05, 0.01, and 0.001, respectively, while ‘ns’ denotes no significant difference. Data are presented as mean ± standard deviation (SD) (n = 3). Graphs and charts were created in Graphpad Prism (v10.2.2). Gene IDs were searched and their corresponding gene names were identified using the SilkDB 3.0 database (https://silkdb.bioinfotoolkits.net/main).

3. Results

3.1. Screening of candidate reference genes based on transcriptome data

Principal component analysis (PCA) revealed that diapause eggs clustered separately from activated eggs, while the larval stages (Instar 1 to Instar 5) formed another distinct group, and pupae and moths were grouped. PC1, PC2, and PC3 explained 42.33%, 16.00%, and 11.33% of the total variance, respectively, together accounting for approximately 69.66% of the transcriptional variation. This suggests clear distinctions between embryonic, larval, and adult life stages while maintaining high similarity within each group (Figure 1(a)).

Figure 1.

Transcriptome Sequencing Analysis of Candidate Reference Genes in Bombyx mori. (a) PCA illustrating the transcriptional separation and sample clustering of distinct developmental stages (Instar 1-5, Pupa D4/D8, Moth, Activated egg, and Diapause egg). (b) Venn diagram showing the overlap of genes with stable expression (FPKM ≥ 190) across all developmental stages, identifying 125 core candidate genes. (c) GO enrichment analysis of the 125 candidate genes, depicting significantly enriched functional categories. (d) KEGG pathway enrichment analysis highlighting the top biological pathways associated with the candidate genes. (e) GO functional annotation classification of the 125 candidate genes, grouped by biological process, cellular component, and molecular function. (f) KEGG pathway annotation classification of the 125 candidate genes, summarizing their involvement in major metabolic and cellular processes.

To identify genes that are stably expressed throughout all the developmental stages with relatively high expression levels, an expression threshold of FPKM ≥ 190 was applied. Using FPKM ≥ 190 as the high expression threshold, a total of approximately 854 highly expressed genes were identified. Among them, 125 genes were stably expressed across all developmental stages, accounting for 14.64% of all highly expressed genes. Venn analysis identified 125 genes meeting this criterion (Figure 1(b)), and Gene Ontology (GO) enrichment analysis indicated that these genes were primarily associated with translation, peptide biosynthesis process, and structural constituent of ribosome. Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis further revealed that these genes were enriched in pathways related to Ribosome, Oxidative phosphorylation, and Parkinson’s disease (Figure 1(c) and (d)).

The GO functional annotation also revealed that in the molecular function category, the genes were associated with structural molecule activity, binding, and catalytic activity; and in the cellular component category, they were linked to cell part, protein-containing complex, and organelle. In the biological process category, they were primarily involved in metabolic processes, cellular processes, and localization (Figure 1(e)). The KEGG pathway annotation revealed that these genes were primarily involved in translation under the Genetic Information Processing module, energy metabolism under the Metabolism module, transport and catabolism under the Cellular Processes module, environmental adaptation under the Organismal Systems module, and neurodegenerative diseases under the Human Diseases module (Figure 1(f)).

3.2. Optimal candidate reference genes selection

To ensure stable expression of reference genes, the following selection criteria were established: (i) The gene’s expression level must exceed 190 in all samples; (ii) The gene’s expression level must show no significant variation across samples; (iii) The coefficient of variation (CV) of gene expression must be low, indicating higher stability.

Genes identified through sequencing were originally named in the format BGIBMGAXXXXXX, which was simplified to XXXXXX for further analysis. Results showed that 125 genes were stably expressed and met the criteria listed in the Methods section, where the top nine genes with the lowest CV values were selected as candidate reference genes. Actin (005576) was also included for comparison, resulting in a total of 10 candidate reference genes. The CV values of these genes are as follows: 002429 (cyclophilinA) is 0.171, 000526 is 0.184, 006219 is 0.192, 007469 (eukaryotic initiation factor 5A) is 0.203, 002493 is 0.213, 010751 is. 0.214, 010604 (Ribosomal protein L40) is 0.232, 006794 (tmbim6) is 0.244, 007043 is 0.248463, 005576 (actin) is 0.896. Among them, the gene with the lowest expression level was 007043, while 010604 (Ribosomal protein L40) exhibited the highest expression. The most stable gene was 002429 (cyclophilinA) (CV = 0.171), whereas actin showed the highest variability (CV = 0.896) (Figure 2).

Figure 2.

Candidate Reference Gene Analysis in Bombyx mori. (a) Cumulative expression levels of nine candidate reference genes and the actin gene across different developmental stages. The line plot represents the coefficient of variation (CV) for each candidate reference gene at different stages. (b) GO enrichment analysis of the nine candidate reference genes and actin. (c) KEGG enrichment analysis of the nine candidate reference genes and actin.

The GO enrichment revealed that these genes were primarily involved in translational frameshifting, regulation of translational termination, positive regulation of translational termination, and positive regulation of translational elongation. The KEGG analysis indicated that these genes were associated with Toxoplasmosis and the cCMP-PKC signalling pathway. Additionally, seven target genes were randomly selected from the complete genome to evaluate the impact of these 10 candidate reference genes on the normalization of gene expression (Table 1).

Table 1.

Basic information of the seven candidate genes.

Gene_id	Gene name	Gene_description
BGIBMGA001173	- -	--
BGIBMGA002350	aurora kinase B-like	--
BGIBMGA005098	--	--
BGIBMGA005841	mitochondrial ribosomal protein S21	--
BGIBMGA006280	cecropin A	Cecropin-A [Source:UniProtKB/Swiss-Prot; Acc: Q27239]
BGIBMGA007363	ribosomal protein S26	Ribosomal protein S26 [Source:UniProtKB/TrEMBL; Acc: Q5UAL7]
BGIBMGA013953	--	--

3.3. Ranking and selection of optimal reference genes

To assess the stability of the 10 selected candidate reference genes across different developmental stages of Bombyx mori, samples were collected from eggs, first to fifth instar larvae (1-day post-moulting), pupae (4 days post-pupation), and adults (1-day post-eclosion). The expression levels of the candidate reference genes were measured using RT-qPCR, and stability was evaluated using the delta CT, BestKeeper, and NormFinder algorithms, resulting in a comprehensive ranking.

The CT values of the 10 candidate reference genes ranged from 14 to 30. The lowest CT value was observed for BGIBMGA010604 (Ribosomal protein L40) at the moth stage (14.603 ± 0.064), while the highest CT value was observed for BGIBMGA006219 at the egg stage (29.022 ± 0.235) (Figure 3(a)).

Figure 3.

Stability Evaluation of Candidate Reference Genes Across All Developmental Stages of Bombyx mori. (a) CT values of nine candidate reference genes and the actin reference gene across the Egg, Instar 1-5, Pupa, and Moth stages. (b-i) Stability of nine candidate reference genes and actin evaluated using delta CT, NormFinder, and BestKeeper algorithms at the Egg, Instar 1-5, Pupa, and Moth stages, indicating the most suitable reference genes at each stage. (j) Comprehensive ranking of the most suitable reference genes across all developmental stages based on delta CT, NormFinder, and BestKeeper algorithms.

The RefFinder results are visible in Table S2. The delta CT and NormFinder results indicated that BGIBMGA006794 (tmbim6) and actin were the most stable reference genes at the egg and first instar stages, whereas BestKeeper identified BGIBMGA000526 and BGIBMGA002429 (cyclophilinA) as the most stable at these stages (Figures 3(b) and 3(c)). At the second and third instar stages, delta CT and NormFinder suggested BGIBMGA007469 (eukaryotic initiation factor 5A) and actin as the most suitable reference genes, respectively (Figures 3(d) and 3(e)).

During the fourth instar stage, delta CT, BestKeeper, and NormFinder recommended BGIBMGA007043 as the optimal reference gene (Figure 3(f)). At the fifth instar stage, delta CT and NormFinder identified BGIBMGA006219 as the most stable, while BestKeeper suggested BGIBMGA007043 (Figure 3(g)). During the pupal stage, actin was universally recommended as the best reference gene (Figure 3(h)). In the moth stage, delta CT recommended BGIBMGA002429 (cyclophilinA), while NormFinder and BestKeeper suggested BGIBMGA010751 and BGIBMGA007469 (Eukaryotic initiation factor 5A), respectively (Figure 3(i)).

When considering all developmental stages, NormFinder recommended BGIBMGA010751, delta CT suggested BGIBMGA007469 (Eukaryotic initiation factor 5A), and BestKeeper identified BGIBMGA002429 (cyclophilinA) as the most suitable reference genes (Figure 3(j)). These findings indicate that the optimal reference genes vary across different developmental stages of Bombyx mori.

The “Comprehensive” ranking combines the results from all three algorithms, with a lower score indicating higher gene stability. BGIBMGA002429 (cyclophilinA), BGIBMGA007469. (eukaryotic initiation factor 5A), BGIBMGA010604 (Ribosomal protein L40), BGIBMGA006794 (tmbim6). The official gene symbols of the above-mentioned genes were retrieved from the Majorbio platform and SilkDB 3.0 database, whereas the remaining genes have not been annotated with definitive gene names to date.

3.4. Validation of selected reference genes across developmental stages

Due to discrepancies in the optimal reference genes recommended by different algorithms across developmental stages, this study selected three genes that remained relatively stable across all stages: BGIBMGA002429 (cyclophilinA), BGIBMGA007469 (eukaryotic initiation factor 5A), and BGIBMGA010751, along with actin. These genes were used to normalize the expression levels of seven target genes across different developmental stages of Bombyx mori. The relative expression levels were calculated with the egg stage as the control (set to 1, indicated by a dashed line in Figure 4), and statistical significance was assessed (Figure 4).

Figure 4.

(a-g) Normalized expression levels of seven target genes using BGIBMGA002429 (cyclophilinA), BGIBMGA007469 (eukaryotic initiation factor 5A), BGIBMGA010751, and actin as reference genes. The dashed line represents the control group (egg stage), where gene expression was set to 1. Official gene symbols for BGIBMGA002350 (aurora kinase B-like), BGIBMGA005841 (mitochondrial ribosomal protein S21), BGIBMGA006280 (cecropin A), BGIBMGA007363 (ribosomal protein S26) were retrieved from the Majorbio platform and SilkDB 3.0 database, whereas the remaining genes have not been annotated with definitive gene names to date.

BGIBMGA001173: The expression patterns remained consistent across the four reference genes, with expression levels exceeding 1. However, significance analysis showed that when normalized to BGIBMGA002429 (cyclophilinA) and BGIBMGA007469 (eukaryotic initiation factor 5A), expression levels in Instar 5, Pupa, and Moth stages were not significantly different (P > 0.05). When normalized to BGIBMGA010751, expression significantly increased in Instar 5 (P < 0.05) but showed no significant changes in Pupa and Moth stages (P > 0.05). In contrast, normalization using actin led to significant differences in Instar 5, Pupa, and Moth stages (P < 0.05) (Figure 4(a)).

BGIBMGA002350 (aurora kinase B-like): Expression patterns were consistent across the four reference genes in all developmental stages except for the Pupa stage, where normalization to BGIBMGA010751 and actin resulted in significantly increased expression (P < 0.05) (Figure 4(b)).

BGIBMGA005098: Expression levels were highly influenced by reference gene selection, and when normalized to BGIBMGA007469 (eukaryotic initiation factor 5A) and BGIBMGA010751, expression levels were lower across all stages. Normalization to BGIBMGA002429 (cyclophilinA) led to significantly higher expression in the Moth stage (P < 0.05), while normalization to actin resulted in significantly higher expression in Instar 5 and Moth stages (P < 0.05) (Figure 4(c)).

BGIBMGA005841 (mitochondrial ribosomal protein S21): Expression trends were consistent across all stages, with no significant differences (P > 0.05) when normalized to BGIBMGA007469 (eukaryotic initiation factor 5A) and BGIBMGA010751. However, normalization using BGIBMGA002429 (cyclophilinA) and actin resulted in expression levels above 1, with a significant increase in the Instar 4 stage (P < 0.05). Additionally, normalization to BGIBMGA002429 (cyclophilinA) showed no significant changes in the Moth stage (P > 0.05), while normalization with actin led to a significant increase in expression at this stage (P < 0.05) (Figure 4(d)).

BGIBMGA006280 (cecropin A): Expression levels in the Instar 2 and Pupa stages were highly influenced by the choice of reference gene. Normalization with BGIBMGA002429 (cyclophilinA) resulted in expression levels above 1, with a significant increase in the Instar 2 stage (P < 0.05). Normalization using actin showed a similar pattern, but it was not statistically significant (P>0.05). When normalized to BGIBMGA007469 (eukaryotic initiation factor 5A) and BGIBMGA010751, expression levels in the Instar 2 and Pupa stages were below 1, with no significant differences (P > 0.05). In the Instar 1, Instar 3, and Instar 5 stages, expression trends were consistent across all four reference genes, with significant differences observed (P < 0.05) (Figure 4(e)).

BGIBMGA007363 (ribosomal protein S26): When normalized to BGIBMGA002429 (cyclophilinA), BGIBMGA007469 (eukaryotic initiation factor 5A), and BGIBMGA010751, expression trends were consistent, with significant downregulation across all developmental stages (P < 0.05). However, normalization using actin resulted in a significant increase in expression in the Instar 5 stage (P < 0.01), while expression significantly decreased in all other stages (P < 0.05) (Figure 4(f)).

BGIBMGA013935: Expression patterns were consistent across all reference genes, with expression levels decreasing in the Instar 1–3 stages, increasing in the Instar 4, Instar 5, and Pupa stages, and significantly increasing in the Moth stage (P < 0.05) (Figure 4(g)).

The selection of reference genes significantly impacts target gene expression across developmental stages, highlighting the importance of precise normalization methods for ensuring the accuracy of RT-qPCR data. This study demonstrates that using a dual-reference gene system can effectively reduce bias resulting from reference gene selection, thereby improving the reliability of gene expressionresults.

3.5. Evaluation of dual reference gene normalization

Using the geNorm algorithm, BGIBMGA006219 and BGIBMGA000526 were identified as the most stable reference gene pair across all the developmental stages of Bombyx mori, with a combined stability score of 0.524. Therefore, these two genes were selected as dual reference genes for subsequent analyses, and their performance was compared with that of actin in normalizing the expression levels of seven target genes across different developmental stages (Figure 5(a)).

Figure 5.

Evaluation of Dual Reference Gene Stability and Expression Normalization. (a) Stability scores of reference genes recommended by the geNorm algorithm for all developmental stages of Bombyx mori. (b-h) Normalized expression levels of seven target genes using BGIBMGA006219, BGIBMGA000526, and actin as reference genes. The dashed line represents the control group, where gene expression in the egg stage is set to 1. The official gene symbols of BGIBMGA002350 (aurora kinase B-like), BGIBMGA005841 (mitochondrial ribosomal protein S21), BGIBMGA006280 (cecropin A), BGIBMGA007363 (ribosomal protein S26) were retrieved from the Majorbio platform and SilkDB 3.0 database, whereas the remaining genes have not been annotated with definitive gene names to date.

Normalization with these dual reference genes resulted in consistent expression patterns for six target genes, BGIBMGA001173, BGIBMGA002350 (aurora kinase B-like), BGIBMGA005841 (mitochondrial ribosomal protein S21), BGIBMGA006280 (cecropin A), BGIBMGA007363 (ribosomal protein S26), and BGIBMGA013953, across all developmental stages (Figure 5(b)–(h)). For BGIBMGA005098, when normalized to BGIBMGA006219, expression decreased significantly during the pupal stage (P < 0.05), whereas normalization with BGIBMGA000526 also resulted in a decrease, but without statistical significance (P > 0.05). However, across the other developmental stages, BGIBMGA005098 exhibited a consistent expression trend following dual-reference normalization (Figure 5(d)).

Normalization solely with actin caused noticeable variations in the expression levels of multiple genes, unlike dual-reference normalization. For example, when normalized to actin, BGIBMGA001173 showed a significant increase in expression during the pupal and moth stages (P < 0.05), whereas dual-reference normalization resulted in no significant changes (P > 0.05) (Figure 5(b)). Similarly, BGIBMGA002350 (aurora kinase B-like) exhibited a significant increase in expression in the pupal stage (P < 0.05) and a non-significant increase in the moth stage (P > 0.05). In contrast, dual-reference normalization resulted in a non-significant increase in the pupal stage (P > 0.05) and a decrease in the moth stage (P > 0.05) (Figure 5(c)).

For BGIBMGA005098, actin-normalized expression significantly increased during the fifth instar and moth stages (P < 0.05), whereas dual-reference normalization led to a decrease in the fifth instar stage and a significant decrease in the moth stage (P < 0.05) (Figure 5(d)). BGIBMGA005841 (mitochondrial ribosomal protein S21) exhibited significantly increased expression in the fourth instar and moth stages (P < 0.05) when normalized to actin, but a reduction in expression under dual-reference normalization (Figure 5(e)). The expression of BGIBMGA006280 (cecropin A) in the pupal stage showed no significant change (P > 0.05) under actin normalization, whereas dual-reference normalization resulted in a significant decrease (P < 0.05) (Figure 5(f)). BGIBMGA007363 (ribosomal protein S26) showed a significant increase in expression in the fifth instar stage (P < 0.05) when normalized to actin, whereas dual-reference normalization led to a significant decrease (P < 0.05) (Figure 5(g)). In the third instar stage, BGIBMGA013953 expression decreased significantly with dualreference normalization (P < 0.05), but not with actin normalization (P > 0.05) (Figure 5(h)).

In this study, three algorithms including ΔCT, BestKeeper and NormFinder were used to evaluate the stability of CT values of genes across all developmental stages of Bombyx mori. The optimal single reference genes screened by the three algorithms were inconsistent, resulting in discrepancies in their ranking results.

Therefore, the single reference genes selected by the above three algorithms were separately applied to calculate the expression levels of target genes in Section 3.4. The results showed that the use of a single reference gene could markedly affect the calculation of gene expression levels and lead to inaccurate quantification results.

On this basis, further analysis was performed in Section 3.5. The geNorm algorithm was adopted to screen and identify BGIBMGA006219 and BGIBMGA000526 as the optimal pair of double reference genes. Comparative verification was carried out between this double reference gene system and the traditional single internal reference actin. The results revealed that the gene expression levels calculated using these two genes jointly as internal references exhibited higher stability and good consistency.

Instead of directly adopting the conflicting ranking results of single reference genes from ΔCT, BestKeeper and NormFinder as the final conclusion, this study reconciled the discrepancies among the three algorithms by introducing and validating a double reference gene correction system. Finally, the combination of BGIBMGA006219 and BGIBMGA000526 was confirmed as an optimal double reference system, which is suitable for the quantification of gene expression at different developmental stages of Bombyx mori.

4. Discussion

With the continuous advancement of insect genome research, Bombyx mori has garnered significant attention as a model insect for physiological and developmental studies, where RNA interference and CRISPR-based gene editing have facilitated functional studies by effectively downregulating target gene expression.^20,21 Nevertheless, qPCR remains an essential tool for gene expression analysis, and its accuracy depends on selecting stable reference genes for data normalization.^1,22 Gene expression in B. mori varies across developmental stages; however, some genes, such as BmTMED6 and BmFKBP12B, exhibit consistently high expression throughout all stages, suggesting their involvement in the entire life cycle.^4–6

During the silkworm’s development-from egg to larva, pupa, and adult-tissue remodelling and cellular renewal pose challenges in identifying reference genes with stable expression across all stages.^3,23,24 In addition, Bombyx mori nucleopolyhedrovirus (BmNPV) and Bombyx mori bidensovirus (BmBDV) are highly prevalent and frequently co-occur in field silkworm populations, causing severe disease and larval mortality, which may further disturb normal physiological status and hamper the selection of stably expressed reference genes. ²⁵ Although actin is commonly used as a reference gene across species, its expression in insects, particularly during developmental transitions, is inconsistent. ¹⁴

In this study, principal component analysis (PCA) and Venn analysis were applied to transcriptome data spanning all developmental stages to identify 125 highly expressed genes with stable expression patterns across allstages. Among these, nine genes with the lowest coefficient of variation (CV), along with the commonly used actin gene, were selected as candidates for stability evaluation. The comprehensive ranking, delta CT, BestKeeper, and NormFinder algorithmsrevealed stage-specific variations in optimal reference gene selection. This highlights the limitations of using a single reference gene across all developmental stages, a trend also observed in honeybee studies. ²⁶ The geNorm algorithm identified an optimal dual reference gene system (BGIBMGA006219 and BGIBMGA000526), which demonstrated more consistent target gene normalization across all stages compared to actin. The gene BGIBMGA006219 has no official gene name, and there are no relevant functional studies or published literature available. According to the annotation from the SilkDB 3.0 database, BGIBMGA000526 is a Bombyx mori Y-box family gene, that encodes a multifunctional nucleic acid-binding protein belonging to the cold shock domain (CSD) protein family, which is extensively involved in core biological processes including transcriptional regulation, translational control, RNA stability maintenance, and DNA repair. Recent research showed that the silkworm Y-box protein (YBP) interacts with C/EBPγ to activate Dnmt1 transcription and regulate embryonic development and ovarian function. ²⁷ Housekeeping genes play fundamental roles in cellular function, for example, GAPDH ensures glycolytic flux through the EntnerDoudoroff pathway, ²⁸ while ribosomal and cytoskeletal proteins, such as actin and tubulin, are consistently expressed across various cell types and physiological conditions. ²⁹ Consistent with previous findings, transcriptome analysis identified 125 candidate reference genes, primarily enriched in translation, metabolism, and ribosome-related pathways. Ribosomal proteins, which are highly conserved across eukaryotes, are integral to intracellular protein translation ³⁰ and prior studies have shown that ribosome-associated genes exhibit minimal variation, making them reliable reference genes for RT-qPCR normalization in both mammals and insects.^31,32

The four algorithms found no single gene stably expressed throughout development. However, based on overall CT values, BGIBMGA002429 (cyclophilinA), BGIBMGA007469 (eukaryotic initiation factor 5A), and BGIBMGA010751 were identified as the most suitable reference genes. BGIBMGA002429 encodes Cyclophilin A, a molecular chaperone involved in protein folding, apoptosis, and immune regulation.^33,34 Complete metamorphosis in insects requires the synthesis of structural proteins such as chitin, ³⁵ and Cyclophilin A likely aids in the proper folding of these proteins, facilitating normal tissue development. ³⁶ Moreover, as apoptosis serves as a primary defence mechanism in insects lacking adaptive immunity, ³⁷ Cyclophilin A interacts with apoptosis-inducing factor (AIF) to regulate cell death, underscoring its role in innate immunity. ³⁸ BGIBMGA007469 encodes eukaryotic translation initiation factor 5A (eIF5A), which also participates in apoptosis regulation. Translation initiation factors have been used as reference genes in previous studies,^16,39,40 and BGIBMGA010751 is homologous to transcription factor BTF3, which is involved in transcription initiation and apoptosis regulation. ⁴¹ Consistent with earlier findings, eIF5A and BTF3 have been validated as reference genes, ⁴² while Cyclophilin A may exhibit variable expression under pathogen challenge, necessitating further validation under different experimental conditions.

Notably, actin, frequently used as a reference gene, exhibited instability across developmental stages (CV = 0.896). This variability may result from tissue remodelling and cellular renewal during metamorphosis, leading to fluctuations in actin expression,^43–45 as seen in previous studies. ⁴⁶

Thus, to address the limitations of single reference genes, the geNorm algorithm identified a dual reference gene system (BGIBMGA006219 and BGIBMGA000526) that demonstrated stable expression across all developmental stages, ensuring consistent normalization of the seven target genes. For instance, BGIBMGA007363 (ribosomal protein S26) exhibited increased expression when normalized to actin, but showed reduced expression at the fifth instar stage when normalized using the dual reference system. Similar discrepancies were observed for BGIBMGA005841 (mitochondrial ribosomal protein S21) and BGIBMGA006280 (cecropin A). These findings align with previous research, which shows that reliance on a single reference gene can introduce significant errors in RT-qPCR analysis. ⁴⁷ Finally, adoption of a dual reference gene system adheres to MIQE guidelines, which advocate the use of multiple reference genes to minimize technical variability. ⁴⁸

5. Conclusion

This study reveals stage-specific variations in reference gene stability during silkworm development and emphasizes the requirement of a dual reference gene system for reliable RT-qPCR normalization. The paired reference strategy detailed herein effectively reduces systematic bias due to by single reference genes, greatly enhancing the accuracy and repeatability of quantification. Current results establish a robust methodological foundation for gene expression analysis in Bombyx mori and other insects. Further validation of this dual reference system under diverse biotic and abiotic stresses will strengthen its application and promote progress in insect molecular biology research.

Supplemental material

Supplemental material - Selection and evaluation of dual reference genes in silkworm (Bombyx mori) based on transcriptome data

Supplemental material for Selection and evaluation of dual reference genes in silkworm (Bombyx mori) based on transcriptome data by Lindan Sun, Jingqi Xv, Binbin Sun, Xvsheng Jin, Yi Jin, Min Tang, Keping Chen in Science Progress