Genome-wide analysis reveals transcriptional and translational changes during diapause of the Asian corn borer (Ostrinia furnacalis)

RNA-seq of ACB head at different stages of diapause revealed a DEGs profile at transcriptional level

To better understand genome-wide changes of transcription and translation in ACB at various stages of diapause, we performed Ribo-seq and RNA-seq to identify gene expressions in ACB at various stages of diapause (Additional file 1: Fig. S1A). After quality control, the RNA-seq had an average Q20 score above 97% and an average Q30 score above 92% (Supplementary Table 5). The average alignment rate was 85%, indicating that we obtained sufficient and high-quality data for subsequent analyses. First, we performed RNA-seq to detect differentially expressed genes at the transcriptional level. The gene density plots of the nine samples used in RNA-seq indicated that the main peak positions in each state are roughly the same, with the pre-diapause state having the most concentrated gene expression density distribution (Additional file 1: Fig. S1B). The principal component analysis (PCA) of the RNA-seq data showed a clear separation between groups (Additional file 1: Fig. S1C). The reliability of the data performance was evaluated through the repeatability study represented by the heat map, which showed the high correlation among the three biological replicates of each sample (Additional file 1: Fig. S1D). Transcriptome alterations were determined through differential expression analysis of RNA-seq data. Comparative analysis of the transcriptome of pre-diapause, diapause, and no-diapause samples of ACB showed DEGs between these samples. Comparing with the no-diapause sample, there are 974 significantly changed genes in the diapause sample, with comparable numbers of genes with increased abundance and genes with decreased abundance (Fig. 1A, D). When compared to pre-diapause period, there are 1455 DEGs in the diapause sample, with more than twice as many at genes with decreased abundance compared with genes with increased abundance (983 vs. 472) at diapause period (Fig. 1B, D). Comparing with the no-diapause sample, there are 201 significantly changed genes in the pre-diapause sample (log2(fold change), adjusted p-value < 0.05), including 168 genes with increased and 33 with decreased abundance (Fig. 1C, D). Comparative analysis showed that 35 genes were significantly changed between every two of the three samples (Fig. 1D).

Fig. 1figure 1

RNA-seq revealed gene profiles that are changed at transcriptional level. AC Volcano plot representing the DEGs of RNA-seq candidates with the threshold of |log2 (fold change) |≥ 1 and adjusted p-value < 0.05 were regarded as significant differentially expressed genes. The genes with significant increased and decreased abundance were marked by orange and blue, respectively. D The upset plot represents the intersection of significantly differentially expressed genes identified in the RNA-seq analysis across three comparisons for AC. The bars indicate the number of significantly differentially expressed genes in each intersection set, with the horizontal bars showing the total number of significantly differentially expressed genes for each comparison. E, F Compare GO and KEGG enrichment analysis of DEGs revealed by RNA-seq. The circle color indicated adjusted p-value of the enrichment and the dot size indicated the number of DEGs in the functional class or pathway. G Clustering and functional enrichment analysis of genome-wide gene expression at transcriptional levels. The heatmap show the eight profiles of the clustered genes, with the number of genes for each cluster were noted in the box just right to the heatmap. The genes in each cluster were presented in the heat map, which were normalized by Z-score. The representative trends of the expression of each cluster were shown in the mid box. The top five GO enrichment term of biological progress and KEGG enrichment pathway are shown to the right of each cluster. The text size was changed according to the p-value (the larger the font, the more significant the pathway)

To understand the biological implications underlying these differences, GO enrichment and KEGG enrichment analysis were carried out. GO enrichment analysis showed that the mitosis related term (such as mitotic spindle organization, microtubule cytoskeleton organization involved in mitosis, spindle midzone), negative chemotaxis, and synaptic target inhibition were specifically enriched in DEGs between diapause and no-diapause samples. The mitochondrial energy metabolism related terms (such as generation of precursor metabolites and energy, ATP metabolic process, mitochondrial protein-containing complex) were specifically enriched in DEGs between diapause and pre-diapause samples. It is interesting to note that the cuticle development, chitin-based cuticle development, and organic acid biosynthetic process were all enriched between each pair of the three samples (Fig. 1E). Specifically, 34, 52, and 13 genes were enriched in the cuticle development process in the comparisons between diapause vs. no-diapause, and diapause vs. pre-diapause and pre-diapause vs. no-diapause, respectively. In the comparisons between diapause and no-diapause, as well as between diapause and pre-diapause, the majority of these genes showed significantly decreased expression. Conversely, in the comparison between pre-diapause and no-diapause, most of these genes (495 of 947 in diapause vs no-diapause; 983 of 1455 in diapause vs pre-diapause; 33 of 210 in pre-diapause vs no-diapause) exhibited significantly increased expression (Additional files 3 and 4, Additional file 1: Fig. S2A-C). The KEGG enrichment analysis indicated that DNA replication, FoxO signaling pathway, and Wnt signaling pathway were specifically enriched between diapause and no-diapause samples. Purine metabolism and Insect hormone biosynthesis was enriched between each pair of the three samples (Fig. 1F). Specifically, 6, 11, and 5 differentially expressed genes were significantly enriched in the insect hormone biosynthesis process in the comparisons between diapause vs. no-diapause, diapause vs. pre-diapause, and pre-diapause vs. no-diapause, respectively. The majority of these genes exhibited decreased expression (Additional files 3 and 4, Additional file 1: Fig. S3A-C). All the gene lists used for analysis were included in Additional file 4. The specific gene lists of the enrichment analysis were included in Additional file 3.

To gain deeper insights into the temporal expression patterns of genes across the three time points represented by the samples in our study, we performed the gene temporal expression analysis by clustering and visualizing the time series of the gene expression data from the RNA-seq. The fuzzy K-means cluster method from Mfuzz package was used to obtain appropriate 8 clusters (Fig. 1G). For cluster 1, there were 818 genes tended to be up regulated from no diapause to pre-diapause to diapause and were primarily enriched in the organic acid biosynthetic process, peroxisome. A similar trend was found in cluster 2. Clusters 5 and 6 exhibit a common gene expression trend characterized by an initial decrease followed by a gradual increase. The genes significantly enriched in these clusters are mainly associated with biological processes including non-coding RNA metabolism, intracellular molecular biology of DNA, RNA, and protein synthesis, transport, repair, and degradation. In cluster 3, there is an initial and substantial decrease in gene expression levels, followed by a sustained, unchanging pattern. These genes are associated with cellular pathways related to proteasomes and ubiquitin-mediated degradation. GO enrichment highlights proteasomal protein catabolic process, protein modification by small protein removal, and mRNA processing. KEGG enrichment includes proteasome, spliceosome, nucleocytoplasmic transport, polycomb repressive complex, and ubiquitin-mediated proteolysis pathways. These suggest roles in proteolysis, protein modification, and RNA processing. We also found that transcript abundance of genes in clusters 7 and 8 initially increases and then decreases. GO enrichment analysis for genes in cluster 7 shows involvement in processes such as chitin-based cuticle development, purine metabolism, and various metabolic processes. KEGG enrichment analysis indicates these genes are involved in pathways related to steroid biosynthesis, lipoic acid metabolism, neuroactive ligand-receptor interaction, amino acid degradation, and fatty acid metabolism. The transcript abundance of genes in cluster 8 initially increases and then decreases. GO enrichment analysis shows that these genes are involved in energy production and metabolism processes, including ATP metabolism and oxidative phosphorylation. KEGG enrichment analysis indicates involvement in major metabolic pathways such as oxidative phosphorylation, carbon metabolism, and fatty acid degradation. These results suggest that the genes in cluster 8 play a crucial role in energy production and overall metabolic regulation.

Similarly, we employed Gene Set Enrichment Analysis (GSEA) to investigate RNA-seq data from different diapause stages. Consistently, our analysis revealed the activation and suppression of pathways related to many cellular processes during the transition from the no-diapause state to diapause (Additional file 1: Fig. S4A-R). This indicates a concerted downregulation of these pathways as part of the diapause onset process. We found that numerous gene groups, including those related to cuticle development and proteasomes exhibit significant dynamic changes (Additional file 1: Fig. S4A-R). For instance, pathways related to cuticle development, oxidative phosphorylation, and glycolysis/gluconeogenesis were significantly activated during the pre-diapause stage, with a marked decrease of many core genes compared to the no-diapause state. However, when transitioning to the full diapause stage, these pathways were significantly inhibited (Additional file 1: Fig. S4A, E, F, G, H, K, L), with a substantial decrease of the core genes relative to the no-diapause state. This trend is also evident when compared to the pre-diapause stage (Additional file 1: Fig. S4M, N, Q, R). This indicates that the cuticle development, oxidative phosphorylation, and glycolysis/gluconeogenesis pathways exhibit a similar dynamic pattern: significant increase followed by a decrease, with the final levels being lower than the no-diapause state. Other pathways, such as the proteasome (Additional file 1: Fig. S4J, P) and mitochondrial energy metabolism-related pathways (Additional file 1: Fig. S4C), were also significantly enriched, corroborating the previous findings from the differentially expressed genes’ GO/KEGG analysis and the clustering results.

These gene dynamics are consistent with the trends observed in the clustering profiles. Observing these trends from multiple perspectives, including significantly differentially expressed genes and dynamically changing genes, validates the accuracy of the observed changes.

Ribo-seq of ACB head at different stages around diapause revealed a DEGs profile at translational level

To understand translation changes during diapause, we used Ribo-seq to identify ribosome-protected footprints (RPFs). The Ribo-seq data had an average Q20 score above 99% and an average Q30 score above 97%, indicating that we obtained high-quality data for subsequent analyses (Additional file 4). Analysis showed the dominant RPF length was 28 nucleotides (nt) across all stages (Additional file 1: Fig. S5A-C). These 28-nt RPFs had higher counts in the dominant reading frame, ensuring significant 3-nt periodicity, a key quality control measure (Additional file 1: Fig. S5A-C). Metagene profiles of RPF counts near translation initiation (TIS) and termination sites (TTS) confirmed proper P-site estimation. During diapause stages, read length distribution consistently showed 28-nt RPFs (Additional file 1: Fig. S5A-D). Pre-diapause and diapause stages had broader RPF length distributions, suggesting ribosomal differences. P-site analysis indicated that 28-nt RPFs predominantly occupied frame 0 of the CDS, over 95%, showing preferential translation within the CDS (Additional file 1: Fig. S5E). Lower P-site signals in 5′ and 3′ UTRs suggested weaker translation activity in these regions (Additional file 1: Fig. S5F).

We compared ORFs and found that over 90% annotated across all conditions. Actively translated ORFs were more common during no-diapause and diapause stages than pre-diapause (Additional file 1: Fig. S5G, H). Meta-gene heatmaps showed distinct three-nucleotide periodicity, emphasizing precise ribosomal positioning (Additional file 1: Fig. S5G). This highlighted dynamic changes in ribosomal activity during diapause, affecting protein synthesis regulation (Additional file 1: Fig. S5I, J).

Comparative analysis of the translation levels of pre-diapause, diapause and no-diapause samples of ACB also showed DEGs between these samples (Fig. 2A, B, C). Compared with the no-diapause sample, there are 383 significantly changed genes in the pre-diapause sample (log2(fold change) ≥ 1, adjusted p-value < 0.05), including 297 genes with increased abundance and 86 downregulated genes with decreased abundance (Fig. 2D). Comparing with the no-diapause sample, there are 1769 significantly changed genes in the diapause sample. When compared to pre-diapause period, there are 1538 DEGs in the diapause sample, the translation of most genes was decreased (647 genes with increased abundance vs. 891 genes with decreased abundance) at diapause period (Fig. 2D). Comparative analysis showed that 134 genes were significantly changed between every two of the three samples (Fig. 2D). To understand the biological processes these genes were involved in, we performed GO and KEGG enrichment analyses. The gene lists used for analysis were included in Additional file 4. The specific gene lists of the enrichment analysis were included in Additional file 3. GO enrichment analysis showed a bunch of terms that were similarly enriched in the three samples. It mainly involves cuticle development and metabolism as well as metabolic processes such as ribonucleotide and amino acids metabolism etc. (Fig. 2E). The KEGG analysis showed that these genes were enriched in different pathways. For example, translation-related genes were mainly enriched in oxidative phosphorylation, Toll and Imd signaling pathway, glycolysis and gluconeogenesis, lysosome, and purine metabolism, when comparing the diapause sample to the no-diapause or pre-diapause samples (Fig. 2F). When compared, the no-diapause and pre-diapause samples, the pathways of drug metabolism, biosynthesis of nucleotide sugars, pyrimidine metabolism, etc., are enriched. The three pathways of glycine, serine and threonine metabolism, amino sugar and nucleotide sugar metabolism, and ascorbate and aldarate metabolism are enriched in differentially translated genes between every two of these three gene lists (Fig. 2F). These results indicate that the translational profiles changed during the different stages of diapause.

Fig. 2figure 2

Ribo-seq revealed gene profile that are changed at translation level. AC Volcano plot representing the DEGs of Ribo-seq Candidates that satisfied the criteria of |log2 (fold change) |≥ 0.265 and adjusted p-value < 0.05 were regarded as significantly expressed genes, which are marked by orange and blue for up and downregulation, respectively. D The upset plot represents the intersection of significantly differentially expressed genes identified in the Ribo-seq analysis across three comparisons for AC. The bars indicate the number of significantly differentially expressed genes in each intersection set, with the horizontal bars showing the total number of significantly differentially expressed genes for each comparison. E, F Compare GO and KEGH enrichment analysis of DEGs in revealed by Ribo-seq. The circle color indicates the enrichment adjusted p-value, and the dot size indicates the number of DEGs in the functional class or pathway. G Clustering and functional enrichment analysis of genome-wide gene expression at translational levels. The heatmap show the eight profiles of the clustered genes, and the number of genes for each cluster were noted in the box just right to the heat map. The genes in each cluster were presented in the heat map, which were normalized by Z-score. The representative trends of the expression of each cluster are shown in the mid box. The top five GO enrichment term of biological progress and KEGG enrichment pathway are shown to the right of each cluster. The text size corresponds to the p-value (the larger the font, the more significant the pathway)

We also conducted a cluster analysis of the change patterns of translation levels genes in the three time points represented by the three samples (Fig. 2G). With the development of diapause, the expression of 1300 genes in cluster1 decreased at the pre-diapause stage and further decreased during diapause. Enrichment analysis showed that these pathways involve various metabolic processes, lipid synthesis, subcellular functions, molecular transport, and neural signal transduction during diapause. Similar patterns of change were observed in 1922 genes associated with encompassing various aspects of energy production, organic molecule synthesis, breakdown, and cellular functions in cluster 2. There are also two gene sets with opposite change trends. For example, genes are primarily involved in various vital biological functions, including protein processing, cellular uptake, signal transduction, photoreception, metabolism, immune response, intracellular molecular transport, and membrane processes. These genes low expressed in pre-diapause process but are highly expressed in no-diapause and diapause process in cluster 6. Cluster 3 exhibits a similar trend. Additionally, many clusters with similar dynamic trends and enrichment analysis results show a high degree of similarity with the previous RNA-seq results, such as cluster 3 (Fig. 2G) and cluster 5 (Fig. 1G).

Combinational analysis of the RNA-seq and Ribo-seq data revealed differences in translational efficiency between different samples

Translation efficiency (TE) is an important aspect of translation process, which standardizes the change in translation level according to transcription level. The genes with significant RPF and mRNA but not significant TE belong to “Forwarded”. A total of 834 genes were “Forwarded” genes when comparing diapause vs. no-diapause samples (Fig. 3A). Similarly, diapause vs. pre-diapause and pre-diapause vs. no diapause had 956 and 102 such genes, respectively (Fig. 3B, C). The second group of genes was classified as “Exclusive”, whose change in RPF was not driven by changes in mRNA (yellow in Fig. 3A–C). Therefore, for these genes TE and RPF had significant changes, but mRNA changes were not significant, which were called differentially translation efficiency genes (DTEGs). There were 155, 103, and 8 in each of the three comparisons (Fig. 3A–C). For example, we found glycolysis/gluconeogenesis related genes showed significantly decreased abundance only at translational level between the diapause and no-diapause by comparative enrichment analysis (Fig. 3D, E bottom right). The third group of genes was classified as “Intensified” that were regulated by both transcription and translation (significantly changed in mRNA, RPFs, and TE) (purple in Fig. 3A–C). In this case, mRNA changes acted with the change in TE. There were 28, 34, and 2 in each of the three comparisons. In this regard, we found ORFs encoding fatty acid synthase activity increases or decreases in both mRNA levels and translational efficiency between the diapause and no-diapause (Fig. 3D, E top left). Finally, the fourth group of genes was classified as “Buffered” whose transcriptional changes were offset by TE changes (blue in Fig. 3A–C). The transcriptional change and TE change of these genes were in the opposite direction. Thus, genes with significant TE and mRNA but not significant RPF are also considered translation buffered. There were 81, 141, and 18 genes in the three comparisons (Fig. 3A–C). In this regard, we found, for example, that developmental decreases in drug metabolism gene RNA levels were effectively buffered by corresponding increases in TE between the diapause and no-diapause (Fig. 3D, E top right).

Fig. 3figure 3

Genome-wide transcriptional and translational regulation at different stages around diapause. AC Scatter plot of fold changes between different stage of diapause for all ORFs in Ribo-seq data and the corresponding gene in RNA-seq data. The genes with DTGs and DTEGs are displayed. D Dot plot of the top enriched KEGG pathway in each regulatory category defined in A. The circle color indicates the enrichment adjusted p-value, and the dot size indicates the generation in the pathway. E Heat map of genes associated with the top KEGG enrichment pathway in each regulatory category identified in A

We also performed the GO enrichment analysis for genes categorized by transcription, translation, translation efficiency, and deltaTE differences between diapause and no-diapause (Additional file 1: Fig. S6A). Furthermore, we conducted functional enrichment analysis for differential expressed genes between diapause and pre-diapause (Additional file 1: Fig. S6B–C) and between pre-diapause and no-diapause (Additional file 1: Fig. S6D–E). Through comparative analysis, we identified several common and unique terms enriched in various gene groups. For instance, genes related to GO terms like vitamin binding, fatty acid synthase activity, and fatty acid elongase activity, as well as KEGG pathways like fatty acid metabolism and fatty acid elongation, were intensified in both diapause vs. no-diapause and diapause vs. pre-diapause comparisons, indicating the significant and widespread role of fatty acid synthesis metabolism in later stages of diapause. Conversely, genes related to the Toll and Imd signaling pathway, oxidative phosphorylation, and lysosome were categorized as forwarded in these comparisons, with lysosome specifically showing significant enrichment at the RNA level in pre-diapause vs. no-diapause, and in the RPF and forwarded groups in the other comparisons. The gene lists used for analysis were included in Additional file 4. The specific gene lists of the enrichment analysis were included in Additional file 3.

The lncRNAs profile were changed in the diapause process

To investigate the lncRNAs during diapause of ACB, we analyzed the transcriptome data of ACB in no-diapause, pre-diapause and diapause stages, and systematically identified and analyzed the lncRNAs of ACB (Fig. 4A). In total, 4561 high-confidence ACB lncRNAs were selected through four different lncRNA prediction software (Fig. 4B). Furthermore, by removing the already annotated lncRNAs in the insectbase2 database, 318 new high-confidence ACB lncRNAs were identified. All subsequent lncRNA analyses were conducted based on these 318 newly annotated and 7349 known lncRNAs. To characterize their genomic features, we compared the putative ACB lncRNAs with known protein-coding mRNAs. Generally, the number of exons in the ACB lncRNAs transcripts is mostly between 2 and 3 (Fig. 4C), and the transcript length is within 2500 bp (Fig. 4D), with a significantly lower GC content compared to mRNA (Fig. 4E).

Fig. 4figure 4

Comprehensive analysis of LncRNAs identification and differential expression in ACB developmental stages. A Pipeline of lncRNAs identity. This part depicts the workflow for the identification of lncRNAs in ACB. B Venn diagram demonstrating the overlap and unique sets of lncRNAs identified through the utilization of filtering software, including CPC2, CAPT, CNCI, and Pfam, subsequent to the initial identification step. C Bar chart visualizes the distribution of exon counts within the identified lncRNAs. The blue bars represent known lncRNAs, the red bars denote novel lncRNAs, and the green bars signify mRNAs. D Density plots revealing the distribution of transcript sizes for known lncRNAs (in blue), novel lncRNAs (in red), and mRNAs (in green). E Box plots of the GC content for known lncRNAs, novel lncRNAs, and mRNAs. FH The volcano plot shows the differential expression analysis results of lncRNAs among the three stages of ACB: no-diapause, pre-diapause, and diapause. The x-axis represents the log2(fold change) value, and the y-axis represents the -log10(adjusted p-value). Red dots indicate significantly upregulated lncRNAs, blue dots indicate significantly downregulated lncRNAs, and grey dots indicate lncRNAs with no significant differential expression. The number of upregulated and downregulated lncRNAs in each comparison is shown in the plot. I The relationship network among lncRNAs, target genes, and KEGG pathways of diapause vs. pre-diapause. The color represents the significance value of lncRNAs differential analysis or mRNA differential analysis. Different shapes represent different molecular types. The size of the shape proportionally maps the log2(fold change) of lncRNAs differential analysis or mRNA differential analysis, with positive values indicating a significant increase during diapause and negative values indicating a significant increase during pre-diapause stage. JL Validation of the relative expression levels of lncRNA using qPCR. Relative expression levels of lncRNA ofu_lnc_1066, ofu_lnc_0096 and ofu_lnc_004514 in pre-diapause vs. no-diapause (from left to right mean ± SEM: 1.0716 ± 0.29626, 3.7935 ± 0.22433, 1.0351 ± 0.17751, 1.5516 ± 0.198, 1.0907 ± 0.30674, 0.61699 ± 0.17486, n = 3). K Relative expression levels of lncRNA MSTRG33, ofu_lnc_004834 and ofu_lnc_005214 in diapause vs. no-diapause (from left to right mean ± SEM: 1.0672 ± 0.271, 7.3254 ± 1.5441, 1.0471 ± 0.2332, 0.0007 ± 0.0001, 1.1665 ± 0.3857, 22.1190 ± 5.5802, n = 3). L Relative expression levels of lncRNA MSTRG16991, ofu_lnc_001159 in diapause vs. pre-diapause (from left to right mean ± SEM: 1.0369 ± 0.20120, 0.0374 ± 0.0045, 1.0164 ± 0.1312, 0.0460 ± 0.0142, n = 3)

After the lncRNAs identification and filtering steps, differential expression analysis was conducted to identify differentially expressed lncRNAs between no-diapause, pre-diapause, and diapause stages of ACB. PCA analysis was performed to visualize the variation between samples, and the results showed that samples from the same stage clustered together (Additional file 1: Fig. S7A-C), indicating high reproducibility and consistency between replicates.

Differentially expressed lncRNAs were identified based on the RNA-seq data (Fig. 4F–H). Comparing pre-diapause vs. no-diapause, there were 35 lncRNAs with increased abundance and 7 with decreased abundance. Comparing diapause vs. no-diapause, there were 159 differentially expressed lncRNAs with increased abundance and 141 with decreased abundance. Finally, comparing diapause vs. pre-diapause, there were 198 differentially expressed lncRNAs with increased abundance and 308 with decreased abundance. These results suggest that lncRNAs may play an important role in the regulation of diapause in ACB.

In order to gain a deeper understanding of the biological functions of lncRNAs, we conducted target predictions and GO and KEGG enrichment analyses of these lncRNA targets (Additional file 1: Fig. S7D, E). GO enrichment analysis showed that several processes were enriched in both the diapause vs. pre-diapause and diapause vs. no-diapause comparisons, including cytoplasmic side of early endosome membrane, structural constituent of cuticle, and integral component of mitochondrial inner membrane, among others. These findings suggest that these processes may play important roles in metabolic regulation during diapause. Additionally, the diapause vs. pre-diapause comparison showed enrichment in endopeptidase activity and actin filament, while the diapause vs. no-diapause comparison showed enrichment in UDP-glycosyltransferase activity and chitin-based extracellular matrix, among others. These results suggest that important changes in metabolism and cellular function occur when insects enter or exit diapause. In the pre-diapause vs. no-diapause comparison, enrichment was found in several processes related to the extracellular matrix, such as chitin-based extracellular matrix and basement membrane, as well as in integral component of plasma membrane and mitochondrial matrix. This suggests that changes in extracellular matrix composition and plasma membrane function may be important in the transition from pre-diapause to no-diapause state. Moreover, enrichment was found in endonuclease activity and axon, among others, indicating potential changes in DNA repair and nervous system function during this transition.

Similarly, KEGG enrichment analysis was conducted to further investigate the pathways associated with the different groups (Additional file 1: Fig. S7E). Both comparison between diapause vs. pre-diapause and diapause vs. no-diapause showed enrichment in many metabolic pathways, including purine metabolism, glycine, serine and threonine metabolism, and the biosynthesis of amino acids. This suggests that these pathways play important roles in metabolic regulation during diapause. Additionally, in comparison to the pre-diapause state, the diapause state also showed enrichment in pathways such as arachidonic acid metabolism and starch and sucrose metabolism, indicating that these pathways may undergo significant regulatory changes upon entering diapause. The comparison between pre-diapause and no-diapause showed significant enrichment in fatty acid metabolism pathways, such as biosynthesis of unsaturated fatty acids, fatty acid metabolism, and fatty acid elongation, as well as other metabolic pathways such as galactose metabolism and tryptophan metabolism. This suggests that significant changes in fatty acid metabolism occur upon entering the pre-diapause state, which may have important implications for insects entering diapause.

To further explore the potential role of lncRNAs and their target genes in the pre-diapause process, we constructed an interaction network between lncRNAs, target genes, and pathways (Additional file 1: Fig. S6F). Our analysis revealed that the target genes of Ofur-lnc-001066 and Ofur-lnc-005150, namely LOC114362590, LOC114352303, and LOC114362590, were enriched in three pathways associated with fatty acid metabolism: fatty acid elongation, biosynthesis of unsaturated fatty acids, and fatty acid metabolism. These pathways have been shown to play important roles in many biological processes, including energy metabolism, cell membrane composition, and signaling pathways, and are also important during insect diapause. Fatty acid elongation is a process that converts medium-chain fatty acids into long-chain fatty acids, which are important components of cell membranes and serve as an energy source during insect diapause. Biosynthesis of unsaturated fatty acids is also crucial for maintaining membrane fluidity and function. In addition, fatty acid metabolism plays a key role in regulating energy metabolism and maintaining lipid metabolism balance. We also conducted interaction analysis of differentially enriched target genes and their interacting lncRNAs in the JH pathway in ACB diapause vs. pre-diapause. We found that out of the 10 interacting mRNAs, 5 showed significant differential expression, while 4 lncRNAs also exhibited significant differential expression. Among them, the change in Ofur-lnc-007063 and its host gene LOC114360142 was remarkably significant, with both showing a significant decrease during diapause (adjusted p-value < 0.01, |log2(fold change)|> 1) (Fig. 4I). This suggests the presence of a certain regulatory interaction between lncRNA and JH pathway-related mRNAs. We also selected 8 LncRNAs for differential expression gene validation in random (Fig. 4K–M). The qPCR results were consistent with the RNA-seq sequencing results, confirming the accuracy of our findings.

In this study, dsRNAi constructs targeting lncRNAs corresponding to five diapause-associated genes were created and used in feeding experiments to assess their impact on diapause rates (Fig. 4J). The dsRNAi for DHAMT was used as a positive control. The effectiveness of the dsRNAs was tested by RT-PCR (Additional file 1: Fig. S8A-E). The dsRNAi for Ofur-lnc-001217 and Ofur-lnc-005152 showed no significant effect, while Ofur-lnc-002794 exhibited a trend of increased diapause rates, although not statistically significant. RNA-seq data indicated a significant decrease in the expression of Ofur-lnc-002794 and its target genes during diapause initiation (Additional file 2: Table S1), suggesting their potential role in promoting diapause. The rise in diapause rates following dsRNAi interference with Ofur-lnc-002794 hints that it may no-diapausely inhibit diapause; thus, its disruption might relieve this inhibition, leading to an increased diapause rate. This observation suggests that Ofur-lnc-002794 and its target genes could be significant in diapause regulation. More extensive tests are needed in the future.

In conclusion, the enrichment of lncRNAs target genes in these pathways suggests their potential involvement in regulating critical metabolic processes during the pre-diapause stage. Further investigation of the intricate molecular mechanisms by which lncRNAs regulate these pathways may provide new insights into the physiological regulation of diapause.

Comparation of the DEGs profile of ACB diapause with other types of diapause revealed the molecular characteristics of larval diapause

In order to compare the regulatory pathways involved in different types of insect diapause, we conducted RNA-seq analysis on public data and performed KEGG enrichment analysis on the significantly DEGs. Through a comparative study of diapause in the ACB with other types of diapause, we revealed the molecular characteristics of larval diapause. Our analysis primarily focused on conducting KEGG enrichment analysis across various diapause types, including egg diapause, larval diapause, and adult diapause (Fig. 5A). The results demonstrated that multiple diapause types exhibited enrichment in common metabolic pathways. Among the enriched pathways, the longevity regulating pathway—multiple species, amino sugar and nucleotide sugar metabolism, glutathione metabolism, carbon metabolism, pyruvate metabolism, drug metabolism—cytochrome P450, mTOR signaling pathway, peroxisome, and purine metabolism were consistently enriched in all three diapause types. This indicates the involvement of common regulatory mechanisms and metabolic pathways across diapause stages. Significantly, these enriched pathways in the ACB’s diapause profile reveal the conservation of essential metabolic processes during larval diapause.

Fig. 5figure 5

Comparison of KEGG pathway enrichment analysis among different diapause types. A Dot plot of comparison of KEGG pathway. Each point or triangle represents a KEGG pathway, with the color intensity indicating the degree of enrichment, where redder the color, the more significant the pathway. The size of the points or triangle corresponds to the GeneRatio, showing a positive correlation. Notably enriched pathways are marked with triangular symbols (p-value < 0.05). B Heatmap of specific pathway genes of the ACB. C Heat map of insect hormone biosynthesis related gene

These pathways, such as autophagy, pentose and glucuronate interconversions, pyruvate metabolism, glycerolipid metabolism, sphingolipid metabolism, inositol phosphate metabolism, and galactose metabolism, are significantly enriched during the reproductive diapause of fruit flies, in contrast to ACB larvae. This exclusivity holds significant implications, indicating that these pathways may play distinct roles in nutrient utilization, energy regulation, and cellular adaptation during the reproductive diapause of fruit flies.

We also identified several pathways that were significantly enriched during diapause in ACB larvae. These pathways include insect hormone biosynthesis, DNA replication, alanine, aspartate, and glutamate metabolism, ECM-receptor interaction, Wnt signaling pathway, nucleotide metabolism and FoxO signaling pathway. These pathways exhibited prominent enrichment specifically during diapause in ACB larvae, highlighting their preferential involvement in this developmental stage.

During the diapause of ACB larvae, a pronounced enrichment of the α-linolenic acid metabolism, pyrimidine metabolism, caffeine metabolism, and sulfur metabolism pathways has been observed, potentially holding significant implications. The α-linolenic acid metabolism is linked to the synthesis of ω-3 fatty acids, which play a crucial role in various physiological processes such as cell membrane structure and signaling molecules. During diapause, larvae experience metabolic slowdown and reduced activity, and the alterations in fatty acid metabolism might reflect the demand for adapting to an energy-saving mode during diapause. The pyrimidine metabolism is associated with the synthesis of nucleotides for building DNA and RNA. Considering the metabolic deceleration during diapause, the enrichment of this pathway may be related to the maintenance of fundamental cellular processes, such as DNA repair and maintenance, even under reduced metabolic activity. The enrichment of caffeine metabolism could indicate a response to environmental stress. Caffeine is believed to have a protective role against oxidative stress and participates in regulating various metabolic pathways. The accumulation of this pathway could signify an adaptive response to potential stress factors during diapause.

In these pathways, we also identified specific pathways and genes with significant changes unique to the diapause process of the ACB (Fig. 5B). The JH biosynthesis pathway is crucial in the insect diapause process. We observed a significant increase in the expression of the upstream gene, farnesyl diphosphate phosphatase (FPPP), promoting JH synthesis when compared to the no-diapause state. Simultaneously, we found a significant inhibition of LOC114361407, which encodes NADP+-dependent farnesol dehydrogenase, an enzyme in the corpora allata involved in JH synthesis. Additionally, the expression of genes encoding aldehyde dehydrogenase (NAD+) (LOC114357983) and JH-III synthase (LOC114365517), involved in JH biosynthesis, was altered, leading to the induction of diapause in the ACB (Fig. 5C).

Exploring the influence of HSPs and the proteasome pathway on diapause in the ACB

We found that the HSP70 gene family showed a significant dynamic change during the diapause process in the ACB, and their expression levels were different at different diapause stages, both at the transcriptional and translational levels. At the transcriptional level, they generally showed a pattern of decreasing first and then increasing, indicating that the HSP70 family genes were suppressed and then facilitated the entry into diapause. At the translational level, this change was very significant (p < 0.05), for both hsp68-like and HSP70 genes (Fig. 6A). We also found that the hsp70 (LOC114359959) gene had a markedly increased transcription level and a markedly decreased translation level after entering diapause (Fig. 6B). When treated with the HSP activator, the diapause percentage significantly increased, while treating them with the HSP inhibitor significantly decreased the diapause percentage (Fig. 6C) [32, 33]. The effectiveness of the activator and inhibitors for HSP protein was verified by examining the expression of HSP proteins (Additional file 1: Fig. S8F–H). These results indicate that HSPs plays an important regulatory role in the diapause process of the ACB.

Fig. 6figure 6

Analysis of heat shock protein expression in relation to different diapause types. A Representative Integrative Genomics Viewer (IGV) plots for the transcription and translation profiles of four HSPs. B Comparative box plots of RNA-seq and Ribo-seq expression for the four corresponding HSPs. These expression patterns have been subjected to statistical analysis using the one-way analysis of variance (ANOVA). C Statistical representation of diapause rate under treatment with HSPs inhibitors and activators. D Comparison of diapause percentage of ACB under HSPs inhibitors and activators. * stands for p < 0.05, ** stands for p < 0.01, *** stands for p < 0.001, ns stands for not significant

This study also revealed the involvement of proteasome pathway in diapause of ACB. In cluster 3, there is an initial and substantial decrease in gene expression levels, followed by a sustained, unchanging pattern. Results shown above indicated that cluster 3 enriched a group of genes that are associated with cellular pathways related to proteasomes and ubiquitin-mediated degradation (Fig. 1G). In GSEA, to investigate RNA-seq data from different diapause stages, our analysis revealed the restriction and suppression of pathways related to proteasomes during the transition from the no-diapause state to diapause (Fig. 7A, B). To validate this process, treatment of ACB using MG132, an inhibitor of proteasome activity, led to an increase in percentage of diapause at a concentration of 100

Comments (0)

No login
gif