1. Introduction
As the most valuable pollinators for modern agriculture and ecological diversity (The Honeybee Genome Sequencing Consortium, 2006; Klein et al., 2007; Belsky and Joshi, 2019), the development of honeybee populations is influenced by many factors (Potts et al., 2010; Belsky and Joshi, 2019; Denlinger and Yocum, 2020), including environmental temperature (Abou-Shaara et al., 2012, 2017). Environmental temperature can affect the internal and external activities of honeybee colonies (Tautz et al., 2003; Abou-Shaara et al., 2012, 2017), such as development (Tautz et al., 2003; Groh et al., 2004; Chuda-Mickiewicz and Samborski, 2015; Medina et al., 2020), phenology (Langowska et al., 2017), learnability (Tautz et al., 2003), behavioral performance (Tautz et al., 2003; Al-Qarni, 2006; Blazyte-Cereskiene et al., 2010), and susceptibility to disease (Groh et al., 2004; Conte and Navajas, 2008; Dalmon et al., 2019). And temperature can affect the biochemical and physiological processes of honeybees, including survival, tolerance, loss of body water, and the activity of antioxidant and detoxification enzymes (Dalmon et al., 2019; Li X. Y. et al., 2019). At the same time, Honeybees have evolved a series of measures to cope with temperature stress, such as behavioral avoidance (Nicolson, 2009; Ohashi et al., 2009; Joshi and Joshi, 2010), morphological changes (Langowska et al., 2017), life histories, adjusting metabolism and nutrient supply, adjusting protein synthesis or decomposition, and increasing the content of protective substances, so as to slow down or repair the damage caused by temperature stress (Ahamed et al., 2010; Li X. Y. et al., 2019; Xu et al., 2022).
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 heat-related 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) (Ma et al., 2019), 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 RNA-seq 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.
2. Materials and methods
2.1. 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 temperature-treatment × 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).
2.2. 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.
2.2.2. 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 (Kim et al., 2015), 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.
2.2.3. 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.
2.2.4. 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.
2.2.5. 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, co-occurrence, fusion, co-expression, experiments, databases, and text mining.
2.3. 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 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® 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.
2.4. 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.
4. 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). 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 stress-induced 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 (Li H. et al., 2019). 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 (Ma et al., 2019; 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 (Stabentheiner and 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.