Integrative Proteomic and Phosphoproteomic Analyses Revealed Complex Mechanisms Underlying Reproductive Diapause in Bombus terrestris Queens

Simple Summary Bombus terrestris is one of the most ideal pollinators in the world and brings great economic benefits through the pollination of fruits and vegetables. The huge demand for pollinators, combined with habitat loss, climate change, and pesticide use over the past decades, has resulted in a dramatic decline in wild bumblebee populations and distribution. Hence, the artificial breeding of bumblebees is of great significance and has a broad market prospect. Diapause is an important process in artificial breeding. Although some studies have been carried out on the diapause-related genes of bumblebees, the precise mechanisms affecting the phosphorylation level remain unclear. Here, we performed a new comparison of three diapause-stage expression profiles using isobaric tandem mass tag (TMT)-labeled proteomics and phosphoproteomics. The results provided abundant resources and contributed to a better understanding of the mechanisms underlying the regulation of reproductive diapause in eusocial insects. Abstract Reproductive diapause is an overwintering strategy for Bombus terrestris, which is an important pollinator for agricultural production. However, the precise mechanisms underlying reproductive diapause in bumblebees remain largely unclear. Here, a combination analysis of proteomics and phosphoproteomics was used to reveal the mechanisms that occur during and after diapause in three different phases: diapause (D), postdiapause (PD), and founder postdiapause (FPD). In total, 4655 proteins and 10,600 phosphorylation sites of 3339 proteins were identified. Diapause termination and reactivation from D to the PD stage were characterized by the upregulation of proteins associated with ribosome assembly and biogenesis, transcription, and translation regulation in combination with the upregulation of phosphoproteins related to neural signal transmission, hormone biosynthesis and secretion, and energy-related metabolism. Moreover, the reproductive program was fully activated from PD to the FPD stage, as indicated by the upregulation of proteins related to fat digestion and absorption, the biosynthesis of unsaturated fatty acids, fatty acid elongation, protein processing in the endoplasmic reticulum, and the upregulation of energy-related metabolism at the phosphoproteome level. We also predicted a kinase–substrate interaction network and constructed protein–protein networks of proteomic and phosphoproteomic data. These results will help to elucidate the mechanisms underlying the regulation of diapause in B. terrestris for year-round mass breeding.


Introduction
Many insect species enter a period of developmental or reproductive stasis, called diapause, as an adaptive strategy to enhance survival in response to unfavorable conditions and seasonal changes [1]. Diapause is a complicated and dynamic developmental process characterized by decreased metabolic activity and increased stress resistance [2]. Various neuroendocrine, molecular, cellular, enzymatic, metabolic, endocrine, and behavioral changes occur during diapause [1][2][3]. Diapause usually occurs at species-specific developmental stages (i.e., embryo, larva/nymph, pupa, and adult) [1,[4][5][6]. Reproductive diapause often occurs in overwintering adults upon the arrest of vitellogenesis and ovulation [7]. Diapause is classified into three major successive stages: prediapause, diapause, and postdiapause. Prediapause is defined as the induction and preparation for diapause, while diapause includes the initiation, maintenance, and termination of diapause, and postdiapause is characterized by quiescence and rapid development [1]. During prediapause, insects perceive changes in the environment and regulate multiple physiological processes. Resources are accumulated, and maturation is delayed. During diapause, insects are in a state of "suspended animation", characterized by a low metabolic rate and enhanced stress resistance. In the postdiapause stage, the organisms return to an active state to complete the life cycle. The biological and cellular processes involved in these major physiological transitions are complex and interesting, especially the transition from diapause to the postdiapause stage, which requires the restarting of different groups of biological processes and needs to be further studied [8,9]. Here, the proteins and phosphoproteins involved in the transitions between the successive stages of diapause in bumblebee (Bombus terrestris) queens were investigated.
B. terrestris is a pollinator that plays important roles in maintaining natural ecosystems and facilitating agricultural production [10]. Over the past 20 years, facility agriculture has rapidly developed along with global agricultural initiatives. Some bumblebee species have been successfully domesticated and raised on a large scale to meet the pollination needs of facility agriculture. As one of the most widely used commercial pollinators globally, B. terrestris is economically beneficial to the facility production of fruits and vegetables [11]. Low temperatures induce reproductive diapause in adult B. terrestris [12]. In the wild, B. terrestris senses temperature changes to terminate diapause and start to development when spring comes. Factory production also mimics natural conditions. Diapause is an important factor restricting colony formation in the artificial mass rearing of B. terrestris. The elucidation of the regulatory mechanisms underlying reproductive diapause is a long-standing goal in organismal biology and will facilitate sustainable year-round mass breeding of key pollinator species.
Multiple regulatory mechanisms have been implicated in diapause. Protein phosphorylation is a key regulatory event that is reported to play very important roles in insect physiology, reproduction, and development. In the fruit fly (Drosophila melanogaster), calcineurin, a ubiquitous serine/threonine protein phosphatase, broadly influences protein phosphorylation during oocyte maturation and egg activation [7]. Phosphorylation regulatory mechanisms are also involved in the acquisition of thermal tolerance [13]. To date, 27 potential circadian kinases and 789 phosphorylation sites with circadian oscillations have been identified in D. melanogaster [14]. In the silkworm (Bombyx mori), several key protein-synthesis-related pathways are regulated by phosphorylation, thereby influencing silk production [10]. In the honeybee (Apis mellifera), phosphoproteomic analysis has been widely applied to assess the transition from nurse to forager bee, the embryo-larva transition, the production of royal jelly and venom, and the development of the salivary glands and hypopharyngeal glands [15][16][17][18][19][20]. Previous studies have demonstrated that phosphorylation modifications also participate in the regulation of insect diapause [21][22][23][24][25]. A prior report showed significant differences in 23 phosphoproteins in the brain of the cotton bollworm (Helicoverpa armigera) in the nondiapause-and diapause-destined pupae [26]. In addition, 27 proteins and phosphoproteins involved in cell proliferation, adult Insects 2022, 13, 862 3 of 23 development, and aging were identified in the brain of the flesh fly (Sarcophaga crassipalpis) during the initiation of pupal diapause [27].
Proteomic and phosphoproteomic analyses are regarded as powerful approaches for large-scale investigations of gene expression patterns at the protein and post-translational levels. Many proteomic studies have been reported on insect development, immune response, metamorphosis, and diapause [28][29][30][31][32]. However, few studies have looked into the post-translational phosphorylation involved in insect diapause, and even fewer have gone beyond the limitations of two-dimensional electrophoresis (2-DE) coupled with mass spectrometry (MS) in order to identify proteins; this has led to a lack of insight into the regulatory pathway involved in diapause [26,27].
Recently, a high-sensitivity and high-resolution proteome technology, known as the tandem mass tags (TMT) labeling strategy paired with MS, has been developed, allowing for the identification of thousands of proteins with increased coverage and accuracy [33]. Here, a TMT-label-based quantitative analysis was applied to identify differentially expressed proteins (DEPs) and phosphoproteins (DEPPs) of B. terrestris queens in the diapause (D), postdiapause (PD), and founder postdiapause (FPD) stages.

Bumblebee Sampling
The B. terrestris queens used in this research were purchased from Shandong Lubao Technology Co. Ltd. (Jining, China), an organization dedicated to the commercial production of bumblebees. The sampling was performed as described in a previous study (Chen et al., 2021). In brief, the sampling of B. terrestris queens was conducted at three different developmental stages: D-mated queens were maintained at 2 • C in 24 h dark and 50-60% humidity for 10 weeks; PD-diapaused queens were treated with CO 2 for 1 min, which was repeated once after 24 h at room temperature, and then fed a diet of 50% sucrose for 2 days at 28 • C in 24 h dark and 50-60% humidity; FPD-the postdiapause queens were collected until they began to lay eggs at 28 • C in 24 h dark and 50-60% humidity. All samples were frozen in liquid nitrogen and stored at −80 • C until use. Each group had three biological replicates, and each replicate had three queen bees.

Protein Extraction and Digestion
A schematic overview of the experiment is depicted in Figure 1. One gram of each sample was fully ground into powder in a mortar, which was precooled with liquid nitrogen, then transferred to centrifuge tubes, mixed with four volumes of phenol extraction buffer (10 mM dithiothreitol, 1% protease inhibitor, and 1% phosphatase inhibitor), and sonicated. Following the addition of an equal volume of Tris equilibrium pheno, the sample was centrifuged at 5500× g at 4 • C for 10 min. Afterward, the supernatant was collected, mixed with five volumes of 0.1 M ammonium acetate/methanol, and allowed to precipitate overnight. The precipitated proteins were washed with methanol and acetone then redissolved in 8 M urea. The bicinchoninic acid assay was used to measure the protein concentration.
Each protein sample (300 µg) was subjected to enzymatic hydrolysis. In brief, 20% trichloroacetic acid was slowly added to equal amounts of protein and urea lysate then precipitated at 4 • C for 2 h. The supernatant was discarded after centrifugation at 4500× g for 5 min, and the precipitation was washed twice with acetone that had been precooled, air-dried, and then mixed with 200 mM triethylammonium bicarbonate buffer and ultrasonicated. Then, trypsin was added at a ratio of 1:50 (protease/protein, m/m), and the mixture was enzymolyzed overnight. The next day, 5 mM dithiothreitol was added, and the mixture was reduced at 56 • C for 30 min, followed by the addition of 11 mM iodoacetamide and incubation at room temperature for 15 min in the dark. Each protein sample (300 μg) was subjected to enzymatic hydrolysis. In brief, trichloroacetic acid was slowly added to equal amounts of protein and urea lysate precipitated at 4 °C for 2 h. The supernatant was discarded after centrifugation at 450 for 5 min, and the precipitation was washed twice with acetone that had been precoo air-dried, and then mixed with 200 mM triethylammonium bicarbonate buffer and u sonicated. Then, trypsin was added at a ratio of 1:50 (protease/protein, m/m), and the ture was enzymolyzed overnight. The next day, 5 mM dithiothreitol was added, and mixture was reduced at 56 °C for 30 min, followed by the addition of 11 mM iodoace ide and incubation at room temperature for 15 min in the dark.

TMT Labeling
After hydrolysis with trypsin as depicted in Figure 1, we desalted the peptides Strata™ X C18 resin (Phenomenex, Torrance, CA, USA), freeze-dried in vacuum, then solved in 0.5 M triethylammonium bicarbonate buffer and labeled using a TMT10plex Isobaric Label Reagent Set (Thermo Fisher Scientific, Waltham, MA, USA) in accord with the manufacturer's instructions. In brief, the labeled reagent was defrosted, solved in acetonitrile, mixed with the peptides, and incubated for 2 h at room tempera The labeled peptides were mixed, desalted, and freeze-dried in vacuum.

TMT Labeling
After hydrolysis with trypsin as depicted in Figure 1, we desalted the peptides with Strata™ X C18 resin (Phenomenex, Torrance, CA, USA), freeze-dried in vacuum, then dissolved in 0.5 M triethylammonium bicarbonate buffer and labeled using a TMT10plexTM Isobaric Label Reagent Set (Thermo Fisher Scientific, Waltham, MA, USA) in accordance with the manufacturer's instructions. In brief, the labeled reagent was defrosted, dissolved in acetonitrile, mixed with the peptides, and incubated for 2 h at room temperature. The labeled peptides were mixed, desalted, and freeze-dried in vacuum.

High-pH Peptide Fractionation and Enrichment of Phosphopeptides
High-pH reversed-phase liquid chromatography (RPLC) with a ZORBAX 300Extend-C18 column (particle size, 5 µm; inner diameter, 4.6 mm; length, 250 mm; Agilent Technologies, Inc., Santa Clara, CA, USA) was used to fractionate the peptides. Then, the components were separated with acetonitrile (gradient, 8-32% at pH 9.0) over 60 min into 60 fractions, which were combined into nine fractions and freeze-dried in a vacuum chamber for subsequent analysis.
Peptide mixtures were incubated with a suspension of prewashed immobilized metal chelate affinity chromatography (IMAC) microspheres in a loading buffer (50% acetonitrile/6% trifluoroacetic acid) while vibrating. To elute nonspecifically adsorbed peptides, the IMAC microspheres were successively washed with 50% acetonitrile/6% trifluoroacetic acid and 30% acetonitrile/0.1% trifluoroacetic acid. Afterward, an elution buffer containing Insects 2022, 13, 862 5 of 23 10% NH 4 OH was added to elute the enriched phosphopeptides from the IMAC microspheres while vibrating. The supernatant containing the phosphopeptides was collected and lyophilized for LC-MS/MS analysis.

LC-MS/MS Analysis
The peptides from the holoprotein and enriched phosphopeptides were dissolved in solvent A (0.1% formic acid in 2% acetonitrile) and solvent B (0.1% formic acid in 90% acetonitrile) at a constant flow rate of 500 nL/min with an EASY-nLC 1200 UPLC system. The peptides from the holoproteins were separated with linear gradients of 5-22% solvent B for 36 min, 22-32% solvent B for 16 min, 32-80% solvent B for 4 min, and 80% solvent B for 4 min. The phosphopeptides were separated with linear gradients of 4-20% solvent B for 40 min, 20-32% solvent B for 12 min, 32-80% solvent B for 4 min, and 80% solvent B for 4 min. Afterward, the peptides were subjected to a nanospray ionization source followed by tandem MS with a Q ExactiveTM HF-X Hybrid Quadrupole-OrbitrapTM MS system (Thermo Fisher Scientific) coupled to the EASY-nLC 1200 system. The electrospray voltage was set at 2.2 kV. For peptides from the holoproteins, the m/z scan range was 350 to 1600 with a resolution of 120,000 and a second-order mass scan resolution of 30,000. For the enriched phosphopeptides, the m/z scan range was 350 to 1400 with a resolution of 60,000 and a second-order mass scan resolution of 30,000. In data-dependent collection mode, the 20 most massive precursors were chosen for MS/MS fragmentation by higher-energy collisional dissociation at 28% normalized collision energy. The automatic gain control was set at 3E6 to avoid overfilling the ion trap. The maximum inject time was set to 50 ms, and the dynamic exclusion time was set to 15 s.

Database Search
The MS/MS data that were generated were analyzed by Jingjie PTM Biolab (Hangzhou, China) using a MaxQuant (v.1.5.2.8) quantitative proteomics software package (https: //maxquant.org/, accessed on 25 December 2020). Tandem mass spectra were searched against the UniProt database (B. terrestris, 17,032 sequences; https://www.uniprot.org/, accessed on 29 December 2020) concatenated with the reverse decoy database. The UniProt Gene Ontology Annotation database (http://www.ebi.ac.uk/GOA, accessed on 13 March 2021) was used to annotate the GO terms assigned to the proteins, which were classified as cellular components, molecular functions, and physiological processes.

Motif Analysis
Motif-X software (http://motif-x.med.harvard.edu/motif-x.html, accessed on 26 September 2021) was used to examine sequence models composed of amino acids modifying the particular 13-mer locations (six amino acids upstream and downstream) in all protein sequences. All of the protein sequences were searched using background database parameters and other default parameters.
Insects 2022, 13, 862 6 of 23 2.9. Kinase Analysis GPS 5.0 software (http://gps.biocuckoo.cn/userguide.php, accessed on 26 September 2021) was used to predict kinase substrate regulation based on the theory that main specificity is provided by short linear motifs around phosphorylation sites. The corresponding kinases were identified by comparison with the kinase sequences in the IEKPD2.0 database (http://iekpd.biocuckoo.org, accessed on 26 September 2021). Potential false-positive hits were filtered in reference to protein-protein interactions (PPIs). A "medium" threshold was chosen for use with the GPS 5.0 software.
The regulatory state of a kinase is reflected by changes to the phosphorylation levels of substrate sites. To predict kinase activity, a gene set enrichment analysis was employed, using log-transformed phosphorylation levels (or ratio values) as a rank file and kinasephosphorylation site regulations structured as a gmt file for each sample or comparable group. Normalized enrichment scores were regarded as kinase activity scores. Each kinase was considered to have positive activity if the predominant change to the substrate was an increase in phosphorylation and vice versa.
Numerous substrates may be controlled by the same kinase, and a phosphorylation site can be regulated by multiple kinases. Based on the complicated regulatory relationships, a kinase-substrate regulatory network was built using kinases that were predicted to have positive or negative activity and significantly differentially phosphorylated sites.

PPI Networks
All DEPs were searched against the STRING database (version 11.0; https://string-db. org/, accessed on 28 September 2021) for PPIs. Only interactions involving the searched dataset were selected, thereby excluding external candidates. STRING defines the metric "confidence score" as interaction confidence. All interactions with confidence scores ≥ 0.7 (high confidence) were retrieved. The PPI networks were visualized with the R package "networkD3" (https://cran.r-project.org/web/packages/networkD3/index.html, accessed on 28 September 2021).

Quantitative Proteomic and Phosphoproteomic Analyses of B. terrestris Queens in the D, PD, FPD Stages
The 10-plex TMT isobaric labeling quantitative proteomic method with high-pH reversed phase (high-pH-RP) chromatography coupled with LC-MS/MS was used to investigate and compare differences in protein the abundance and phosphorylation status of the three diapause stages of B. terrestris. The principal component analyses (PCA) of both omics datasets showed good clustering of the biological replicates as well as clear separation among the three treatments (Figure 2A,B). In total, 34,760 unique peptides from 4655 proteins were identified, with average sequence coverage of 22% ( Figure 2C, Supplementary Tables S1 and S2). The phosphopeptide dataset was first normalized based on the total abundance and further normalized to the corresponding protein abundance for each phosphopeptide. Then, 10,072 unique phosphopeptides were identified from 3339 phosphoproteins with 10,600 phosphorylation sites ( Figure 2D, Supplementary Tables S3 and  S4). In addition, 1927 proteins overlapped between the proteome and phosphoproteome ( Figure 2E). Venn diagrams demonstrate that there were more DEPs in FPD vs. PD than PD vs. D ( Figure 2F), with opposite results for DEPPs ( Figure 2G). Frequency distributions show that the number and magnitude of changes (log2 ratio) were more significant to the phosphoproteome than the proteome ( Figure 2H). The distribution patterns of statistical significance (vertical axis) and log2 fold-change (horizontal axis) for all proteins (Supplementary Figure S1A) and phosphoproteins (Supplementary Figure S1B) identified between the different groups are visualized using volcano plots. Collectively, these results confirmed that phosphorylation plays a crucial role in all stages of diapause.
show that the number and magnitude of changes (log2 ratio) were more significant to the phosphoproteome than the proteome ( Figure 2H). The distribution patterns of statistical significance (vertical axis) and log2 fold-change (horizontal axis) for all proteins (Supplementary Figure S1A) and phosphoproteins (Supplementary Figure S1B) identified between the different groups are visualized using volcano plots. Collectively, these results confirmed that phosphorylation plays a crucial role in all stages of diapause.

Expression Profile of Proteins during Different Diapause Stages
Most of the peptides were comprised of 7-20 amino acid residues, and the charge numbers were between 2 and 5 (Supplementary Figure S2A), while the molecular weights

Expression Profile of Proteins during Different Diapause Stages
Most of the peptides were comprised of 7-20 amino acid residues, and the charge numbers were between 2 and 5 (Supplementary Figure S2A), while the molecular weights varied (Supplementary Figure S2B). A heatmap of Pearson's correlation coefficients indicated that the proteome was highly reproducible (Supplementary Figure S2C). Collectively, these results indicated that the MS data matched the quality control requirements. The relative quantity of the identified proteins was considered up-or downregulated when the ratio was >1.3 or <0.77 (p < 0.05). The DEPs of the three groups are summarized in Supple- To elucidate the protein expression patterns during the three stages of diapause, an Mfuzz clustering analysis was performed based on the abundance levels. In order to screenout proteins with significant changes, the relative expression levels were first converted to log2 values, and those with SD > 0.4 were screened. The Mfuzz expression pattern clustering analysis (clustering ambiguity "m" = 2) divided the remaining 851 proteins into six clusters (Supplementary Table S5). The results of the GO/KEGG/domain enrichment analyses of proteins in different clusters are shown in Supplementary Figure S6. The representative entries in the GO/KEGG/domain enrichment results are displayed on the right side of the corresponding cluster ( Figure 3). The expression levels of 153 and 205 proteins in clusters 1 and 2 were upregulated from D to PD, respectively ( Figure 3). The KEGG analysis showed that most of the upregulated proteins correlated to pathways of ribosomes, protein digestion and absorption, steroid hormone biosynthesis, and starch and sucrose metabolism (Supplementary Figure S6A). The biological processes of these proteins are mainly associated with the regulation of transcription, ribosome assembly and biogenesis, cytoplasmic translation, and the post-transcriptional regulation of gene expression (Supplementary Figure S6B). The cellular component analysis showed that these proteins were mainly related to ribosomes, nuclei, and chromosomes (Supplementary Figure S6C). These different genes revealed the reactivation of gene expression and protein translation that emerges from a state of developmental arrest, low metabolic rates, and depressed levels of mRNA. Meanwhile, 89 and 123 proteins in clusters 3 and 5 were downregulated from D to PD, respectively ( Figure 3). Most of the downregulated proteins were associated with valine/leucine/isoleucine/tryptophan metabolism and degradation, fatty acid degradation, and glycolysis/gluconeogenesis (Supplementary Figure S6A). Collectively, these results indicated that protein synthesis and energy metabolism increased as the B. terrestris queens released from diapause. The most significantly up-and downregulated proteins from PD to FPD were grouped into clusters 6 and 4, respectively ( Figure 3). The proteins of cluster 6 were mainly related to fat digestion and absorption, the biosynthesis of unsaturated fatty acids, fatty acid elongation, and protein processing in endoplasmic reticulum, and most of these proteins were located in endoplasmic reticulum (Supplementary Figure S6A,C). Meanwhile, the KEGG analysis showed that the proteins in cluster 4 were mainly associated with oxidative phosphorylation, the thyroid hormone signaling pathway, and retrograde endocannabinoid signaling (Supplementary Figure S6A). The cellular component analysis showed that these proteins are mainly related to mitochondria (Supplementary Figure S6C). Together, these results suggested that enhanced synthesis of unsaturated fatty acid and protein processing in combination with decreased oxidative phosphorylation may contribute to ovarian maturation in B. terrestris.
There were 21 (7 upregulated and 14 downregulated) and 16 upregulated proteins with a fold change > 5 (p < 0.05) in the PD/D and FPD/PD groups, respectively. The upregulated DEPs are listed in Table 1. The greatest increase in abundance in the PD/D group (11.5-fold) was observed for a receptor protein associated with the development and maintenance of the mammalian nervous system (LOC100648620: leucine-rich repeat and immunoglobulin-like domain-containing nogo receptor-interacting protein 2) [35]. The remining proteins were annotated with putative roles in the cell cycle (LOC100642207: CDKactivating kinase assembly factor MAT1), cell adhesion and proliferation (LOC100647337: adenomatous polyposis coli), serum Pi level regulation (LOC100647822: inorganic phosphate cotransporter), actin cytoskeleton regulation (LOC100650848: disheveled-associated activator of morphogenesis), protein processing (LOC100643542: E3 ubiquitin-protein ligase Rnf220), and immunity (LOC110119488: transmembrane protease serine). The protein with the greatest increase in the FPD/PD group (10.6-fold) was observed for a calcium-activated chloride channel protein (LOC100646632: anoctamin-1) [36]. Certain immune-associated proteins included serine protease inhibitor 3/4 (LOC100652301), p53 protein (LOC100646232), CLIP domain-containing serine protease (LOC100652157), and Insects 2022, 13, 862 9 of 23 TNF receptor-associated factor 5 (LOC105665677). Other proteins with increased expression of at least 5-fold included two lipid metabolic proteins (LOC100645800: pancreatic lipase-related protein 2 and LOC100651730: phospholipase A2), two venom proteins (LOC100642484: venom protease and LOC100647178: venom acid phosphatase Acph-1), two glycometabolism-related proteins (LOC100643608 and LOC100644978), and a major royal jelly protein (LOC100648898). There were 21 (7 upregulated and 14 downregulated) and 16 upregulated proteins with a fold change > 5 (p < 0.05) in the PD/D and FPD/PD groups, respectively. The upregulated DEPs are listed in Table 1. The greatest increase in abundance in the PD/D group (11.5-fold) was observed for a receptor protein associated with the development and maintenance of the mammalian nervous system (LOC100648620: leucine-rich repeat and immunoglobulin-like domain-containing nogo receptor-interacting protein 2) [35]. The remining proteins were annotated with putative roles in the cell cycle (LOC100642207:  The expression levels of 14 proteins were significantly reduced in the PD/D group (Supplementary Table S6). Two of these proteins play structural roles (LOC100642449: titin and LOC100651429: basement membrane-specific heparan sulfate proteoglycan core protein). The remaining proteins were associated with gene expression regulation (LOC100651892: trithorax; LOC100644496: LIM domain-containing protein jub; and LOC100648151: ataxin-2-like protein) and membrane or channel trafficking (LOC100647837 and LOC100650948). Other proteins with expression levels reduced by at least 5-fold included BMP-binding endothelial regulator (LOC100652150), putative serine/threonine-protein kinase samkC (LOC100649678), immunity-related protein CDV3 (LOC100642536), a storage protein (LOC100650745: hexamerin), vitellogenin-like precursor (LOC100643258), and a protein of unknown function (LOC100644429). In addition, 14 cuticle proteins were significantly decreased in the PD or FPD stage (Supplementary Table S12).

Functional Analysis of the Identified Phosphoproteins
The length and charge distribution of the identified phosphoproteins met the quality control requirements (Supplementary Figure S7A). A heatmap of the Pearson's correlation coefficients indicated that the phosphoproteomes were highly reproducible (Supplementary Figure S7B). Among these, 2241 phosphopeptides corresponding to 1218 proteins, 1021 phosphopeptides corresponding to 692 proteins, and 2152 phosphopeptides corresponding to 1262 proteins were upregulated by more than 1.3-fold between PD/D, FPD/PD, and FPD/D, respectively (Supplementary Figure S7C). To elucidate functional correlations, the DEPPs were enriched by the KEGG pathway, domain, and GO classifications analyses, and the results are showed in Supplementary Figures S8-S10, respectively.
An Mfuzz clustering analysis was also performed to elucidate the abundance patterns of the phosphoproteins. In total, 3572 phosphoprotein sites were screened (Supplementary Table S7). The results of the GO/KEGG/domain enrichment analyses of the phosphoproteins among the different classifications are shown in Supplementary Figure S11. The proteins assigned to clusters 2 and 3 exhibited increased phosphorylation from D to the PD stage and were selectively enriched for several KEGG terms, including signal transduction (e.g., the calcium signaling pathway, the cyclic guanosine monophosphate (cGMP)-protein kinase G (PKG) signaling pathway, and the hedgehog signaling pathway), neural signal transmission (e.g., the glutamatergic synapse, the GABAergic synapse, and the serotonergic synapse), hormone biosynthesis and secretion (e.g., gonadotropin-releasing hormone secretion, insulin secretion, thyroid hormone synthesis, and cortisol synthesis and secretion), and energy-related metabolism (e.g., the tricarboxylic acid (TCA) cycle) (Figure 4, Supplementary Figure S11A). The detailed GO terms revealed that the biological processes of proteins classified to clusters 2 and 3 were mainly associated with chromosome organization and the regulation of transcription by RNA polymerase II, which is consistent with the changes to the proteome that suggested the reactivation of gene expression (Supplementary Figure S11B). In addition, the enrichment of muscle thin filament assembly and the citrate metabolic process indicated the reactivation as the B. terrestris queens released from diapause (Supplementary Figure  S11B). The proteins assigned to clusters 1, 4, and 5, which exhibited dephosphorylation from D to the PD stage, were enriched in the KEGG term (mitogen-activated protein kinase) MAPK signaling pathway (Figure 4, Supplementary Figure S11A). Proteins with increased phosphorylation levels from PD to FPD were over-represented in cluster 6 and involved in RNA transport, and the biological process was related to the cell cycle (Figure 4, Supplementary Figures S5B and S11A). The KEGG analysis of all DEPPs in Supplementary  Figure S8 showed that the upregulated phosphorylated proteins from PD to FPD were involved in energy-related metabolism, including glycerolipid metabolism, glyoxylate and dicarboxylate metabolism, and glycan degradation. Meanwhile, the downregulated phosphorylated proteins were associated with hormone synthesis and secretion, muscle contraction, and the calcium signaling pathway. In total, 86 (19 upregulated and 67 downregulated) and 10 (9 upregulated and 1 downregulated) phosphorylation sites had fold changes >5 in the PD/D and FPD/PD groups, respectively. The specific upregulated phosphorylation sites are shown in Table  2. Of the 19 phosphorylation sites significantly increased by at least 5-fold in the PD/D groups, 5 phosphorylated proteins were annotated with roles in cytoskeleton and muscle proteins (LOC100648462: catenin; LOC100647948: microtubule-actin cross-linking factor 1; LOC100646259: trichohyalin; LOC100647343: myosin heavy chain; and LOC100643650: titin). Three proteins were annotated as chromatin regulators (LOC100652289: RNA In total, 86 (19 upregulated and 67 downregulated) and 10 (9 upregulated and 1 downregulated) phosphorylation sites had fold changes >5 in the PD/D and FPD/PD groups, respectively. The specific upregulated phosphorylation sites are shown in Table 2. Of the 19 phosphorylation sites significantly increased by at least 5-fold in the PD/D groups, 5 phosphorylated proteins were annotated with roles in cytoskeleton and muscle proteins (LOC100648462: catenin; LOC100647948: microtubule-actin cross-linking factor 1; LOC100646259: trichohyalin; LOC100647343: myosin heavy chain; and LOC100643650: titin). Three proteins were annotated as chromatin regulators (LOC100652289: RNA polymerase-associated protein LEO1; LOC100644519: nucleolin; and LOC100648929: lysinespecific demethylase lid). Additional proteins with significantly increased phosphorylation levels included two kinases (LOC100646154: tau-tubulin kinase homolog Asator and LOC100645231: serine/threonine-protein kinase MARK2), two transporter proteins (LOC100649391: solute carrier family 25 member 38 and LOC105665975: golgin subfamily A member 4), transcription-and translation-associated proteins (LOC100647181: glycosyltransferase subunit STT3B; LOC100643117: splicing factor 3B subunit 2; LOC100642874: 60S ribosomal protein L5; and LOC100648934: 60S ribosomal protein L15), and a key component of the ubiquitin ligase protein (LOC100644051: DDB1-and CUL4-associated factor 8). The 67 downregulated phosphorylation sites in the PD/D groups are shown in Supplementary Table S8. These 67 phosphorylation sites of 62 proteins were annotated with putative roles in signal transduction mechanisms (n = 14), the cytoskeleton (n = 8), intracellular trafficking and vesicular transport (n = 7), RNA processing and modification (n = 5), transcription (n = 4), carbohydrate transport and metabolism (n = 4), post-translational modification, protein turnover, and chaperones (n =4) as well as proteins of unknown function and unannotated proteins (n = 17) (Supplementary Table S8). Of the nine phosphorylation sites significantly increased by at least 5-fold in the FPD/PD groups, four phosphorylated proteins were annotated with roles in neuroendocrine function (LOC100642774: serine/arginine repetitive matrix protein 2; LOC105666848: jerky protein; LOC100650297: synaptic vesicle membrane protein; and LOC100645647: Katnb1_0 protein), two proteins had roles in sugar metabolism (LOC100644377: sugar transporter ERD6-like 8 and LOC100644978: facilitated trehalose transporter Tret1-like), two proteins had roles in ovarian development (LOC100650436: vitellogenin and LOC100645056: TBC1 domain family member 30), and one protein was annotated as heparanase (LOC100646845).
It is generally believed that the regulatory state of a kinase is reflected by a change in the phosphorylation levels of its substrate sites. Phosphokinase activity was predicted, and the statistics are shown in Supplementary Figure S12. A predictive analysis of the kinase activity showed that 11 kinases, including LOC100647871 (serine/threonine kinase 11 (STK11/LKB1)), LOC100647145 (CDK7), and LOC100631091 (JNK), were strongly activated in the D stage, while four and seven kinases were activated in the PD and FPD stages, respectively ( Figure 5A). LOC100647871 (STK11/LKB1) was predicted to phosphorylate A0A6P3U9P6 (Forkhead box protein O (FoxO)) at residues Thr 53 and Thr 205, A0A6P3DMX8 (5 adenosine monophosphate-activated protein kinase (AMPK) subunit γ-2) at residues Ser 639 and Ser 691, A0A6P3DLX9 (AMPK subunit β-1) at residue Ser 55, and A0A6P3UAT7 (AMPK subunit α-2) at residue Thr 175 ( Figure 5C). From D to PD, the activity of LOC100647871 (STK11/LKB1) was blocked, and the phosphorylation level was downregulated.
It is generally believed that the regulatory state of a kinase is reflected by a change in the phosphorylation levels of its substrate sites. Phosphokinase activity was predicted, and the statistics are shown in Supplementary Figure S12. A predictive analysis of the kinase activity showed that 11 kinases, including LOC100647871 (serine/threonine kinase 11 (STK11/LKB1)), LOC100647145 (CDK7), and LOC100631091 (JNK), were strongly activated in the D stage, while four and seven kinases were activated in the PD and FPD stages, respectively ( Figure 5A). LOC100647871 (STK11/LKB1) was predicted to phosphorylate A0A6P3U9P6 (Forkhead box protein O (FoxO)) at residues Thr 53 and Thr 205, A0A6P3DMX8 (5′ adenosine monophosphate-activated protein kinase (AMPK) subunit γ-2) at residues Ser 639 and Ser 691, A0A6P3DLX9 (AMPK subunit β-1) at residue Ser 55, and A0A6P3UAT7 (AMPK subunit α-2) at residue Thr 175 ( Figure 5C). From D to PD, the activity of LOC100647871 (STK11/LKB1) was blocked, and the phosphorylation level was downregulated.

Motif Analysis
Protein kinases prefer specific substrates with conserved motifs. An analysis with Motif-X software predicted 21 upregulated and 23 downregulated conserved motifs among all of the identified phosphorylation sites in the PD vs. D group and 12 upregulated and eight downregulated motifs in the FPD vs. PD group (Supplementary Table S11).
Among the upregulated phosphorylation motifs in PD vs. D, the enriched motif contained the acidic amino acid sequence Ser-Asp-x-Glu (the serine was the phosphorylated amino acid residue), which had the highest motif score of 32.00 and the highest fold increase of 11.0 (Table 3). Among the downregulated phosphorylation motifs in PD vs. D, the prolinedirected sequence Thr-Pro-Pro (the threonine was the phosphorylated amino acid residue) was strikingly decreased, with a motif score of 27.44 and a fold decrease of 38.2. According to the predicted kinase-substrate regulations, the " . . . TPP . . . " sequence was predicted to be recognized by CDK7 (LOC100647145) ( Table 3). For the upregulated phosphopeptides in FPD vs. PD, the most statistically significant and enriched motif was also Thr-Pro-Pro, with a motif score of 24.20 and a fold increase of 39.9, which was recognized by CDK1 (LOC100651604) ( Table 3). The most significantly downregulated motif in FPD vs. PD was the basic amino acid sequence Lys-x-x-Ser-Pro (Table 3). Heatmaps showing the enrichment and depletion of specific amino acids neighboring the serine/threonine phosphosites were also generated (Supplementary Figure S13). Among the upregulated phosphorylation motifs in PD vs. D, the enriched motif contained the acidic amino acid sequence Ser-Asp-x-Glu (the serine was the phosphorylated amino acid residue), which had the highest motif score of 32.00 and the highest fold increase of 11.0 (Table 3). Among the downregulated phosphorylation motifs in PD vs. D, the proline-directed sequence Thr-Pro-Pro (the threonine was the phosphorylated amino acid residue) was strikingly decreased, with a motif score of 27.44 and a fold decrease of 38.2. According to the predicted kinase-substrate regulations, the "…TPP…" sequence was predicted to be recognized by CDK7 (LOC100647145) ( Table 3). For the upregulated phosphopeptides in FPD vs. PD, the most statistically significant and enriched motif was also Thr-Pro-Pro, with a motif score of 24.20 and a fold increase of 39.9, which was recognized by CDK1 (LOC100651604) ( Table 3). The most significantly downregulated motif in FPD vs. PD was the basic amino acid sequence Lys-x-x-Ser-Pro (Table 3). Heatmaps showing the enrichment and depletion of specific amino acids neighboring the serine/threonine phosphosites were also generated (Supplementary Figure S13).

PPI and Module Analysis of the DEPs and DEPPs
To further clarify the complicated biological processes regulated in the D, PD, and FPD stages, PPI networks were established based on the identified DEPs and DEPPs. PPI Among the upregulated phosphorylation motifs in PD vs. D, the enriched motif con-tained the acidic amino acid sequence Ser-Asp-x-Glu (the serine was the phosphorylated amino acid residue), which had the highest motif score of 32.00 and the highest fold increase of 11.0 (Table 3). Among the downregulated phosphorylation motifs in PD vs. D, the proline-directed sequence Thr-Pro-Pro (the threonine was the phosphorylated amino acid residue) was strikingly decreased, with a motif score of 27.44 and a fold decrease of 38.2. According to the predicted kinase-substrate regulations, the "…TPP…" sequence was predicted to be recognized by CDK7 (LOC100647145) ( Table 3). For the upregulated phosphopeptides in FPD vs. PD, the most statistically significant and enriched motif was also Thr-Pro-Pro, with a motif score of 24.20 and a fold increase of 39.9, which was recognized by CDK1 (LOC100651604) ( Table 3). The most significantly downregulated motif in FPD vs. PD was the basic amino acid sequence Lys-x-x-Ser-Pro (Table 3). Heatmaps showing the enrichment and depletion of specific amino acids neighboring the serine/threonine phosphosites were also generated (Supplementary Figure S13).

PPI and Module Analysis of the DEPs and DEPPs
To further clarify the complicated biological processes regulated in the D, PD, and FPD stages, PPI networks were established based on the identified DEPs and DEPPs. PPI tained the acidic amino acid sequence Ser-Asp-x-Glu (the serine was the phosphorylated amino acid residue), which had the highest motif score of 32.00 and the highest fold increase of 11.0 (Table 3). Among the downregulated phosphorylation motifs in PD vs. D, the proline-directed sequence Thr-Pro-Pro (the threonine was the phosphorylated amino acid residue) was strikingly decreased, with a motif score of 27.44 and a fold decrease of 38.2. According to the predicted kinase-substrate regulations, the "…TPP…" sequence was predicted to be recognized by CDK7 (LOC100647145) ( Table 3). For the upregulated phosphopeptides in FPD vs. PD, the most statistically significant and enriched motif was also Thr-Pro-Pro, with a motif score of 24.20 and a fold increase of 39.9, which was recognized by CDK1 (LOC100651604) ( Table 3). The most significantly downregulated motif in FPD vs. PD was the basic amino acid sequence Lys-x-x-Ser-Pro (Table 3). Heatmaps showing the enrichment and depletion of specific amino acids neighboring the serine/threonine phosphosites were also generated (Supplementary Figure S13).

PPI and Module Analysis of the DEPs and DEPPs
To further clarify the complicated biological processes regulated in the D, PD, and FPD stages, PPI networks were established based on the identified DEPs and DEPPs. PPI tained the acidic amino acid sequence Ser-Asp-x-Glu (the serine was the phosphorylated amino acid residue), which had the highest motif score of 32.00 and the highest fold increase of 11.0 (Table 3). Among the downregulated phosphorylation motifs in PD vs. D, the proline-directed sequence Thr-Pro-Pro (the threonine was the phosphorylated amino acid residue) was strikingly decreased, with a motif score of 27.44 and a fold decrease of 38.2. According to the predicted kinase-substrate regulations, the "…TPP…" sequence was predicted to be recognized by CDK7 (LOC100647145) ( Table 3). For the upregulated phosphopeptides in FPD vs. PD, the most statistically significant and enriched motif was also Thr-Pro-Pro, with a motif score of 24.20 and a fold increase of 39.9, which was recognized by CDK1 (LOC100651604) ( Table 3). The most significantly downregulated motif in FPD vs. PD was the basic amino acid sequence Lys-x-x-Ser-Pro (Table 3). Heatmaps showing the enrichment and depletion of specific amino acids neighboring the serine/threonine phosphosites were also generated (Supplementary Figure S13).

PPI and Module Analysis of the DEPs and DEPPs
To further clarify the complicated biological processes regulated in the D, PD, and FPD stages, PPI networks were established based on the identified DEPs and DEPPs. PPI

PPI and Module Analysis of the DEPs and DEPPs
To further clarify the complicated biological processes regulated in the D, PD, and FPD stages, PPI networks were established based on the identified DEPs and DEPPs. PPI networks for the PD vs. D group are shown in Figure 6, and those for the FPD vs. PD and FPD vs. D groups are shown in Supplementary Figure S14. The use of the MCODE plugin identified four highly interconnected clusters based on a proteomic analysis of the PD vs. D network that were related to ribosomes, spliceosomes, fructose and mannose metabolism, and glycolysis/gluconeogenesis. Three clusters were identified among the phosphorylated proteins, including proteins related to ribosomes/ribosome biogenesis, oxidative phosphorylation, and metabolic pathways. Taken together, the networks inferred from the proteomic and phosphoproteomic analyses suggest disturbances in transcription, translation, and glycometabolism. D network that were related to ribosomes, spliceosomes, fructose and mannose metabolism, and glycolysis/gluconeogenesis. Three clusters were identified among the phosphorylated proteins, including proteins related to ribosomes/ribosome biogenesis, oxidative phosphorylation, and metabolic pathways. Taken together, the networks inferred from the proteomic and phosphoproteomic analyses suggest disturbances in transcription, translation, and glycometabolism.

Discussion
The aim of this study was to create a global overview of the mechanisms underlying diapause, which allows organisms to adapt to unfavorable conditions, with the use of tag-based multiplex proteomic/phosphoproteomic quantitative analyses. The results demonstrated that the termination of diapause from the D to PD stages was closely associated with the neuroendocrine system, energy metabolism, steroid hormone biosynthesis, and various signaling transduction pathways. The transition from PD to the FPD stage was characterized by various biosynthesis and catabolism processes. Kinases and interaction networks were also characterized. According to our best knowledge, this is the first large-scale proteome/phosphoproteome analysis of diapause in bumblebees, and it provides a comprehensive overview of a dynamically intertwined protein and phosphoprotein interaction network.
There have been several studies of the transcriptome and proteome of the different diapause stages of B. terrestris. Comparative analyses indicated that diapause in bumblebees is associated with nutrient storage, core metabolic pathways, stress resistance, insect hormone biosynthesis, and proteins involved in cuticle maintenance [37][38][39]. In the current study, comparisons of the protein expression profiles during the three stages of diapause were performed using TMT-labeled proteomics to expand the current understanding of changes at the phosphorylation level.
In arthropods, growing evidence indicates the roles of neurotransmitters in the regulation of diapause. For example, glutamate is reportedly involved in triggering diapause memory in the induction phase of the cotton bollworm (H. armigera) [40]. Key proteins involved in neuroactive ligand-receptor interactions and the glutamatergic synapse pathway related to glutamate neurotransmission were upregulated in early diapause in the female red spider mite (Tetranychus urticae) [41]. The expression levels of glutamate decarboxylase, which plays a role in neurotransmitter synthesis, were increased in the hemolymph of bumblebee queens at 48 h compared to 6 h postdiapause [39]. In this study, many genes associated with neurotransmitter synthesis were regulated at both the protein and phosphorylation levels during different stages of diapause. KEGG and Mfuzz clustering analyses both revealed that genes involved in neuroactive ligand-receptor interactions (map04080) were over-represented among the upregulated proteins in the PD vs. D group. Compared to the PD to D stage phosphorylation level, the glutamatergic synapse (map04724) and retrograde endocannabinoid signaling (map04723) were over-represented among the upregulated phosphorylated proteins. Indeed, in this study, the phosphorylation levels increased for glutamine synthetase, which is critical for effective neurotransmission; two G-protein subunits, guanine nucleotide-binding protein G (Q) subunit alpha and G(I)/G(S)/G(T) subunit beta-1, which are necessary for G-protein coupling to metabotropic glutamate receptors; and one excitatory amino acid transporter, which is the major transport mechanism for extracellular glutamate in the nervous system [42]. In addition, the GABAergic synapse and the serotonergic synapse were enriched in cluster 2 and upregulated from D to the PD stage, as determined by the phosphoproteomic analysis. Proteins with notable changes in the neuroendocrine system are potential candidates for future mechanistic studies of their upstream roles in insects.
MAPK plays an important role in insect growth, development, and immunity. Rapid condensation, which may be related to MAPK, is important for insets to quickly adapt to low temperatures. A previous study reported that p38, a main member of the MAPK family, regulated the rapid onset and termination of temperature-induced diapause in the flesh fly (S. crassipalpis) [43]. At 0 • C, p38 is rapidly activated, consistent with the temperature required for rapid condensation. However, at 25 • C, the phosphorylation of p38 is decreased, indicating that p38 is critical for the rapid response of S. crassipalpis to low temperatures [44]. Extracellular signal-regulated kinase (ERK), another member of the MAPK family, is activated by exposure to cold temperatures and regulates the termination of diapause in the silkworm (B. mori) [45]. The expression levels of phosphorylated ERK remain high in the prediapause period, then decrease in the diapause period of the false melon beetle (Atrachya menetriesi). Phosphorylated ERK decreases as the temperature increases from 7.5 • C to 25 • C [24]. The results of the present study showed that the phosphorylation level of the MAPK signaling pathway of bumblebees decreased when the temperature was increased from 2 • C (D) to 28 • C (PD), consistent with previous studies. The phosphorylation of RAC (LOC100643457) at residues Ser 157 and Thr 362, PKC (LOC100643492) at residue Thr 57, PAK1 (LOC105666875) at residues Ser 62 and Ser 192, PAK3 (LOC100652006) at residues Ser 2 and Ser 256, MAPKKK15 (LOC100644945) at residue Ser 906, MAPKK (LOC100631092) at residues Ser 689, Ser 419, Ser 701, Ser 245, Ser 714, and Thr 754, and RAF (LOC100645437) at residue Ser 338 were significantly decreased from D to the PD phase. Hence, all these sites are worthy of molecular mechanism studies in the future.
Generally, highly active mitochondria and strong oxidative phosphorylation are hallmarks of growing organisms. However, metabolic quiescence leads to a switch to glycolysis and related gluconeogenesis [46,47]. In the present study, the TCA cycle was enriched in the upregulated phosphorylation group from D to the PD stage, while glycolysis/gluconeogenesis was enriched in the downregulated protein group. These results are consistent with those of previous studies and suggest that increased phosphorylation indicates the activation of oxidative phosphorylation and the TCA cycle. From PD to the FPD stage, oxidative phosphorylation and the TCA cycle were over-represented in both the downregulated proteins and phosphoproteins, which suggests decreased energy demand.
In addition, cuticle proteins were also identified among the DEPs and DEPPs. Of the 24 identified proteins, 14 were significantly decreased at the PD or FPD stage at the protein level. Among the 24 cuticle proteins, three proteins were phosphorylated at 15 sites, including 7 that were decreased at the PD or FPD stage. The characteristic expression profiles of cuticle proteins and phosphorylation modification suggest roles in the recovery from diapause. Other changes were observed to the calcium signaling pathway and the cGMP-PKG signaling pathway, which are reported to play important roles in the reproductive diapause regulation of the red spider mite (T. urticae), pupal diapause regulation in the onion fly (Delia antiqua), and the termination of diapause in the Chinese citrus fly (Bactrocera minax) [41,48,49]. The phosphorylation levels of the calcium signaling pathway and the cGMP-PKG signaling pathway were upregulated from D to PD in the present study. Moreover, other molecular pathways associated with vitamin metabolism were reported to regulate the developmental trajectory associated with diapause in the butterfly (Phalera bucephala L.), annual killifish (Austrofundulus limnaeus), and American mink (Neovison vison) [50][51][52]. In the present study, proteins associated with the metabolism of retinol (vitamin A), ascorbate (vitamin C), aldarate, and vitamin B6 were over-represented at the protein level in the PD vs. D up group, both the PD vs. D and FPD vs. PD up groups, and the FPD vs. PD up group, respectively. These results demonstrate the important roles of vitamin metabolism in the development and diapause of B. terrestris. STK11/LKB1, which is among the most interesting kinases identified, was strongly activated in the D stage but was weakly activated in the PD and FPD stages. LKB1 is a tumor suppressor that regulates multiple biological pathways, including energy metabolism, cell cycle control, and cell polarity by the direct phosphorylation of 14 different AMPK family members [53][54][55][56]. AMPK is a multicomponent enzyme complex that acts as a metabolic stress sensor. To maintain the energy balance, AMPK turns off a number of ATP-using processes once it is active. AMPK is activated by the allosteric binding of AMP, which encourages the binding of upstream kinases to increase activity [57]. Therefore, the downphosphorylation of AMPK from D to the PD stage by LKB1 facilitates the utilization of ATP to provide energy for development. This is the first report of the phosphorylation of FoxO by LKB1. In contrast, the transcription factor FoxO is reportedly involved in regulating the transcription of LKB1 [58,59]. The transcriptional activity of FoxO is tightly regulated by post-translational modifications, including methylation, acetylation, glycosylation, ubiquitination, and phosphorylation [60,61]. Many kinases (i.e., AKT, PKA, DYRK1, and GSK3) reduce the transcriptional activity of FoxO via phosphorylation [62][63][64]. In response to extracellular stimuli (e.g., insulin, growth factors, and cytokines), AKT phosphorylates and inhibits the transcriptional activity of FoxO via export from the nucleus to the cytoplasm. However, other kinases (i.e., AMPK, MST1, and p38) activate the transcriptional activity of FoxO [65][66][67]. AMPK phosphorylates FoxO at residues Ser 22, Ser 383, and Thr 649 [65,68]. The decreased phosphorylation level of FoxO at residues Thr 53 and Thr 205 from D to the PD stage by LKB1 suggested that these two sites could contribute to the transcriptional activity of downstream diapause-related genes. In addition, the down-phosphorylation of AMPK from D to the PD stage also contributed to the decreased transcriptional activity of FoxO, indicating that LKB1 may act as a switch for the diapause transition.
Our group established a platform based on nanocarrier-mediated RNA interference and the CRISPR/Cas9 system to study the molecular mechanisms of genes involved in reproductive diapause. These DEPs, DEPPs, and kinases could be candidates for RNA interference and CRISPR/Cas9 to manage the diapause of B. terrestris, such as shortening or prolonging the diapause time according to the market demand in commercial breeding. The KEGG pathway, enriched with DEPs and DEPPs, such as vitamins and insect hormones, could also be applied to B. terrestris to test whether they have a regulatory effect on diapause termination and the preoviposition time.

Conclusions
In summary, TMT-labeled quantitative proteomics and phosphoproteomics were employed to study dynamic changes in B. terrestris at three diapause stages. GO, KEGG pathway enrichment, and Mufzz clustering analyses identified various DEPs and DEPPs related to diapause termination and postdiapause development. Important kinases were identified in the D, PD, and FPD stages. The results of this study will lead to a better understanding of the molecular regulation of diapause in B. terrestris.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/insects13100862/s1, Figure S1: Volcano plots from different group comparisons in proteome (A) and phosphoproteome (B). Red dots represent fold change > 1.3 and p < 0.05, green dots represent fold change <0.77 and p < 0.05, and grey dots represent insignificantly changed proteins or phosphorylation sites; Figure Table S1: The summary information in proteomic analyses; Table S2: The list of proteins identified in proteomic analyses; Table S3: The summary information in phosphoproteomic analyses; Table S4: The list of phosphorylation sites identified in phosphoproteomic analyses; Table S5: The 851 screened proteins for Mfuzz clustering analysis; Table S6: The 14 downregulated DEPs in PD/D with fold changes higher than 5; Table S7: The 3572 screened phosphoprotein sites for Mufzz clustering analysis; Table S8: The 67 downregulated phosphorylation sites in PD/D with fold changes higher than 5; Table S9: Kinase annotation; Table S10: Kinase-substrate interaction networks in each comparable group; Table S11: Phosphorylation motifs enriched by motif-x and putative protein kinases in each comparable group; Table S12: Cuticle proteins and phosphorylation sites identified in omics.  Institutional Review Board Statement: Ethics approval was not required for this study.