Corona or hydrochloric acid modulates embryonic diapause in silkworms by activating different signaling pathways

Silkworm materials

In this study, the silkworm strain 7532 was bred by the Guangxi Zhuang Autonomous Region Sericulture Technology Promotion Station. Domesticated silkworms were reared until moths hatched. Male and female moths were allowed to self-cross for 4 h, after which female moths were placed on silkworm egg paper to lay eggs. The eggs were maintained at 25 °C for three days and then stored at 5 °C for further analysis.

Silkworm egg treatment

After 50 days, silkworm eggs were removed from refrigeration and left at room temperature for 4 h before treatment. The eggs were divided into three groups (each group contains 120–150 eggs): (1) acid treatment at 47.2 °C for 5 min and 30 s (HCl relative density = 1.092) (Zhao et al., 2012), (2) corona treatment (12 KV, 1 min) (Zhang et al., 2022), and (3) control, which received no treatment. Samples were collected at 1, 6, and 20 h after treatment (HAT) and immediately frozen in liquid nitrogen for further analysis. Two replicates were prepared for each treatment group.

RNA extraction, library construction, and sequencing

Total RNA was extracted using a TRIzol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) via RNAse-free agarose gel electrophoresis. mRNA was enriched using oligo(dT) beads, fragmented, and reverse-transcribed into cDNA using NEBNext Ultra RNA Library Prep Kit for Illumina (NEB #7530; New England Biolabs, Ipswich, MA, USA). The cDNA fragments were purified, ligated with Illumina sequencing adapters, and size-selected using agarose gel electrophoresis. The ligation products were then subjected to polymerase chain reaction (PCR) amplification, and the resulting cDNA library was sequenced using the Illumina NovaSeq 6000 platform.

RNA-seq data analysis

Adapter sequences and low-quality reads were removed from each dataset using fastp (version 0.18.0) (Chen et al., 2018). An index of the reference genome was constructed, and paired-end clean reads were mapped to this reference genome using HISAT2 (version 2.2.4) via the “-rna-strandness RF” option and other default parameters (Kim et al., 2019). Clean reads were then used for gene assembly and abundance calculations. The silkworm reference genome and annotation files were downloaded from the National Center for Biotechnology Information genome database (accession no. ASM15162v1) (Xia et al., 2004). The fragments per kilobase of transcript per million mapped reads (FPKM) values were calculated via StringTie (version 1.3.1) and HISAT2 (Langmead & Salzberg, 2012; Pertea et al., 2016). The statistical power of this experimental design, calculated in RnaSeqSampleSize was 0.84 (https://cqs-vumc.shinyapps.io/rnaseqsamplesizeweb/) (Zhao et al., 2018). Differentially expressed gene (DEG) analysis was conducted using DESeq (Love, Huber & Anders, 2014), and transcripts showing an absolute log2 (fold change) value of > 1 and a p-value of < 0.05 were considered differentially expressed. Time series expression patterns of the three groups were analyzed using Mfuzz (Kumar & Futschik, 2007). DEGs were further subjected to Gene Ontology functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Boyle et al., 2004; Kanehisa et al., 2007; Kanehisa & Goto, 2000). RNA-seq data are available at NCBI BioProject (accession number PRJNA1165327).

Validation of gene expression patterns

Ten genes were selected for validation using quantitative reverse transcription-PCR (qRT-PCR). First-strand cDNA was synthesized using a TUREscript 1st Strand cDNA Synthesis Kit (Aidlab, Beijing, China). Gene-specific primers are listed in Table S1. Actin3 was used as an internal control. qRT-PCR was conducted using 2 × SYBR^® Green Master Mix (DBI, Langsing, MI, USA) and IQ5 (Bio-Rad, Hercules, CA, USA) with a cycling temperature of 60 °C and a single peak on the melting curve to obtain a single product according to the manufacturer’s instructions. The 20 µL reaction volume consisted of forward and reversed primers (1 µL), 2 × SYBR^® Green Master Mix (10 µL), ddH2O (6 µL), and cDNA (2 µL). Three technical replicates were prepared for each gene. Relative gene expression levels were calculated using the 2^−ΔΔCt method and visualized graphically. All the numerical data in figures are presented as mean ± standard deviation (SD) of three independent experiments. The obtained data were subjected to unpaired a two-tailed Student’s t-tests using GraphPad Prism software (version 8).

Effects of corona and HCl treatments on diapause termination in silkworm eggs

The results showed that silkworm eggs began hatching on the ninth day (Fig. 1A). By the tenth day, over 80% of the corona-treated eggs had hatched, whereas approximately 5% of the HCl-treated eggs had hatched. Moreover, none of the untreated eggs had hatched. Although the hatching rate of HCl-treated eggs eventually reached that of corona-treated eggs, the incubation period for corona-treated eggs was relatively shorter (Fig. 1B). To explore the molecular mechanisms through which corona treatment terminates diapause, we conducted further transcriptomic analyses of early-stage corona-treated silkworm eggs.

RNA-seq of corona-treated, HCl-treated, and control groups

RNA samples from the control and corona-treated groups were collected at 1, 6, and 20 HAT. For comparison, RNA samples from the HCl-treated group were collected at 6 and 20 HAT. These samples were defined in Table 1 and sequenced using Illumina sequencing (Table 1). In total, approximately 752 million raw reads were obtained, with an average of 47 million reads per sample (Table S2). After filtering, >94% of the reads were classified as clean (Table S3). The clean reads were mapped to the silkworm reference genome, with a mapping rate of > 84% (Table S4). Over 64% of the mapped reads corresponded to genes, with > 95% mapped to exons and approximately 30% mapped to intergenic regions (Table S5).

DEGs in corona-treated, HCl-treated, and control groups

First, we analyzed the correlation among corona-treated, HCl-treated, and control samples. The results showed a high degree of correlation between corona-treated samples at 1 h and control samples (at 1, 6, and 20 h), indicating that there was almost no difference between the corona-treated eggs of 1 HAT and control samples. E6 showed higher similarity to H20 and E20 than to H6 and CK20. The high similarity between E20 and H20 indicated that both the corona and HCl treatments after 20 h can activate the expression of genes required for embryonic development, potentially initiating the development process (Fig. 2A).

Comparisons between the treated and control groups showed relatively low numbers of DEGs in CK1 vs. E1 (109 DEGs, with 70 upregulated and 39 downregulated genes) and H20 vs. E20 (62 DEGs, with 34 upregulated and 28 downregulated genes) (Tables S6 and S7). However, at 6 HAT, 1,168 DEGs (502 upregulated and 666 downregulated genes) were observed in CK6 vs. E6 and 809 DEGs (503 upregulated and 306 downregulated genes) in CK6 vs. H6 (Tables S8 and S9). At 20 HAT, 1,116 DEGs (468 upregulated and 648 downregulated genes) were detected in CK20 vs. E20 and 1,351 DEGs (531 upregulated and 820 downregulated genes) in CK20 vs. H20 (Tables S10 and S11). The similarity between E6 and H20 or E20 was even higher than that between E6 and H6, indicating substantial differences between E6 and H6. In particular, 947 DEGs (293 upregulated and 654 downregulated genes) were observed in H6 vs. E6 (Fig. 2B).

Functions of DEGs in multiple comparison groups

As samples treated after 1 h and control samples showed minimal differences and the control group showed limited variation across the three time points, our study focused on DEGs between the treated groups (at 6 and 20 h) and the control group. Functional analysis of DEGs in the four comparison groups showed that DEGs in CK6 vs. E6 were mainly enriched in FoxO signaling pathway and vitamin B6 metabolism (Fig. 3A). Conversely, DEGs in CK6 vs. H6 were mainly enriched in pathways related to longevity regulation, protein processing in the endoplasmic reticulum, and arginine and proline metabolism (Fig. 3B). At 20 HAT, DEG enrichment in CK20 vs. E20 and CK20 vs. H20 showed certain similarities. The main enriched pathways in both comparison groups included glycolysis/gluconeogenesis, vitamin B6 metabolism, and pentose phosphate pathway (Figs. 3C and 3D).

Identification of DEGs involved in corona and HCl treatments

Gene function enrichment analysis revealed differences in gene expression were observed between corona and hot HCl treatments at 6 h. A comparison of DEGs at multiple time points between the treated and control groups showed only nine common DEGs across five comparison groups and 129 common DEGs across four comparison groups (Figs. 4A and 4B). This suggests that corona and HCl treatments may modulate different pathways to promote silkworm embryonic development.

At 6 HAT, DEG analysis showed that two comparison groups—CK6 vs. E6 and CK6 vs. H6—shared 295 DEGs, and their functions were mainly related to cellular senescence. The CK6 vs. E6 comparison group had 873 specific DEGs, which were mainly enriched in purine metabolism, whereas the CK6 vs. H6 comparison group had 514 specific DEGs, which were mainly related to the longevity regulating pathway and protein processing in the endoplasmic reticulum (Fig. 4C).

At 20 HAT, the number of common DEGs between CK20 vs. E20 and CK20 vs. H20 increased to 888, with their functions enriched in the Toll and lmd signaling pathways and insect hormone biosynthesis. The CK20 vs. E20 comparison group had 228 specific DEGs, which were mainly associated with the p53 signaling pathway, whereas the CK20 vs. H20 comparison group had 463 specific DEGs, which were enriched in the purine metabolism pathway (Fig. 4D).

FoxO signaling pathway response to corona treatment

After corona treatment, the FoxO signaling pathway was notably enriched. Corona treatment likely induces DNA damage, which may lead to the upregulation of CDK2 and Ras1. Ras1 promotes the expression of FoxO genes, whereas CDK2 has been shown to inhibit the expression of downstream FoxO genes. The interplay between these factors may result in a complex regulatory effect on the overall expression of FoxO genes (Figs. 5A and 5B). FoxO genes ultimately influence the expression levels of cell cycle–related genes, such as Cyclin B2/B3 (significantly upregulated after treatment) and Cyclin G2 (downregulated after treatment) genes (Figs. 5A and 5B). These results suggest that corona treatment activates cell cycle regulation through the FoxO signaling pathway. Further validation of corona-treated samples using qRT-PCR showed that the expression patterns of most genes, except for Ras1, were consistent with the RNA-seq data. In addition, several genes in the control group showed changes in expression levels across the three time points, which may be related to the process of adjusting the silkworm eggs to room temperature after removing them from refrigeration (Fig. 5C).

Heat shock proteins and insect hormone biosynthesis may be involved in embryonic development via HCl treatment

The results showed that the enriched pathways in the HCl-treated group differed significantly from those in the corona-treated group. Notably, HCl treatment increased the expression of multiple heat shock protein (HSP)-related genes, potentially activating downstream pathways such as longevity regulation and ER-associated degradation (Figs. 6A and 6B). In addition, hormone-related pathways may be activated, with upregulated expression of ALDH and CYP18A1 possibly promoting the production of juvenile and molting hormones (Figs. 6A and 6B). Examination of HCl-treated samples through qRT-PCR showed that the expression levels of HSP70 and ALHD were significantly upregulated at 1 and 6 HAT. Conversely, the corona-treated group showed relatively minimal changes. Similarly, the expression levels of sHSP20 and CYP18A showed significant upregulation at 6 HAT compared with those at 1 HAT. Fluctuations in gene expression levels in the control group were also observed, possibly due to temperature changes after removing the eggs from refrigeration (Fig. 6C). Although 1-h data from RNA-seq were unavailable, gene expression trends were generally consistent with qRT-PCR results.