Comparative transcriptome analysis of adult worker bees under short-term heat stress

High temperature affects behavior, physiology, survival, and the expression of related genes in adult honeybees. Apis mellifera is the common pollinator in greenhouse and is susceptible to high temperature stress. To further explore the molecular basis related to heat stress, we compared the transcriptome profiles of adult worker bees at 25 and 45°C, and detected the expression patterns of some differentially expressed genes (DEGs) in different tissues by q RT-PCR. Differential expression analysis showed that 277 DEGs were identified, including 167 genes upregulated and 110 genes downregulated after heat stress exposure in adult worker bees. In GO enrichment analysis, DEGs were mostly enriched for protein folding, unfold protein binding, and heme binding terms. Protein processing in endoplasmic reticulum and longevity regulating pathway-multiple species were significantly enriched in KEGG. The expression levels of 16 DEGs were consistent with the transcriptome results. The expression patterns of 9 DEGs in different tissues revealed high levels in the thorax, which was supposed that the thorax may be the most important part in the response to heat stress. This study provided valuable data for exploring the function of heat resistance-related genes.

The foraging behavior of honeybees usually occurs within the temperature range of 10-40 • C and is inhibited by high temperatures (Joshi and Joshi, 2010). 20∼25 • C is the optimal temperature range for the flight of forager bees (Zeng, 2007;Rameshwor and Resham Bahadur, 2014;Cui et al., 2022). When the external temperature exceeds 40 • C, worker bees stop foraging (Joshi and Joshi, 2010) and low the hive temperature by collecting water and rapidly fanning their wings (Nicolson, 2009;Ohashi et al., 2009). Further, heat stress triggers the expression of heatrelated genes, which are involved in a series of biological processes such as cell metabolism, protein folding, and degradation (Koo et al., 2015). At present, among the many thermal-resistance genes, the heat shock proteins (HSPs) gene family plays an extremely important role (Jerbi-Elayed et al., 2015;Bhattacharyya et al., 2019;Shih et al., 2021). In addition to HSPs, there are zinc finger proteins (ZFPs), serine/threonine protein kinases (STKs) , antimicrobial peptide genes and endocrine regulatory genes (Li et al., 2022a,b). Ambient temperature has a negative effect on foraging behavior (Al-Qarni, 2006;Blazyte-Cereskiene et al., 2010), further affecting pollinator decisions and crop pollination effects. The internal environment of facility agriculture is different from the natural environment, and there will be some problems such as large temperature differences and high temperatures (Su et al., 2017). In the greenhouse, honeybee colonies are prone to adult dysplasia and colony decline, which is unfavorable to reproducing honeybee colonies (Nguyen et al., 2013;Ma et al., 2020). Transcriptome technology has played an increasingly important role in various aspects of honeybee research, including studies of behavioral classification, oxidative stress (Southwick and Heldmaier, 1987), immune responses (Mao et al., 2013), and growth and development (Santos et al., 2021). Therefore, we compared the transcriptome of Apis mellifera workers in control groups (25 • C) and in high temperature groups (45 • C) were examined using RNAseq technology and detected the expression patterns of some differentially expressed genes (DEGs) in different tissues by q RT-PCR. Our results provided basic data to better understand the regulatory molecular mechanisms of honeybee adaptation to high temperature.

Experimental insects
Three healthy colonies (A. mellifera) were obtained from the experimental apiary of Shanxi Agricultural University (Taigu, China). Two capped brood combs were selected from each colony (two combs × three colonies), then put in an incubator at 34 ± 0.5 • C and 75 ± 5% relative humidity (RH) to hatch. The hatched bees were marked on their backs with a non-toxic dye, recorded as 1 day old (1 d), and released back to their original colony. The marked bees were collected as samples at 20 d. The marked bees were collected as samples at 20 days old (20 d). The 20 d old worker bees were used as foragers.
According to the temperature data reached over 40 • C from the greenhouse in summer and previous experiments on temperature gradients of 25, 30, 35, 40, and 45 • C (Li X. Y. et al., 2019). Two treatments were designed: normal temperature, as control (CK: 25 • C, 30% RH), and high temperature (HT: 45 • C, 30% RH). All samples were put in a 15 mL diameter container and were maintained at 30% RH for 2 h at different temperatures (25 and 45 • C). After treatment, 60 bees (10 × two temperaturetreatment × three replicates) were randomly selected for RNA-seq, 540 bees (80 × two temperature-treatment × three replicates) for expression patterns of DEGs in different tissues, and then these samples were immediately frozen in liquid nitrogen and stored at −80 • C. The head, thorax, abdomen, legs, and wings were dissected (all operations were performed on ice).

RNA-seq experiment 2.2.1. Library construction and sequencing
Total RNA was extracted using Trizol Reagent according to the manufacturer's instructions (Ambion, Foster City, CA, USA). The concentration of total RNA was assessed using an RNA NanoDrop ND-2000 spectrophotometer (NanoDrop Inc., Wilmington, DE, USA). Non-denaturing agarose gel electrophoresis and the Agilent 2100 LabChip GX (Agilent, Santa Clara, CA, USA) were used to determine the integrity of extracted nucleic acids for quality control.
Three individuals were used for each biological replicate for transcriptomic analysis. Six individual cDNA libraries from the 25 • C (T01, T02, and T03) and 45 • C (T04, T05, and T06) groups were constructed. The mRNA was fragmented and randomized hexamer primers were used to copy the full-length mRNA; reverse transcriptase was used as a template for first strand cDNA synthesis. After purification, the cDNA was connected with adenine and sequencing splitter at the 3 -ends and amplified using PCR to establish a cDNA library. The Illumina NovaSeq 4000 platform (Illumina, San Diego, CA, USA) was used for sequencing. Raw data were stored as FASTQ files and all raw sequence data were deposited in the NCBI Sequence Read Archive (SRA) under Submission ID: SUB12053792, BioProject ID: PRJNA880001.

Mapping and transcriptome annotation
Reads with ploy-N (with a ratio of "N" > 10%) and low-quality (with quality scores < 20) were removed to get the clean reads. The clean reads were mapped to the reference genome of A. mellifera using HISAT2 (Li et al., 2014). Based on the selected reference genome sequence, the mapped reads were spliced using StringTie , and the potential novel transcripts were predicted using Cufflinks (Cole et al., 2012). BLAST was used to compare the discovered new genes with NR (Deng et al., 2006), Swiss-Prot (Rolf et al., 2017), Gene Ontology (GO) (Ashburner et al., 2000),  Cluster of Orthologous Groups of proteins (COG) (Tatusov et al., 2000), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2004) databases to obtain annotation information of the new genes.

Differentially expressed gene analysis
Fragments per kilobase of exon model per million (FPKM) were estimated as a measure of transcript levels or gene expression (Love et al., 2014). The density distribution plot and box-whisper plot of FPKM was presented with R software. Differentially expressed genes (DEGs) were identified by DESeq2 (Love et al., 2014). The FDR is the corrected P-value. Genes with | fold change| > 2 and FDR < 0.05 were considered as DEGs. According to the fold change in relative expression levels between the two groups, the DEGs could be classified as upregulated or downregulated. A volcano plot was generated to visualize the results with R software. A hierarchical clustering analysis of the DEGs was performed (Yu et al., 2012), and genes with the same or similar expression patterns were clustered.

Functional enrichment analysis
To determine the functional category, the DEGs were mapped to the GO database by DAVID (Dennis et al., 2003). The annotation was refined and enriched by using the top-GO (v2.40.0) package in R software (Yu et al., 2012). The biological process (BP), molecular function (MF), and cellular components (CC) were particularly upregulated or downregulated in DEGs and were extracted and visualized via DAVID, which is generally referred to as GO analysis. The KEGG pathway was analyzed by using the KOBAS (2.0). P-value < 0.05 was the threshold value for significant enrichment results.

Protein-protein interaction (PPI) analysis
A protein-protein interaction (PPI) network of the DEGs was constructed by STRING (Franceschini et al., 2013). The constructed protein interaction network was imported into Cytoscape (Shannon et al., 2013) for visualization. The PPI consists of nodes and lines. Network nodes represent proteins, and the lines between proteins represent interactions. There are seven protein-protein interactions including neighborhood, cooccurrence, fusion, co-expression, experiments, databases, and text mining.

Quantitative real-time PCR verification
Quantitative real-time PCR (q RT-PCR) was used to verify the accuracy of RNA-seq (Livak and Schmittgen, 2001). Sixteen genes Volcano plot of differentially expressed genes (DEGs) between control (CK) and high temperature (HT).
were randomly selected for verification. Nine DEGs verified were further quantified in different tissues (i.e., head, thorax, abdomen, legs, and wings). According to the CDS region of selected genes, specific primers were designed using Primer 5 (Supplementary Table 1). The primers were synthesized by Shanghai Sangon Co., Ltd (Shanghai, China). First-strand cDNA was synthesized from total RNA using the PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Kusatsu, Japan), according to the manufacturer protocol. The 15 µL reaction consisted of 100 ng/µL diluted cDNA (6.3 µL of double distilled water dilution), 7.5 µL of SYBR R Premix Ex Taq II (TaKaRa, Japan), and 0.6 µL of each primer. The q RT-PCR thermal cycle procedure was: 95 • C predefined for 30 s, 45 cycles of denaturation at 95 • C for 5 s, annealing at 60 • C, and extension for 30 s. The fold changes for RNA-seq were estimated as base 2 logarithms of the FPKM ratio, and the fold changes for q RT-PCR were based on Cq.

Expression patterns of DEGs
Nine DEGs verified were further quantified in different tissues (i.e., head, thorax, abdomen, legs, and wings) by q RT-PCR. The reaction system was consistent with the q RT-PCR verification. The mRNA quantification was conducted following the 2 − Cq method (Livak and Schmittgen, 2001). The one-way Anova test and Duncan's multiple range test (SPSS 17.0) were used for significance analysis. Data were shown mean ± SEM. The results were plotted using GraphPad Prism 8.0.

Reference genome alignment results
The sequencing data of the RNA-seq samples were shown in Table 1. Among six samples, 109.52 Gb of clean data in total was obtained, 15.89 Gb per sample with Q30 values of 93.89% and above. Filtering to remove low-quality reads, there were an average of 134,316,160 and 113,520,461 clean reads in CK and HT ( Table 1). The clean reads for each sample were compared with the specified reference genome. The length distribution of single genes was also determined and most of the genes were > 3000 bp (Supplementary Figure 2). A total of 1,743 new genes were discovered, 148 of which were functionally annotated (Supplementary Table 2).

Differential expression gene analysis
Fragments per kilobase of exon model per million values spanned six orders of magnitude from 10 −2 to 10 3 ( Figure 1A). The distribution of the six samples was basically the same ( Figure 1B). 277 DEGs were identified between CK and HT using DESeq (| log2 FC| ≤ 2, FDR ≤ 0.05). Of these, there were 167 upregulated DEGs and 110 downregulated DEGs (Figure 2). Cluster analysis showed the expression profiles of CK and HT (Figure 3). In the upregulated DEGs, the most enriched and expressive DEGs were the HSPs (HSC70-4, HSC70-5, HSP70-Ab, HSP70-Cb, HSP90, HSP40, and L (2) efl), fatty acid synthase (FAS), BAG domain-containing protein Samui (BAG), as well as several transcripts of unknown function. In the downregulated DEGs, the most enriched DEGs were E3 ubiquitin-protein ligase RNF8 (UbE3), and insulin-like growth factor 2 mRNA-binding protein 1 (IGF2BP1). Identification and annotation information for all DEGs were provided in the Supplementary Table 3.

Enrichment analysis
We performed GO and KEGG enrichment analysis of DEGs (Figures 4, 5). In a GO enrichment analysis of 277 DEGs, 143 genes were annotated, including 25 process categories at P < 0.05 ( Table 2). In response to heat stress, DEGs were enriched for various terms in the BP, CC, and MF categories (Figure 4). DEGs were classified into three GO categories: biological processes (33.5%), cellular components (24.1%), and molecular functions (41.15%). DEGs were mostly enriched for protein folding, unfold protein binding, and heme binding terms. The DEGs coregulated under heat stress according to GO categorization were biological process [signal transduction, response to stress, and metabolic process (of cellular amino acid, lipid, and protein)], cellular process (nucleus, cell membrane, and transcription factor complex), and molecular function (ATP binding, unfolding binding, activities of gate channels, and kinases). DEGs were mainly concentrated in the nucleus and cell membrane and were involved in protein folding, response to stress, and some signaling pathways.
The KEGG enrichment analysis of these DEGs was shown in Figure 5. In a KEGG enrichment analysis of 277 DEGs, 79 genes were annotated. KEGG results revealed the main cellular processes (endocytosis and apoptosis), environmental information processes (Fox O signaling pathway, MAPK signaling pathway, Wnt signaling pathway, and mTOR signaling pathway), genetic information processing (protein processing in endoplasmic reticulum, spliceosome, and Ubiquitin mediated proteolysis), metabolism (tyrosine, amino sugar and nucleotide, fatty acid, and amino acid), and organismal systems (longevity regulating pathway) enriched in A. mellifera under heat stress. KEGG enrichment analysis showed that DEGs were significantly enriched in genetic information processing and metabolism, among which protein processing in endoplasmic reticulum and longevity regulating pathway were significantly enriched pathways.

PPI analysis
The PPI analysis was conducted by using String database and visualized using Cytoscape. The result revealed that 18 proteins were involved in heat stress in PPI analysis (Figure 6). 60S Ribosomal protein L7 (ENSBTAG00000020139), Zinc finger proteins (ENSBTAG00000032003, ENSBTAG 00000025246, ENSBTAG00000039523), Heat shock 70 kDa protein (ENSBTAG 00000025441), Ankyrin repeat domain-containing protein (ENSBTAG00000004912), and serine/threonine-protein kinase (ENSBT AG00000004514) had strong interactions. Moreover, the PPI network showed that HSPs form a complex molecular network with other proteins, which promoted protein folding and protected cellular proteins from heat damage.

Validation of RNA-seq results
The results of the RNA-seq data showed that hundreds of genes were differentially expressed after high temperature treatment. We randomly selected 16 DEGs for q RT-PCR analysis to verify the results of RNA-seq identification. The results showed that the expression trend detected by q RT-PCR was consistent with the RNA-seq (Table 3). This ensured that the RNA-Seq results were quite reliable for the identification of DEGs.

Expression patterns of DEGs
The expression levels of 9 heat-related genes in different tissues of A. mellifera were evaluated (Figures 7A-I). The 9 DEGs were L (2) elf genes (GB45910, GB45912, GB45913, GB47475), HSP   genes (GB46774, GB50609, GB42297, GB40976), and BAG gene (GB54219). The expression of the DEGs was detected in the head, thorax, abdomen, legs, and wings, and the expression levels of most genes differed among tissues (P < 0.05). The four L (2) efl genes were similarly expressed in different tissues, all being expressed in the thorax (P < 0.05), but with no significant expression in the remaining tissues (P > 0.05) (Figures 7A-D). The four HSP genes were most highly expressed in the thorax and had low expression in other tissues (Figures 7E-H). The BAG gene was most highly expressed in the thorax, followed by the legs, with the lowest expression in the wings. All genes were highly expressed in the thorax with similar expression patterns (P < 0.05) (Figure 7I).

Discussion
Global warming and the pollination of greenhouse crops result in honeybee exposure to heat stress. Honeybees face many adverse effects of heat stress (Medina et al., 2020;Zhao et al., 2021). Ambient temperature can greatly affect foraging activity (Al-Qarni, 2006;Blazyte-Cereskiene et al., 2010), then affect crop pollination. To further understand the molecular mechanism of honeybees exposed to heat stress, a comparative analysis of gene expression among different temperatures was performed using RNA-seq, and the expression patterns of some DEGs were detected in different tissues. This study provided a set of important pathways and DEGs associated with high temperatures.
In our study, 277 DEGs were identified between CK and HT. In all, there were significantly more upregulated DEGs (167) than downregulated DEGs (110). Of these upregulated DEGs, there were 15 HSPs, including HSP40, HSP70, HSP90, L (2) efl, and so on (Supplementary Figure 2). HSPs, as highly conserved molecular chaperones, can influence the folding and functional conformation of proteins (Purschke et al., 2017;Alqarni, 2020). Both GO and KEGG enrichment analysis were associated with the protein processes, which confirmed that heat stress could induce a rapid response by the organism to synthesize resistance proteins and protect against damage (Scheffer et al., 2022). This also provided evidence for the hypothesis that elevated temperature accelerates protein unfolding and initiates molecular mates (Day et al., 2002). Protein-protein interaction (PPI) networks of DEGs. The nodes represent proteins and the lines represent interactions. The size of node represents the degree. The color of node represents the cluster coefficient, values vary from low to high with the color gradient from green to red. The green line indicates conserved neighborhood evidence, the red line indicates fusion evidence, and the yellow-green line indicates text mining. Compared with previous studies, the HSP family, such as sHSP, HSP40, HSP70, and HSP90, also participated in honeybee heat stress responses (Liu et al., 2012;Alqarni et al., 2019). We could not find any sHSPs, however, we found 6 L (2) efl genes. Under heat stress conditions, the intracellular environment is destroyed, which leads to cell apoptosis. The apoptosis process is regulated by caspase protease, and the upregulation of HSP90, HSP70, and other stress-related genes could inhibit the activity of caspase protease (Hitomi et al., 2004). We identified increased expression of HSP genes involved in multiple processes, but a striking finding was that L (2) efl genes were dominant, indicating that the genes contributed the most to 20 d old worker bees in 45 • C. The expression of HSPs may be a potential target related to the survival and performance of honeybees in harsh climatic conditions.
In addition, these DEGs were found to be associated with the zinc finger protein family (ZFPs, 5 genes), BAGs family (3 genes), and ubiquitin family (5 genes). ZFPs can bind to biological macromolecules and participate in the regulation of resistance mechanisms in response to various biological and abiotic stress factors at the transcriptional level (Andrew et al., 2004;Giri et al., 2011;Droll et al., 2013). BAG proteins may serve as bridging molecules that recruit HSP70/HSC70 to specific target proteins (Takayama et al., 1995;Rhoads and Friedberg, 1997). Thus, the BAG family functionally regulated multiple cellular pathways, such as programmed cell death (PCD) and stress responses (Kang et al., 2006). The association between the ubiquitin family and HSP90 chaperone complex has been reported (Ehrlich et al., 2009). Ubiquitin-protein activity inhibited unfolded proteins in stressinduced cell death (Chen et al., 2015), and ubiquitin mediated proteolysis was one of the recorded KEGG enriched pathways that was exposed at high temperatures in our findings. The results suggested that HSPs, BAGs, and ubiquitin family together form a complex molecular network with ZFPs. More interestingly, HSPs were also associated with immune responses (Parcellier et al., 2001;Lanneau et al., 2010;Zhang et al., 2022). Of the upregulated DEGs, toll-like receptor protein was associated with immune responses. The KEGG pathway showed enrichment of the endocytosis pathway, which was also involved in cellular immunity (Zhu et al., 2020). The endocytosis pathway may be related to the transport of proteins, thus exerting an immune function for cellular protection. Previous transcriptomic studies have indicated the relationship between heat exposure and immunity ). In addition, 5 genes related to antioxidants were expressed at heat temperatures, suggesting that heat stress can trigger oxidative stress in bees. A. mellifera may promote the body's defense against external high temperatures through immune and antioxidant responses.
We also discovered 40 DEGs related to metabolism. There were both 26 upregulated and 14 downregulated genes. The KEGG results revealed that metabolic processes were the main pathways (Figure 5). Lipids and amino acids provide energy for organisms via the tricarboxylic acid cycle. Metabolism of energy was essential to maintain homeostasis and growth (Zhu et al., 2020). Heat stress could impair the energy reserve and metabolism of insects (Slimen et al., 2016). This result was consistent with the result of physiological experiments (Yang et al., 2009;Neal, 2015), showing that high temperatures could lead to disruption in metabolism in honeybees. Transcriptome profiles also suggested that high temperatures could cause disruption at metabolic levels Shen et al., 2021). In addition, DEGs were found to be associated with stress resistance and protein kinases. A total of 8 protein kinase-related genes were identified, of which 5 were upregulated and 3 were downregulated. Protein kinases can phosphorylate proteins by transferring phosphate groups to specific substrate proteins, usually ATP or GTP as donors and serine, threonine, and tyrosine as receptors (Krupa et al., 2004;Faucher et al., 2008). The regulation of fatty acid synthase and serine protease is part of the diapause program (Robich and Denlinger, 2005;Vatanparast and Park, 2021). The heat stress caused metabolic disorders which cause diapause in honeybees. Simultaneously, adult honeybees had an ability to adapt to heat stress (Stabentheiner et al., 2010). It may be related to the higher content of metabolic enzymes.
Signal transduction pathways could link environmental stress to the organism, thus inducing most of the above-mentioned defense responses (Kawasaki et al., 2002;Ludwig et al., 2005). When organisms were subjected to heat stress, various signal transduction pathways could be activated, such as Notch, MAPK, Fox O, mTOR, Hippo, and Wnt (Badr et al., 2018). In the current experiment, some signal transduction pathways were involved, such as the Fox O, MAPK-fly, mTOR, Wnt, and TGF-beta. In the Wnt signaling pathway, downregulated of the t-cell factor/lymphoid enhancing factor (TCF/LEF) in the nucleus potentially contribute to cell cycle arrest (Kiuchi et al., 2008), and upregulated of the siah-1 interacting protein (SIP) might lead to proteolysis (Supplementary Figure 3). In the TGF-bate singling pathway, upregulated of the Samd2/3 potentially contribute to G1 arrest and lead to angiogenesis extracellular matrix neogenesis, immunosuppression, and apoptosis induction (Supplementary Figure 4). Studies have consistently shown that both Fox O (Martins and Link, 2016) and TGF-beta (Shaw et al., 2007) signaling pathways were involved in life-cycle regulation. Together with Hippo, the Wnt, TGF-beta, and Fox O signaling pathways affected body segment formation, pigmentation, appendage development, and the development of wings in insects (Barry and Camargo, 2013), which may be a key factor in diapause. Thus, the signaling transduction pathways might be the central convergence point in the A. mellifera response to heat stress.
In this experiment, all nine heat-related genes were higher expressed in the thorax of A. mellifera than in other tissues. The thorax temperature is regulated at different levels depending on food quality and demand in the colony Kovac, 2014, 2016). The regulation of the thorax temperature plays a crucial role in honeybee flight ability (Souza-Junior et al., 2020). The thorax temperature can exceed 40 • C during flight (Kovac et al., 2018). In A. mellifera, the main mechanisms for the prevention of thorax overheating during the flight are reduced metabolic and increased evaporative when the temperatures rise from 33 to 45 • C (Kovac et al., 2018). This may be related to the high expression of heat-responsive genes in the thorax. Therefore, it was supposed that the thorax may be the most important part of the response to heat stress.

Conclusion
High temperatures can influence the physiological metabolism in honeybees, and the expression profiles of stress-related genes can also be induced in response to extreme environmental temperatures. Overall, 1,743 new genes were obtained, and 277 DEGs (167 upregulated DEGs and 110 downregulated DEGs) were identified in A. mellifera under heat stress by RNA-seq. Both GO and KEGG pathway enrichment analysis showed that these DEGs were involved in protein folding, binding, response to stress, and signal transduction pathways. DEGs related to metabolism, protein processing, immune response, and signal transduction may work together to protect against heat stress in A. mellifera. This study provided a set of important pathways and DEGs associated with high temperatures. Nine DEGs were highly expressed in the thorax of A. mellifera than in other tissues, which was supposed that the thorax may be the most important part in the response to heat stress. These results will be helpful to the function research of heat resistance-related genes.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in this article/ Supplementary material.

Author contributions
WM and XL: conceptualization. WM and BZ: methodology and revision. BZ: software, data analysis, and writing the manuscript. YJ and JL: resources. BZ and JZ: visualization. All authors have read and agreed to the published version of the manuscript.

Funding
This research was supported by the earmarked fund for the Agriculture Research System of China (CARS-44-KXJ22).