Transcriptomic Analyses Suggest the Adaptation of Bumblebees to High Altitudes

Simple Summary Determining the adaptive mechanisms by which bumblebees adapt to high altitudes can help us to better understand their distribution, providing a basis for the future protection and utilization of bumblebee resources. For this study, the adaptive mechanisms of two dominant bumblebee species in the northeastern Qinghai-Tibet Plateau—Bombus kashmirensis and B. waltoni—were studied through transcriptomics methods. For each species, enrichment analysis of the differentially expressed genes and gene set enrichment analysis were carried out between samples collected at different altitudes (4000 m, 4500 m, and 5000 m). The results indicate that these bumblebees tend to up-regulate energy metabolism-related genes when facing extremely high-altitude environments. Of the enriched pathways up-regulated in higher altitudes, the pentose and glucuronate interconversions pathway presented the most severe up-regulation in multiple comparisons of different altitudes for B. kashmirensis, as well as the AMPK signaling pathway, which was found to be up-regulated in both species. Notably, limited by the extreme high altitudes in this study, oxidative phosphorylation was found to be down-regulated with increasing altitude, which is uncommon in studies on bumblebee adaptation to high altitudes. Abstract Determining the adaptive mechanisms by which bumblebees adapt to high altitudes can help us to better understand their distribution, providing a basis for the future protection and utilization of bumblebee resources. For this study, the adaptive mechanisms of two dominant bumblebee species in the northeastern Qinghai-Tibet Plateau—Bombus kashmirensis and B. waltoni—were studied through transcriptomics methods. For each species, enrichment analysis of the differentially expressed genes and gene set enrichment analysis were carried out between samples collected at different altitudes (4000 m, 4500 m, and 5000 m). The results indicate that these bumblebees tend to up-regulate energy metabolism-related genes when facing extremely high-altitude environments. Of the enriched pathways up-regulated in higher altitudes, the pentose and glucuronate interconversions pathway presented the most severe up-regulation in multiple comparisons of different altitudes for B. kashmirensis, as well as the AMPK signaling pathway, which was found to be up-regulated in both species. Notably, limited by the extreme hypoxic conditions in this study, oxidative phosphorylation was found to be down-regulated with increasing altitude, which is uncommon in studies on bumblebee adaptation to high altitudes.


Introduction
As one of the most important wild pollinators in both cultivated and natural ecosystems in the northern temperate world, bumblebees (Bombus spp.) play a critical role in sites were less than 100 km, in order to minimize genetic background differences. Transcriptomic differences between samples from different altitudes were analyzed in order to explore the adaptive mechanisms of these bumblebees to high altitudes. The results indicated that these bumblebees tended to up-regulate partial pathways of the energy metabolism but that the oxidative phosphorylation may unavoidably go down in extremely high-altitude environments.

Transcriptome of Samples
Sampling sites were chosen in the northeastern QTP, as this region is next to the global biodiversity center of bumblebees and a large number of species are distributed. On the other hand, the study area is close to the plateau edge and rich in high mountains, with high average altitude and severe altitude changes over a small range, making it a good choice for the study of high-altitude adaptation. Based on a previous investigation in the QTP, two indicator species, B. kashmirensis and B. waltoni [27], were chosen as the objects of this study. In order to minimize the potential influence of genetic difference caused by great altitude changes or large geographical distances, both species were collected at three collection sites along an elevation gradient including 4000 m, 4500 m, and 5000 m ( Figure 1). The distance between sites was less than 100 km. The total altitude range of B. waltoni is about 2500~5200 m [5]. The lower limit of altitude range of B. kashmirensis is about 2100 m [28], and its upper limit of altitude range is not higher than 5400 m [29]. In this study, the upper part of both species' total altitude range was covered to study the two species' adaptation to extremely high-altitude environments. In total, seventeen workers of B. kashmirensis, with five, five, and seven at 4000 m, 4500 m, and 5000 m, respectively, and thirteen workers of B. waltoni, with five, three, and five at 4000 m, 4500 m, and 5000 m, respectively, were sampled. All of the samples were active foragers with a similar body size and were sampled from flowers on sunny days. The total RNA of each sample was extracted from the body without the abdomen, and then checked by electrophoresis in agarose gel, using a Nanodrop 2000 and Agilent 2100 after removing the abdomen of species. The qualified RNA (with mass ≥ 1 µg and concentration ≥ 35 µg/µL) was used to extract mRNA using Oligo (dT) magnetic beads (Invitrogen, Carlsbad, CA, USA). The extracted mRNA was used for synthesis of cDNA by reverse transcription and random hexamers, after treatment with fragmentation buffer (Illumina, CA, USA). Then, the ends of the cDNA were repaired using End Repair Mix (Illumina, CA, USA). During this process, a single A base and adapters were added to each end of the cDNA. Then, cDNA sequencing was carried out on an Illumina NovaSeq 6000 platform after PCR amplification. Raw data was gained from the Illumina platform by CASAVA, and has been deposited into the NCBI SRA database (Bombus kashmirensis: from SRR19913284 to SRR19913300; Bombus waltoni: from SRR19913147 to SRR19913159). Clean data were obtained from raw data by removing the adapter sequence and low-quality reads. Transcriptome assembly was accomplished by Trinity v2.8.5 (k-mer = 31) [30]. Genes were annotated based on the National Center for Biotechnology Information non-redundant protein (NCBI_NR), Swiss-Prot, Protein family (Pfam), clusters of orthologous groups of proteins (COG), Gene ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using TransDecoder v5.5.0. The gene expression levels were evaluated by RSEM v1.3.1 [31], based on the transcript per million (TPM). In the following analysis, transcriptome comparisons of the two bumblebee species between different altitudes were conducted separately, and the analyses were performed using the online Majorbio Cloud Platform (www.majorbio.com (accessed on 10 March 2022)).

Identification of Adaptation-Related GO Terms
According to the annotation results, genes involved in GO terms were used to perform the gene set enrichment analysis (GSEA) using the Majorbio Cloud Platform (https:// cloud.majorbio.com/page/tools/ (accessed on 10 March 2022)). For each species, the transcriptomic differences of 4500 m vs. 4000 m, 5000 m vs. 4500 m, and 5000 m vs. 4000 m were explored through GSEA. In each comparison, the samples from the lower altitude were set as the control group and enriched GO terms were identified at the higher elevation. The enriched GO terms with the top 20 largest absolute normalized enrichment scores (NES) and p-values (adjusted by the Benjamini-Hochberg approach) lower than 0.05 were focused in the following analysis.

Identification of Adaptation-Related GO Terms
According to the annotation results, genes involved in GO terms were used to perform the gene set enrichment analysis (GSEA) using the Majorbio Cloud Platform (https://cloud.majorbio.com/page/tools/ (accessed on 10 March 2022)). For each species, the transcriptomic differences of 4500 m vs. 4000 m, 5000 m vs. 4500 m, and 5000 m vs. 4000 m were explored through GSEA. In each comparison, the samples from the lower altitude were set as the control group and enriched GO terms were identified at the higher elevation. The enriched GO terms with the top 20 largest absolute normalized enrichment scores (NES) and p-values (adjusted by the Benjamini-Hochberg approach) lower than 0.05 were focused in the following analysis.

KEGG Enrichment of Differentially Expressed Genes
For each species, the KEGG enrichment of differentially expressed genes (DEGs) for 4500 m vs. 4000 m, 5000 m vs. 4500 m, and 5000 m vs. 4000 m was assessed, wherein we identified up-and down-regulated genes at the higher elevation. Differentially expressed genes (DEGs) in each analysis were screened using DESeq 2 v1.24.0 [32]. The genes with a fold change more than two and p-value (adjusted by the Benjamini-Hochberg approach) less than 0.05 were considered as DEGs. The KEGG enrichment of differentially up-or down-regulated genes was assessed using a script developed by Majorbio (https://www.majorbiogroup.com (accessed on 10 March 2022)). The enriched KEGG pathways with p-values (corrected by the Benjamini-Hochberg approach) lower than 0.05 were considered as key pathways related to altitude adaptation.

RNA-Seq and Transcriptomic Annotation
Total RNAs of all samples of B. kashmirensis and B. waltoni were qualified. Among

KEGG Enrichment of Differentially Expressed Genes
For each species, the KEGG enrichment of differentially expressed genes (DEGs) for 4500 m vs. 4000 m, 5000 m vs. 4500 m, and 5000 m vs. 4000 m was assessed, wherein we identified up-and down-regulated genes at the higher elevation. Differentially expressed genes (DEGs) in each analysis were screened using DESeq 2 v1.24.0 [32]. The genes with a fold change more than two and p-value (adjusted by the Benjamini-Hochberg approach) less than 0.05 were considered as DEGs. The KEGG enrichment of differentially up-or down-regulated genes was assessed using a script developed by Majorbio (https: //www.majorbiogroup.com (accessed on 10 March 2022)). The enriched KEGG pathways with p-values (corrected by the Benjamini-Hochberg approach) lower than 0.05 were considered as key pathways related to altitude adaptation.

RNA-Seq and Transcriptomic Annotation
Total RNAs of all samples of B. kashmirensis and B. waltoni were qualified. Among these RNA, one of them presented an OD260/280 value of 1.98, and the others had OD260/280 values greater than two. The Q30 values of clean reads of all samples were above 90%, and the error rates were about 0.03%. The average GC content was 42.71% for B. kashmirensis samples and 46.32% for B. waltoni. The numbers of clean bases of B. kashmirensis and B. waltoni were 6.54 G and 6.74 G, respectively. The GC contents of unigenes for these two species were 40.55% and 41.67%, respectively. The mapping ratios of each B. kashmirensis sample computed by hisat2 v2.1.0 ranged from 69.74% to 83.32%, while that for B. waltoni ranged from 74.02% to 82.28% (Table 1). In the sample labels, K and W represent samples of B. kashmirensis and B. waltoni, respectively. The second letters H, M, and L in the sample labels indicate the 5000 m, 4500 m, and 4000 m sampling altitudes respectively. Q30 denotes the number of clean reads.

The Gene Set Enrichment Analysis of Gene Expression
In the GSEA of B. kashmirensis, there were 29 enriched GO terms with adjusted p-values lower than 0.05 in samples from 4500 m compared to 4000 m. The terms with the top 20 largest absolute NES (normalized enrichment scores) were mainly enriched by down-regulated genes, except for protein lipidation and steroid dehydrogenase activity terms, which were mainly enriched by up-regulated genes. The protein lipidation term has a positive effect on effective protein processing. In samples from 5000 m, 120 and 50 enriched GO terms with adjusted p-values lower than 0.05 were discovered compared to samples from 4500 m and 4000 m, respectively. In the comparison of 5000 m vs. 4500 m, all terms with the top 20 largest absolute NES were mainly enriched by up-regulated genes while for 5000 m vs. 4000 m, only six of the terms with the top 20 largest absolute NES were mainly enriched by up-regulated genes ( Table 2).
In the GSEA of B. waltoni, there were 45 enriched GO terms with adjusted p-values lower than 0.05 in samples from 4500 m compared to 4000 m. All of the terms with the top 20 largest absolute NES were mainly enriched by down-regulated genes. In the analysis of 5000 m vs. 4500 m, no GO terms presented an adjusted p-value lower than 0.05. In the analysis of 5000 m vs. 4000 m, based on the adjusted p-values, 11 enriched GO terms were significant in samples from 5000 m compared to 4000 m. Among these enriched GO terms, only the demethylation term was mainly enriched by up-regulated genes (Table 3). The samples from the lower altitude were set as the control group in each analysis. The enriched gene ontology (GO) terms with the top 20 largest absolute normalized enrichment score (NES) are shown here. A positive NES indicates that the genes associated to this term were mainly upregulated in the higher altitude while a negative NES indicates the opposite. The terms shared in the three comparisons are in bold. The genes that contributed significantly to these terms are in Table S1. The samples from the lower altitude were set as the control group in each analysis. The enriched Gene ontology (GO) terms with the top 20 largest absolute normalized enrichment score (NES) are shown here. Columns for 5000 m vs. 4500 m are absent, as no significant GO terms were found in this comparison. The meanings of the NES values are same as in Table 2. The terms shared in the two comparisons are in bold. The genes that contributed significantly to these terms are in Table S1.
In the comparison of 4500 m vs. 4000 m, the terms of chitin binding, motile cilium and odorant binding were mainly enriched by down-regulated genes in both species. In the comparison of 5000 m vs. 4000 m, the terms of odorant binding and olfactory receptor activity were mainly enriched by up-regulated genes of B. kashmirensis but down-regulated by genes of B. waltoni.

KEGG Enrichment of DEGs
In the analysis of B. kashmirensis, 23,796 DEGs in samples from 4500 m were detected, among which 19,588 genes were up-regulated and 4208 genes were down-regulated compared to 4000 m. There were also 23,796 DEGs in samples from 5000 m, among which 563 genes were up-regulated and 23,233 genes were down-regulated compared to 4500 m. There were 10,379 DEGs in samples from 5000 m identified, among which 1669 genes were up-regulated and 8710 genes were down-regulated, compared to 4000 m. As for B. waltoni, there were 3721 DEGs in samples from 4500 m, with 2967 genes up-regulated and 754 genes down-regulated, compared to 4000 m. A total of 13,625 DEGs, including 4455 up-regulated genes and 9170 down-regulated genes, were detected in samples from 5000 m compared to 4500 m. There were 14,832 DEGs, including 10,859 up-regulated genes and 3973 down-regulated genes, recognized in samples from 5000 m compared to 4000 m.
In the pathway analysis of B. kashmirensis, for 4500 m vs. 4000 m, the Fanconi anemia pathway, a pathway related to DNA repair was enriched by up-regulated DEGs. There were no pathways enriched by down-regulated DEGs in this altitude comparison. As the altitude increased from 4500 m to 5000 m, the pentose and glucuronate interconversions pathway, which promotes carbohydrate metabolism, was enriched by up-regulated DEGs. Furthermore, the oxidative phosphorylation pathway-an aerobic glycolysis process-was enriched by down-regulated DEGs. For the 5000 m vs. 4000 m comparison, the pentose and glucuronate interconversions, gap junction, arginine and proline metabolism, and AMP-activated protein kinase (AMPK) signaling pathways were enriched by up-regulated DEGs. These pathways play roles in carbohydrate metabolism, neuroprotection in hypoxic tissues [33,34], amino acids metabolism, and cellular energy homeostasis, respectively. The oxidative phosphorylation pathway was enriched by down-regulated DEGs, as was the situation for 5000 m vs. 4500 m. As shown by the results, when the altitude increased from 4000 m to 4500 m, the DNA repair process was enhanced. When the altitude continued to increase to 5000 m, carbohydrate metabolism was enhanced while the aerobic respiration was weakened. When compared to B. kashmirensis from 4000 m, in the transcriptomes of bumblebees from 5000 m, the amino acids metabolism, carbohydrate metabolism, and cellular energy homeostasis processes were all promoted, while aerobic respiration was weakened ( Figure 2).  In the pathway analysis of B. waltoni, for 4500 m vs. 4000 m, six pathways-the pentose and glucuronate interconversions, apoptosis, gap junction, AMPK signaling, phagosome, and ascorbate and aldarate metabolism pathways-were enriched by up-regulated DEGs. The pentose and glucuronate interconversions pathway was also enriched by downregulated DEGs. A similar phenomenon was observed in the analysis of 5000 m vs. 4500 m. The pentose and glucuronate interconversions pathway was also enriched by both up-and down-regulated DEGs, all with more DEGs and to a greater degree than the 4500 m vs. 4000 m comparison. Another two pathways-the spliceosome and the ascorbate and aldarate metabolism pathways-were also detected to be enriched by down-regulated DEGs in this altitude shift. In the 5000 m vs. 4000 m comparison, the pentose and glucuronate interconversions pathway was again enriched by downregulated DEGs, while there were no other pathways enriched by up-or down-regulated DEGs in this comparison. There were some pathways that were up-regulated with increasing elevation in the results of both species, but the up-regulation of these pathways occurred in different altitude-change situations for the two species. For B. kashmirensis, the pentose and glucuronate interconversions, gap junction, and AMPK signaling pathways were up-regulated at the high altitude of 5000 m, while they were already detected to be enhanced at 4500 m for B. waltoni (Figure 3).

Discussion
Determining the mechanisms by which bumblebees adapt to high altitudes can help us to better understand their distribution, allowing for further prediction of future range shifts driven by climate change, which may serve as a basis for the protection and utilization of bumblebee resources. In this study, the northeastern QTP was selected as the study area, and two dominant and widespread bumblebee species-B. kashmirensis and B. waltoni-were selected as the study subjects. Differing from previous bumblebee transcriptome studies, here, the distances and altitude differences between collection sites were respectively controlled to no more than 100 km and 1000 m, in order to minimize the influence of genetic differences possibly caused by great geographical distances or altitude differences on the transcriptomic results.
As shown in our results, bumblebees tended to up-regulate energy metabolism-related genes in extremely high-altitude environments. Many pathways closely related to energy production were significantly up-regulated in both studied bumblebee species with increasing altitude. The pentose and glucuronate interconversions pathway, which is an important component of carbohydrate metabolism, was enriched by up-regulated DEGs in B. kashmirensis at the high altitude of 5000 m, compared to 4000 m and 4500 m (Figure 2). The AMPK signaling pathway, which has positive effects on cellular energy homeostasis by maintaining sufficient adenosine triphosphate (ATP) resources, was observed to be enhanced in both species with increasing altitude: from 4000 m to 5000 m for

Discussion
Determining the mechanisms by which bumblebees adapt to high altitudes can help us to better understand their distribution, allowing for further prediction of future range shifts driven by climate change, which may serve as a basis for the protection and utilization of bumblebee resources. In this study, the northeastern QTP was selected as the study area, and two dominant and widespread bumblebee species-B. kashmirensis and B. waltoni-were selected as the study subjects. Differing from previous bumblebee transcriptome studies, here, the distances and altitude differences between collection sites were respectively controlled to no more than 100 km and 1000 m, in order to minimize the influence of genetic differences possibly caused by great geographical distances or altitude differences on the transcriptomic results.
As shown in our results, bumblebees tended to up-regulate energy metabolism-related genes in extremely high-altitude environments. Many pathways closely related to energy production were significantly up-regulated in both studied bumblebee species with increasing altitude. The pentose and glucuronate interconversions pathway, which is an important component of carbohydrate metabolism, was enriched by up-regulated DEGs in B. kashmirensis at the high altitude of 5000 m, compared to 4000 m and 4500 m (Figure 2). The AMPK signaling pathway, which has positive effects on cellular energy homeostasis by maintaining sufficient adenosine triphosphate (ATP) resources, was observed to be enhanced in both species with increasing altitude: from 4000 m to 5000 m for B. kashmirensis ( Figure 2) and from 4000 m to 4500 m for B. waltoni (Figure 3). Additionally, enrichment by up-regulated DEGs from the arginine and proline metabolism pathways-both of which belong to the amino acids metabolism category-were detected in B. kashmirensis in the 4500 m vs. 4000 m ( Figure 2) and 5000 m vs. 4000 m (Figure 2) comparisons, respectively. As implied by the above results, the energy production of bumblebees was obviously strengthened with increasing altitude. The reason why energy production should be enhanced may be related to the increasing needs for heat and flight power with increasing altitude. From the perspective of heat, as the ambient temperature will consistently decrease with the rising elevation, increased convective heat loss is inevitable when bumblebees live in colder environments at higher altitudes. Thus, the bumblebees must produce more heat through metabolism and utilization of ATP [27,30,35]. Enhancement of energy substance metabolism and the AMPK signaling pathway can meet the increasing heat demand of bumblebees with increasing altitude. From the perspective of flight power, the forces produced by insects when flapping their wings are directly proportional to air density [22]. The low air density at higher altitudes is a great disadvantage of flight [36,37]. At higher altitudes, insects have to use more energy to flap their wings more frequently, in order to the aerodynamic forces required to offset their body weight [38,39]. Thus, bumblebees up-regulate energy metabolism-related genes as the altitude increases.
In addition to the low temperature and low air density mentioned above, bumblebees at high altitudes are still faced with another great challenge: hypoxia. Liu et al. [25] have considered high-altitude adaptation under a tremendous altitude difference of over 4000 m, and reported that the high-altitude bumblebees had to enhance their aerobic respiration (e.g., through the oxidative phosphorylation pathway), in order to provide enough energy to ensure normal flight, compared to the low-altitude bumblebees. In this study, we seemed to observe a completely different phenomenon. Compared to 4000 m and 4500 m, oxidative phosphorylation was enriched by down-regulated DEGs in B. kashmirensis at the high altitude of 5000 m (Figure 2). The paradox regarding the adaptive change of oxidative phosphorylation to the hypoxia may be attributed to the difference in the level of altitude shift considered in the two studies and the limitation of the adaptation of respiration processes and oxidative phosphorylation. In responding to hypoxia, to some extent, insects can alter their respiratory system [36]. They can increase the number of terminal tracheal branches [40,41] and the diameter of primary tracheae, in order to compensate for the hypoxic environment [42][43][44]; however, the reasonable allocation of muscle fibers, mitochondria, tracheae, and sarcoplasmic reticulum in the thorax volume limits such an increased respiratory system [45]. Thus, to a certain degree, insects can adapt to the hypoxic environment by strengthening their respiratory system and oxidative phosphorylation, in line with the result of Liu et al. [25]. However, when the increased respiratory system still fails to provide enough oxygen under extreme hypoxia at extremely high altitudes, the adaptation mechanism may change and the aerobic respiration is abated to avoid redundancy as much as possible. Here, as our study was aimed at determining the adaptation to extremely high altitudes in the alpine region in the northeastern QTP, all three sampling sites were designed at altitudes no less than 4000 m, which may be close to the limit of respiratory system and oxidative phosphorylation adaptation [46]. Oxidative phosphorylation was enriched by down-regulated DEGs at the high altitude of 5000 m.
As indicated by the above differences in the results of the different comparisons, there were some gradual changes as the altitude increased. For example, in the analysis of B. kashmirensis, the up-regulation of arginine and proline metabolism and the up-regulation of the AMPK signaling pathways were not significant in the comparisons of 4500 m vs. 4000 m and 5000 m vs. 4500 m but were significant in the comparison of 5000 m vs. 4000 m (Figure 2). This implied the gradual changes in the two processes when the elevation changes from 4000 m to 5000 m of the two processes. Meanwhile, some processes demonstrated transcriptomic changes to quite different degrees under different situations of increased altitude. This was also well indicated by the different comparisons considered in this study; for example, in the analysis of B. kashmirensis, the pentose and glucuronate interconversions pathway was greatly up-regulated in the comparison of 5000 m vs. 4500 m but not in the comparison of 4500 m vs. 4000 m (Figure 2). The down-regulation of oxidative phosphorylation demonstrated a similar trend (Figure 2). This suggests that there are more greater changes in these processes in B. kashmirensis for them to respond to the more extreme altitude increasing. In the analysis of B. waltoni, the AMPK signaling pathway was up-regulated in the comparison 4500 m vs. 4000 m but not regulated in comparisons 5000 m vs. 4500 m and 5000 m vs. 4000 m (Figure 3). This might be relative to the adaptation limit of B. waltoni. Although B. waltoni could survive at 5000 m, this extreme high-altitude environment might pose a great challenge and even exceed its adaptation limit. Energy regulations could not meet its original energy-intensive activity. Then, B. waltoni has to limit its own activity, with the AMPK signaling pathway being down-regulated indistinctively in the comparison of 5000 m vs. 4500 m.
It can also be inferred that the mechanisms of extremely high-altitude adaptation for the two bumblebee species (i.e., B. kashmirensis and B. waltoni) differed. Compared to cospecies at 4000 m, B. kashmirensis at 5000 m showed up-regulated pentose and glucuronate interconversions pathway, arginine and proline metabolism pathway, and AMPK signaling pathway ( Figure 2) while B. waltoni at 5000 m showed a down-regulated pentose and glucuronate interconversions pathway and no up-regulated energy metabolism (Figure 3). This may be due to the difference between the two species, regarding their need for energy in response to altitude changes. Aside from the difference regarding the regulation of energy metabolism, there were also differences between the adaptation mechanisms of these species on other aspects. Peroxidation stress can be seen as a common challenge caused by the hypoxia at high altitudes, and DNA repair and apoptosis have been interpreted as two physiological processes related to the response to peroxidation stress [47]. When the altitude increased from 4000 m to 4500 m, B. kashmirensis enhanced its DNA repair function ( Figure 2) while B. waltoni strengthened the apoptosis process ( Figure 3). Furthermore, for B. kashmirensis, the gap junction and AMPK signaling pathways were up-regulated at the high altitude of 5000 m ( Figure 2) while they were already detected to be enhanced at 4500 m in B. waltoni (Figure 3). These results may provide a clue about the different adaptations of two bumblebee species to extremely high altitudes. The different adaptation might be related to genetic background. Bumblebees belonging to different subgenera had different genetic strategies for adapting to high altitudes [23]. The different genetic strategies made them have different transcriptional regulation processes to respond to high altitudes [25]. The genetic differences between the two species considered here resulted in the inherent differences in their transcriptional regulation. The extreme high altitude further resulted in their transcriptional regulation being more targeted in order to meet their own needs more effectively. Therefore, there were many differences regarding the adaptation demonstrated via transcriptional regulation between the two species. However, to obtain a definite explanation regarding these phenomena, larger sample sizes and more comprehensive sampling covering the whole altitude range of the two species, as well as investigations of their distribution patterns, are necessary.

Conclusions
We showed, according to the results of our study, that the considered bumblebees tend to up-regulate energy metabolism-related genes when adapting to extremely high-altitude environments. Of the enriched pathways up-regulated in higher altitudes, the pentose and glucuronate interconversions pathway presented the most severe up-regulation in multiple comparisons at different altitudes for B. kashmirensis, while the AMPK signaling pathway was up-regulated in both species. However, due to the extreme hypoxia in this study, oxidative phosphorylation was detected to be down-regulated as the altitude kept increasing, which is uncommon in studies on bumblebee adaptation to high altitudes.