Time-Series Transcriptomic Analysis Reveals the Molecular Profiles of Diapause Termination Induced by Long Photoperiods and High Temperature in Chilo suppressalis

Survival and adaptation to seasonal changes are challenging for insects. Many temperate insects such as the rice stem borer (Chilo suppressalis) overcome the adverse situation by entering diapause, wherein development changes dynamically occur and metabolic activity is suppressed. The photoperiod and temperature act as major environmental stimuli of diapause. However, the physiological and molecular mechanisms that interpret the ecologically relevant environmental cues in ontogenetic development during diapause termination are poorly understood. Here, we used genome-wide high-throughput RNA-sequencing to examine the patterns of gene expression during diapause termination in C. suppressalis. Major shifts in biological processes and pathways including metabolism, environmental information transmission, and endocrine signalling were observed across diapause termination based on over-representation analysis, short time-series expression miner, and gene set enrichment analysis. Many new pathways were identified in diapause termination including circadian rhythm, MAPK signalling, Wnt signalling, and Ras signalling, together with previously reported pathways including ecdysteroid, juvenile hormone, and insulin/insulin-like signalling. Our results show that convergent biological processes and molecular pathways of diapause termination were shared across different insect species and provided a comprehensive roadmap to better understand diapause termination in C. suppressalis.


Introduction
The capacity to enter diapause is pervasive in insects which enables them to survive adverse seasonal environments such as harsh winters and synchronise their life cycle with favourable conditions for development and reproduction [1,2]. Diapause is not a cessation of development; instead, it is a dynamic alternate developmental process involving physiologically distinct phases: pre-diapause, diapause, and post-diapause [3]. Furthermore, diapause is a succession of three subphases: initiation, maintenance, and termination [3,4]. Insect diapause can be obligatory, with no external environmental stimuli and occur regardless of environmental conditions, or facultative, wherein specific environmental cues trigger a response and induce diapause [3]. Insects have diverse diapause timings and may occur in the egg, larva or nymph, pupa, or adult stage [4,5]. The photoperiod and temperature are the primary environmental signals for facultative diapause initiation [6]. Especially the photoperiod is the major signal through which most insects anticipate seasonal changes [7].
Environmental signals can also stimulate insects to terminate diapause. This occurs when environmental conditions are favourable with the stimuli passing a critical threshold, allowing insects to continue ontogenetic development [8]. Termination of diapause under field conditions (horotelic termination) is a relatively slow process [9]. In contrast, termination of diapause under laboratory conditions (tachytelic termination) can be rapid [9,10]. Insects can perceive various environmental signals as stimuli for diapause termination. For example, low temperature is a necessary signal to terminate diapause for some insects such as Chymomyza costata [4], whereas moisture is more important for other insects [11]. The photoperiod is the most reliable and steady environmental signal amongst signals that seasonally change. Therefore, many insects have evolved to utilise the photoperiod for not only inducing of diapause but also for terminating it [3]. Although not as steady as the photoperiod and temperature, as mentioned earlier, it also plays an important role in the termination of diapause [4,12,13].
The rice stem borer Chilo suppressalis is one of the most serious pests of rice in Asia [14] and causes significant annual damage to crop yield [15]. C. suppressalis larvae are capable of undergoing facultative diapause when the day length is shorter than the critical duration in autumn [16], even though the average temperature is still above 27 • C in the field [16]. Diapausing C. suppressalis larvae are able to endure cold and lack of food, which enables them to survive upcoming harsh winters [17]. Laboratory studies showed that C. suppressalis can be induced to enter diapause by a short photoperiod, and night length is a key component for initiating diapause [16,18,19]. In contrast, diapause can be terminated when field-collected C. suppressalis larvae (from the winter period) are exposed to long photoperiods (light:dark, 16:8) and high temperatures (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) • C) in the laboratory [16].
Many studies have focused on revealing the physiological and molecular mechanisms of diapause under field or laboratory conditions. Transcriptomic analysis of field-collected pre-diapausing C. suppressalis larvae revealed that pathways involving neuroactive ligandreceptor interaction, lysosomes, and glycerolipid metabolism were associated with prediapause [20]. Comparative transcriptomic analysis of 11 insect species revealed that a set of core genes involved in circadian rhythmicity, insulin signalling, and Wnt signalling were differentially regulated during diapause [21]. However, our current understanding of the molecular mechanisms of diapause termination is unclear.
In this study, diapausing C. suppressalis larvae and diapause-terminating larvae induced by long photoperiods (light:dark 16:8) and high temperature (25 • C) were collected at three different time points and their transcriptomes sequenced. Examination of the patterns of gene expression in C. suppressalis revealed the biological processes and molecular pathways involved in diapause termination in a panoramic manner based on over-representation analysis, short time-series expression miner, and gene set enrichment analysis.
A total of 85.75 Gb of clean bases (number of clean reads multiply read length, saved in Gb unit) containing 580.78 million clean reads (reads that passed quality control) was generated from these four groups with 12 samples in all, and each sample contained ≥6.42 Gb of clean bases (43.41 million clean reads) ( Table 1). Between 39.93 and 51.21 million clean reads were mapped to the reference genome [22] for each sample, with total map read rates from 84.18 to 93.50%. Approximately 4.22 to 7.14% of total reads were mapped to multiple locations (Table 1). Sequence assembly and functional annotation revealed 25,913 genes, including 14,799 known genes and 11,114 newly detected genes. Together with the transcriptome coverage calculated (Table 1), all results showed that the sequencing had good depth and coverage [23].

Detection of Differentially Expressed Genes (DEGs)
DEGs were identified by comparing the genes in a total of six possible pairs of diapausing C. suppressalis sample groups (WSG1) and three other treated sample groups (WSG2, WSG3, WSG4). A total of 3734 genes were differentially expressed in at least one of the six comparisons ( Figure 1A). Among them, 2425 DEGs were identified in WSG1_vs_WSG2 (1385 upregulated genes and 1040 downregulated genes). This was far more than the number identified in WSG2_vs_WSG3 (70 genes) and WSG3_vs_WSG4 (54 genes). Meanwhile, 2037 DEGs (1199 upregulated genes and 838 downregulated genes) were identified in WSG1_vs_WSG4, which is close to WSG1_vs_WSG3 (2056 DEGs; 1096 upregulated genes and 960 downregulated genes). There were 1056 shared or common DEGs in the following comparisons: WSG1_vs_WSG2, WSG1_vs_WSG3, and WSG1_vs_WSG4 ( Figure 1B), which represented the DEGs between the state of diapause and terminating diapause. Shared DEGs were classified on the basis of Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) annotations to reveal the biological processes and pathways associated with diapause termination. GO classification was performed for three major categories: biological process, molecular function, and cellular component. Most DEGs were assigned to metabolic and cellular process, followed by biological regulation in biological processes. Meanwhile, most DEGs were categorised into catalytic activity, binding, and transporter activity in the molecular function group. Finally, most DEGs were classified to membrane-and cell parts in terms of cellular components ( Figure 2, Supplementary Table S1). KEGG orthology classification was applied to the shared DEGs in all top categories of the KEGG pathway maps. Lysosomes had the most DEGs in the cellular process category, while the most predominant DEGs were found in the cAMP signalling pathway, neuroactive ligand-receptor interaction, and AMPK signalling pathway in the environmental information-processing category. Protein pro- Shared DEGs were classified on the basis of Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) annotations to reveal the biological processes and pathways associated with diapause termination. GO classification was performed for three major categories: biological process, molecular function, and cellular component. Most DEGs were assigned to metabolic and cellular process, followed by biological regulation in biological processes. Meanwhile, most DEGs were categorised into catalytic activity, binding, and transporter activity in the molecular function group. Finally, most DEGs were classified to membrane-and cell parts in terms of cellular components ( Figure 2, Supplementary Table S1). KEGG orthology classification was applied to the shared DEGs in all top categories of the KEGG pathway maps. Lysosomes had the most DEGs in the cellular process category, while the most predominant DEGs were found in the cAMP signalling pathway, neuroactive ligand-receptor interaction, and AMPK signalling pathway in the environmental information-processing category. Protein processing in the endoplasmic reticulum contained the most DEGs for the genetic information processing category, while drug metabolism and pancreatic secretion DEGS were mostly found in metabolism and organismal systems, respectively (

Short Time-Series Expression Miner
A short time-series expression miner (STEM) was applied to all DEGs from all six possible comparisons of WSG1, WSG2, WSG3, and WSG4 to detect temporal differences in gene expression profiles associated with diapause termination.
Twenty profiles were clustered according to gene expression. Eight of these profiles showed significant expressional patterns in the process of diapause termination (p < 0.05) ( Figure 4). Profiles 0, 2, and 6 contained clusters of genes exhibiting downregulated expression trends, whereas the genes of all other profiles showed upregulated expression trends. Profile 0 (737 genes) and profile 19 (832 genes) clustered more DEGs than the other profiles. Subsequently, GO enrichment analysis was performed on each of these clustered profiles to identify overrepresented biological pathways or processes. Profiles 0, 2, 6, 12, 13, 17, 18, and 19 had significantly enriched GO terms (q value < 0.05) ( Table 2). Profile 0 showed significantly enriched genes related to biosynthetic processes (including small molecules, fatty acids, carboxylic acids, and monocarboxylic acids) and catabolic processes (including organic substances and organonitrogen compounds) (

Short Time-Series Expression Miner
A short time-series expression miner (STEM) was applied to all DEGs from all six possible comparisons of WSG1, WSG2, WSG3, and WSG4 to detect temporal differences in gene expression profiles associated with diapause termination.
Twenty profiles were clustered according to gene expression. Eight of these profiles showed significant expressional patterns in the process of diapause termination (p < 0.05) (Figure 4). Profiles 0, 2, and 6 contained clusters of genes exhibiting downregulated expression trends, whereas the genes of all other profiles showed upregulated expression trends. Profile 0 (737 genes) and profile 19 (832 genes) clustered more DEGs than the other profiles. Subsequently, GO enrichment analysis was performed on each of these clustered profiles to identify overrepresented biological pathways or processes. Profiles 0, 2, 6, 12, 13, 17, 18, and 19 had significantly enriched GO terms (q value < 0.05) ( Table 2). Profile 0 showed significantly enriched genes related to biosynthetic processes (including small molecules, fatty acids, carboxylic acids, and monocarboxylic acids) and catabolic processes (including organic substances and organonitrogen compounds) (Table 2). Profile 19 only showed significantly enriched GO terms in small molecule catabolic processes and extracellular regions.    KEGG enrichment analysis was also implemented to each profile. Profiles 0, 13, 17, 18, and 19 were identified as significant KEGG pathways (q-value < 0.05) ( Table 3). Profile 0 showed significantly enriched genes involved in the PPAR signalling pathway, AMPK signalling pathway, Rap1 signalling pathway, GTP-binding proteins, and proteasome. Meanwhile, steroid biosynthesis, insect hormone biosynthesis, pancreatic secretion, and PPAR signalling pathways were enriched in profile 19 ( Table 3). The pathways of insulin secretion and neuroactive ligand-receptor interactions were enriched in the other profiles (Table 3).

Changes during the Onset of Diapause Termination
Gene set enrichment analysis (GSEA) using GO gene sets and KEGG pathway gene sets was applied to WSG1 vs. WSG2 to reveal the functional interaction within all expressed genes at the onset of diapause termination. There were 148 negatively enriched gene sets (normalised enrichment score (NES) with a minus sign) and 59 positively enriched gene sets (NES with a positive sign) out of the total 431 detected gene sets from the KEGG pathway (Supplementary Table S3). The top 10 negatively enriched KEGG orthology (KO) terms of each category (except Brite hierarchies and human disease) ordered by the absolute value of NES are illustrated in Figure 5. Most of these pathways (such as the Rap1 signalling pathways and MAPK signalling pathway) were also found in KEGG annotations of shared DEGs, which supported and validated the results.
GO gene set analysis showed that 137 out of the 2786 detected gene sets were significantly negatively enriched and 92 positively enriched (Supplementary Table S4). The top 10 negatively enriched GO terms ordered by the absolute value of NES are shown in Figure 5F-H and included protein transport, intracellular signal transduction, and GTP binding. The negative sign of NES demonstrates that the gene sets were activated in state of WSG1.  Table S4). The top 10 negatively enriched GO terms ordered by the absolute value of NES are shown in Figure 5F-H and included protein transport, intracellular signal transduction, and GTP binding. The negative sign of NES demonstrates that the gene sets were activated in state of WSG1.

Changes Occurring in Late-Stage Diapause Termination
The first pupa appeared after 20 d of exposure to long-day photoperiods. Thus, the time point WSG4 (12th day) is likely to be in the late stage of diapause termination. Therefore, the dominant biological processes and physiological characteristics of late-stage diapause termination can be revealed by GSEA analysis applied to WSG1 and WSG4.
Seventy-four gene sets out of the total detected KEGG gene sets (431) (Supplementary Table S5) and 11 gene sets out of the total detected GO gene sets (2876) (Supplementary Table S6) were significantly negatively enriched, while these positively enriched KEGG and GO gene sets were 84 and 80, respectively (Supplementary Tables S5 and S6). The top 10 negatively enriched terms of each category for KEGG (except Brite hierarchies and human disease) ordered by the absolute value of NES are illustrated in Figure 6. There were much fewer negatively enriched KEGG gene sets or GO gene sets than those of WSG1 vs. WSG2 (148 and 137, respectively), and no gene set was negatively enriched under the GO category of cellular component. Cellular processes, autophagy, and cell cycle pathways were predominant ( Figure 6A). The MAPK signalling pathway, ribosome biogenesis, terpenoid backbone biosynthesis, and circadian rhythm were ranked ahead in other KEGG categories ( Figure 6B-E).

Discussion
C. suppressalis larvae exhibit facultative diapause, which is mainly affected by the photoperiod [16]. Larvae enter diapause in autumn (short photoperiods) and terminate it the following spring (long photoperiods co-occur with high temperature) when undisturbed in the field. A state of diapause allows C. suppressalis to endure and survive harsh winter conditions and prolonged periods without food to continue self-perpetuating populations when conditions are favourable. This study characterised the gene expression patterns across the developmental time course of diapause termination in C. suppressalis. Field-collected larvae were exposed to long photoperiods and high temperatures in the laboratory to terminate diapause, which contrasts with previous studies involving extreme endocrine or pharmacological manipulation. Previous studies mainly focused on comparisons of non-diapausing and diapausing states [20,[24][25][26], or clarifying the threshold of environmental stimuli to induce or terminate diapause [16,19,27]. The results of our study provide a roadmap to understand the physiological and molecular mechanisms associated with insect diapause termination, especially in Lepidoptera.
When taking the order of events in diapause termination into account, it is reported that the increasing metabolic rates signal the onset of diapause termination and resumption of development. Studies in Rhagoletis pomonella demonstrated that the timing of thousands of DEGs co-occurred with increased metabolism during diapause termination [28], and an increased metabolic rate is a reliable indicator of this event [29]. In other Error bars indicate the standard error among biological replicates. NS is not significant, * p-adjusted < 0.05, ** p-adjusted < 0.01 using Student's t-test and multiple testing corrected with Benjamini-Hochberg method. (B) Gene expression fold change of selected genes using WSG1 as a control. Significance was determined by DEseq2, * p-adjusted < 0.05, ** p-adjusted < 0.01.

Discussion
C. suppressalis larvae exhibit facultative diapause, which is mainly affected by the photoperiod [16]. Larvae enter diapause in autumn (short photoperiods) and terminate it the following spring (long photoperiods co-occur with high temperature) when undisturbed in the field. A state of diapause allows C. suppressalis to endure and survive harsh winter conditions and prolonged periods without food to continue self-perpetuating populations when conditions are favourable. This study characterised the gene expression patterns across the developmental time course of diapause termination in C. suppressalis. Field-collected larvae were exposed to long photoperiods and high temperatures in the laboratory to terminate diapause, which contrasts with previous studies involving extreme endocrine or pharmacological manipulation. Previous studies mainly focused on comparisons of non-diapausing and diapausing states [20,[24][25][26], or clarifying the threshold of environmental stimuli to induce or terminate diapause [16,19,27]. The results of our study provide a roadmap to understand the physiological and molecular mechanisms associated with insect diapause termination, especially in Lepidoptera.
When taking the order of events in diapause termination into account, it is reported that the increasing metabolic rates signal the onset of diapause termination and resumption of development. Studies in Rhagoletis pomonella demonstrated that the timing of thousands of DEGs co-occurred with increased metabolism during diapause termination [28], and an increased metabolic rate is a reliable indicator of this event [29]. In other words, differential expression of thousands of genes could be an indicator of the onset of diapause termination. Our study showed that the most differentially expressed genes were found in the comparison of WSG1 vs. WSG2 ( Figure 1A). This suggests that the metabolic rates in C. suppressalis larvae dramatically increased after exposure to four long and warm days. Changes in the expression of WSG2_vs_WSG3 and WSG3_vs_WSG4 were small compared to WSG1_vs_WSG2 (Figure 1A), indicating that no dramatic change occurred between day 5-12 of the long photoperiod, high temperature treatment. This pattern suggests that C. suppressalis terminated diapause within four days of long light and high temperature at the transcriptional level. However, the precise timing of diapause termination in C. suppressalis requires further study.
A total of 1056 shared DEGs were identified when WSG2, WSG3, and WSG4 were compared to WSG1 ( Figure 1B). These DEGs reflected the difference in gene expression between diapausing and diapause-terminating C. suppressalis. KEGG classification of shared DEGs showed that multiple pathways play essential roles in diapause termination ( Figure 3). Many differentially expressed KEGG pathways found in C. suppressalis diapause termination were also found in R. pomonella (Diptera family member that diapauses in pupa), such as mTOR signalling, purine metabolism, and steroid biosynthesis [28]. Metabolites such as glucose, pyruvate, succinate, and malate increased during diapause termination in Helicoverpa armigera triggered by injecting 20E [30]. Interestingly, pathways such as carbohydrate metabolism, glycan biosynthesis and metabolism, and pyruvate metabolism were classified as shared DEGs (Figure 3). Enzyme activities associated with glycogen metabolism (such as phosphorylase and glycogen synthetase) were significantly altered during diapause in C. suppressalis [31], confirming the credibility of our experimental results. Furthermore, the similarity in diapause termination metabolic processes among different insect species suggests that they may share analogous mechanisms.
GSEA can provide insights into biological mechanisms in gene sets (i.e., KEGG pathways and GO categories) level between two states of interest, particularly that are subtle at the level of individual genes [32]. Accordingly, the onset and late stage of diapause termination of C. suppressalis were analysed with GSEA in this study. Since the WSG1 was set as control, negatively enriched pathways and biological processes are these that activated in the state of WSG1 and positively enriched ones are these that activated in the state of WSG2 or WSG4. In other words, negatively enriched pathways and biological processes are functioning to maintain insects in state of diapause. Therefore, they are important to explain the molecular changes of diapause termination.
Time-series analysis confirmed that major gene expression changes occurred in the first period of the time course involving diapause termination, especially profiles 0 and 19, which contained nearly half of the genes analysed ( Figure 4). Many GO terms related to biosynthetic and catabolic processes were enriched in these two profiles ( Table 2), suggesting consistency in the increasing metabolic rate during the initiation of diapause termination [28].
Extensive research has shown that changes in environmental cues initiate a cascade of endocrine events that control diapause [5,33]. Steroid biosynthesis and insect hormone biosynthesis pathways enriched in profile 19 showed consistency in the endocrine system in response to changes in photoperiod length (Table 3). Moreover, these endocrine events occurred in the early stages of diapause termination (Figure 4).
The ecdysteroid-, juvenile hormone (JH)-, and insulin/insulin-like signalling pathways regulate responses to environmental changes and are involved in diapause in many insects [5,34]. In pre-adult diapause, ecdysteroids (such as 20-hydroxyecdysone in insects) are generally regarded as diapause terminators [5,34]. Meanwhile, JH induces and maintains diapause in C. suppressalis [35]. In Culex pipiens, JH triggers diapause termination, and JH production is regulated by the insulin/insulin-like signalling pathway and the downstream FoxO gene [36]. Insulin-like peptides and its downstream FoxO signalling pathway increase trehalose anabolism and trigger diapause termination in Antheraea pernyi [37]. Hydrogen peroxide-induced diapause termination in Artemia cyst embryos showed that DEGs are involved in the Hippo-and MAPK signalling pathways [38]. Extracellular signalregulated kinase (ERK) (a member of the mitogen-activated protein kinase (MAPK) family) is activated during the egg diapause-terminating process in Locusta migratoria [39]. ERK/MAPK regulates ecdysteroid and sorbitol metabolism in Bombyx mori through the transcription of key enzymes (sorbitol dehydrogenase-2 and EPPase) for embryonic diapause termination [40,41]. These pathways were also identified in this study during C. suppressalis diapause termination using STEM-or GSEA analysis. These consistencies validate the reliability of our results and demonstrate the convergence of mechanisms of diapause termination between different insects.
Currently, the links between environmental information such as the photoperiod and endocrine signals are poorly understood [42]. However, numerous studies have indicated that circadian clock genes and their regulators contribute to diapause regulation. A period gene knockout line of B. mori shows that the GABAeric pathway mediates the photoperiod information to the diapause hormone pathway [43]. Circadian transcription factors vrille and Par domain protein 1 may regulate signalling pathways underlying arrested eggs follicle development and fat accumulation in diapausing Culex pipiens females [44]. Circadian genes period and timeless promote lipid accumulation during diapause preparation in the diapause destined female Colaphellus bowringi [45]. N-acetyltransferase is considered a critical conjunction of photoperiodism between the circadian system and endocrine axis in Antheraea pernyi [46]. Our GSEA analysis of WSG1 and WSG4 pairs showed that the circadian rhythm pathway (ko04710) was the top enriched gene set (by NES) in the organismal system category (Figure 6), suggesting that this pathway also plays an important role in diapause termination of C. suppressalis. It is notable that the circadian rhythm pathway was not enriched in GSEA analysis in the comparison of WSG1_vs_WSG2. This likely indicates that only a few core genes participated in the early process of photoperiod information transmissions. This finding should be validated in future studies.
The identification of candidate genes and signalling pathways for diapause termination in C. suppressalis provides a foundation for understanding how environmental cues are transformed into endocrine events that regulate the diapause process. The key genes and pathways may be utilised to disturb the diapause process of C. suppressalis as a green strategy to control this pest insect. For instance, if C. suppressalis terminates diapause in an untimely manner, it cannot survive harsh conditions. Generally, exploring the mechanisms of diapause regulation will facilitate comparative research on regulatory factors conserved across insects [21]. Our results show that convergent biological processes and molecular pathways are shared across different insect species, although the specific genes regulating diapause termination may not be the same.

Insects and Collection of Samples
Diapausing C. suppressalis larvae were collected from rice straw left at the paddy side from Huai'an, Jiangsu, in January 2021. A random sample of diapausing larvae was immediately collected and designated WSG1. The remaining diapausing larvae were reared in the laboratory at 25 ± 1 • C under long photoperiod light:dark cycles (16:8 h) and a relative humidity of 70 ± 5% to terminate diapause. One sample group was collected every four days (at 9:00 a.m.) and a total of three groups (WSG2-WSG4) were collected. The remaining larvae were reared until pupation stage. The interval between the last sampling day and the day where the pupa was observed was noted. Each sample group contained three replicates, and each replicate consisted of four C. suppressalis larvae. All collected samples were frozen in liquid nitrogen and stored at −80 • C.

Differential Expression Gene Analysis
DEGs were identified between sample groups by calculating the expression level of each transcript according to the transcript per million (TPM) mapped reads method. RNAsequencing by expectation-maximization (RSEM, v. 1.3.3, http://deweylab.biostat.wisc. edu/rsem/ (accessed on 29 October 2021)) was used to quantify gene abundance. The R package DESeq2 (differential expression analysis for sequence count data 2, v. 1.24.0) was used to analyse the DEGs. Genes with an absolute value fold-change ≥2 and p-adjusted value < 0.05 (using the Benjamini-Hochberg method for multiple testing correction) were considered differently expressed.

GO and KEGG Enrichment Analysis
GO and KEGG enrichment analyses were performed based on over-representation analysis (ORA) with the R package ClusterProfilter (v. 4.2.1) [47]. GO terms and KEGG pathways with p-adjusted values < 0.05 were used for subsequent analyses. Multiple testing was performed with the Benjamini-Hochberg method.

Short Time-Series Expression Miner
Short time-series expression miner (STEM) software (v. 1.3.11) was employed to cluster all DEGs from pairwise comparisons of WSG1-WSG4 to detect temporal gene expression patterns. Gene expression values were transformed into log-normalised data which served as input data. STEM analysis was performed using the STEM clustering method with default settings, except that the maximum unit change in the model profiles between time points was changed to 4. Each expression profile represented a gene set associated with the time-course of diapause termination in C. suppressalis. Gene sets with a p value < 0.05 were considered significant.

Gene Set Enrichment Analysis
Gene set enrichment analysis using GO gene sets and KEGG pathway gene sets were applied to WSG1 vs. WSG2 and WSG1 vs. WSG4, respectively, wherein WSG1 served as a control when calculating the log2 fold change of gene expression. This was performed to reveal the functional interaction within all expressed genes in different time courses involved in diapause termination. The WSG2 time point reflects transcriptional changes in the initiation of diapause termination, while the WSG4 time point shows transcriptional changes in the posterior process of diapause termination.
Genes in paired comparisons were ranked according to log2 transformed expression fold change, and the ranked gene list was analysed using GSEA package (v. 4.2.1, https: //www.gsea-msigdb.org/gsea/index.jsp (accessed on 6 January 2022)). The analysis was performed using the GSEAPreranked tool with default settings except that the enrichment statistic was set to classic and min size was set to 1. KEGG and GO annotations of the C. suppressalis transcriptome were used as enrichment gene set data. Enriched gene sets with a nominal p-value < 0.05, FDR q-value < 0.25, and absolute value of normalised enrichment score (NES) > 1 were considered significant.

qPCR Verification
Six randomly selected DEGs were evaluated by qPCR to verify their expression levels in RNA sequencing. Specific primers (Supplementary Table S7) were designed using Primer3 (https://bioinfo.ut.ee/primer3-0.4.0/ (accessed on 21 February 2022)). The qPCR template was prepared from the RNA used for transcriptomic sequencing. Total RNA (1000 ng) was reverse-transcribed to cDNA using PrimeScript TM RT master mix (Perfect Real Time) (TaKaRa, Kusatsu, Japan). Each qPCR reaction was prepared as follows: 10 µL TB Green premix EX Taq II (Tli RNaseH plus) (TaKaRa), 2 µL cDNA template, 7.2 µL ddH 2 O, 0.4 µL forward primer (10 µM), and 0.4 µL reverse primer (10 µM). Each gene was analysed using three biological replicates and each biological replicate had at least three technical replicates. The reactions were carried out using a program that consisted of an initial heat activation step of 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s and 60 • C for 30 s in a LightCycler 480 II (Roche, Basel, Switzerland). Relative quantification was performed using the 2 −∆∆Ct method [48] with histone 3 as a reference gene [49]. The relative expression levels of genes were compared using Student's t-test instead of the multiple comparison test across all four time-points since all of these DEGs were obtained by pairwise comparison of the control (WSG1) and treatment (WSG2, WSG3, WSG4, respectively). The t-test p-values were corrected with the Benjamini-Hochberg method for multiple testing.

Conclusions
Our results provide a panoramic view of the molecular mechanisms of diapause termination in C. suppressalis in a time-series manner. The identified genes and molecular pathways associated with diapause termination validate the convergent molecular mechanisms of different insects. The new candidate genes and molecular pathways identified in this study should help integrate the understanding of the mechanisms of diapause termination.