Cuticular hydrocarbons of alpine bumble bees (Hymenoptera: Bombus) are species-specific, but show little evidence of elevation-related climate adaptation

Alpine bumble bees are the most important pollinators in temperate mountain ecosystems. Although they are used to encounter small-scale successions of very different climates in the mountains, many species respond sensitively to climatic changes, reflected in spatial range shifts and declining populations worldwide. Cuticular hydrocarbons (CHCs) mediate climate adaptation in some insects. However, whether they predict the elevational niche of bumble bees or their responses to climatic changes remains poorly understood. Here, we used three different approaches to study the role of bumble bees’ CHCs in the context of climate adaptation: using a 1,300 m elevational gradient, we first investigated whether the overall composition of CHCs, and two potentially climate-associated chemical traits (proportion of saturated components, mean chain length) on the cuticle of six bumble bee species were linked to the species’ elevational niches. We then analyzed intraspecific variation in CHCs of Bombus pascuorum along the elevational gradient and tested whether these traits respond to temperature. Finally, we used a field translocation experiment to test whether CHCs of Bombus lucorum workers change, when translocated from the foothill of a cool and wet mountain region to (a) higher elevations, and (b) a warm and dry region. Overall, the six species showed distinctive, species-specific CHC profiles. We found inter- and intraspecific variation in the composition of CHCs and in chemical traits along the elevational gradient, but no link to the elevational distribution of species and individuals. According to our expectations, bumble bees translocated to a warm and dry region tended to express longer CHC chains than bumble bees translocated to cool and wet foothills, which could reflect an acclimatization to regional climate. However, chain lengths did not further decrease systematically along the elevational gradient, suggesting that other factors than temperature also shape chain lengths in CHC profiles. We conclude that in alpine bumble bees, CHC profiles and traits respond at best secondarily to the climate conditions tested in this study. While the functional role of species-specific CHC profiles in bumble bees remains elusive, limited plasticity in this trait could restrict species’ ability to adapt to climatic changes.


Introduction
Mountains are hotspots of biodiversity with many endemic and cold-adapted species (Rahbek et al., 2019;Trew and Maclean, 2021). Yet climate change affects biodiversity by changing spatial ranges, phenology and species interaction (Bellard et al., 2012;Halsch et al., 2021). Since the temperature increase associated with climate change is likely to be more pronounced at higher elevations than at lower elevations (Pepin et al., 2015), our urgent concern should be to understand the responses of these particularly threatened ecosystems, which can also serve as early warning system.
Alpine bumble bees (Apidae: Bombus) are the most efficient and widespread pollinators in temperate mountains (Bingham and Orthner, 1998;Gorenflo et al., 2017). Within the genus Bombus, species differ substantially in their preferred elevational niche (Rasmont et al., 2015;Minachilis et al., 2020;Sponsler et al., 2022a), with few species being restricted to the cool and exposed conditions above the tree line. Within a species, populations and individuals can face substantial differences in their climatic environments along mountain slopes or during foraging flights under the typically fluctuating weather conditions of the mountains (Sponsler et al., 2022a), occasionally, when the cloud cover shifts abruptly, within minutes. Despite their broad temperature tolerance (Peters et al., 2016), many bumble bee species are sensitive to climate change. Many species are globally declining Arbetman et al., 2017) and shift their range to higher elevations or toward the poles to track their preferred temperature niche and to avoid heat and desiccation stress (Kerr et al., 2015;Pyke et al., 2016;Marshall et al., 2020). The magnitude of such shifts varies greatly among species and is likely associated with species capacity to withstand new environmental conditions (Kerr et al., 2015;Pyke et al., 2016;Marshall et al., 2020;Maihoff et al., 2022b). One factor known to influence adaptation to climate change is desiccation resistance. Although studies focusing on gene expression and heat resistance sugguest that the ability to tolerate desiccation, caused by high temperatures and lower percipitation, is relevant for determining bumble bees' response to climate change (Jackson et al., 2020;Maebe et al., 2021;Martinet et al., 2021) and, consequently, the maintenance of their pollination service, mechanisms that determine desiccation resistance can be versatile. Desiccation resistance in insects depend on the ability of insects cuticle to prevent water loss (Chown et al., 2011). Thus, differences in physiological traits modulating water loss are likely to contribute to their uneven ability to cope with environmental conditions along the elevational gradient (Williams et al., 2010;Aguirre-Gutiérrez et al., 2016;Wong et al., 2019;Stemkovski et al., 2020;Maebe et al., 2021). However, such traits have not yet been identified in bumble bees.
The cuticular hydrocarbons (CHCs) composition on the insects' cuticle is such a trait. CHCs provide waterproofing, and comprise a complex mixture of n-alkanes, methyl-branched alkanes and unsaturated hydrocarbons that cover the cuticle (Gibbs et al., 1991;Gibbs, 1998;Blomquist and Bagnères, 2010). Two features are particularly relevant in the context of waterproofing (Gibbs and Pomonis, 1995): the proportion of saturated components and the chain length of hydrocarbons. Saturated hydrocarbons (n-alkanes and monomethyl-branched alkanes) aggregate more tightly than cis-configurated unsaturated hydrocarbons (n-alkenes) due to increased van-der Waals forces. Thus, a higher proportion of saturated components in the CHC profile enhances waterproofing (Gibbs and Pomonis, 1995). Aggregating forces between hydrocarbons also increase with chain length; CHC profiles composed of on average longer hydrocarbons thus provide a better protection against desiccation than profiles with shorter CHCs (Menzel et al., 2017b).
Various studies on insects confirm that the composition of CHC profiles can reflect climate adaptation on an inter-and intraspecific level (Rajpurohit et al., 2017;Menzel et al., 2017a;Michelutti et al., 2018;Sprenger et al., 2019;Mayr et al., 2021). For example ant species from habitats with high rainfall produce various alkenes, alkadienes and methylbranched alkenes, i.e., substance classes with reduced protection against desiccation stress (Menzel et al., 2017a), while Drosophila populations from warmer regions produce longer CHC chains than population from colder regions (Rajpurohit et al., 2017). Thus, it can be hypothesized that bumble bee species that differ in their preferred environmental niche have CHC profiles that reflect the degree of desiccation stress. Importantly, in some insect species CHC profiles alter plastically in a short-term under changing climatic conditions Sprenger et al., 2018). This intraspecific acclimatization of CHC changes can be elicited even within hours (Stinziano et al., 2015) and can effectively reduce heat stress . Thus, the capacity for a short-term acclimatization response, alongside long-term adaptation, may be crucial in determining the climatic range along the elevational gradient of a species and for its survival under changing climatic conditions.
In addition to the function of water balance regulation, CHCs are involved in insect communication. In social insects communication is crucial for colony maintenance. CHCs are important for nestmate recognition and task signalization, in, e.g., wasps (Polistes) honey bees (Apis), stingless bees and ants (Akino et al., 2004;Dani et al., 2005;Nunes et al., 2008;Leonhardt et al., 2016;Maihoff et al., 2022a) and also in caste recognition and signalizing of health status (Leonhardt et al., 2016;Beani et al., 2019). Thus, CHCs that are involved in species communication and recognition of nest affiliations should be under strong pressure to be maintained, as any deviations may have negative effects on fitness and colony survival. Even though the role of CHCs in inter-and intraspecific bumble bee communication remains unclear, the potential dual function of methyl-alkanes, alkanes and alkenes in communication and waterproofing could potentially limit species capacity to adapt or acclimatize to new climatic conditions (Dani et al., 2005;Colazza et al., 2007;Lacey et al., 2008;Sprenger et al., 2019;Awater-Salendo et al., 2020). Furthermore, like many other functional traits, CHCs may be constrained by phylogeny, which should be considered when analyzing their variation within a clade (De Oliveira et al., 2011;Flynn et al., 2011;Kellermann et al., 2012;Kather and Martin, 2015;Menzel et al., 2017a). In general, less is known about CHC profile variation within the clade of bumble bees, but in the context of studies addressing climate change, genetic variation at cuticle formation suggests a likely role in climatic adaptation (Jackson et al., 2020;Straub et al., 2022).
In this study, we used three different approaches to study the role of CHCs in climatic adaptation across and within bumble bee species. We hypothesized that the CHC profiles of bumble bees exposed to higher desiccation stress at higher temperatures have longer chains and a higher proportion of saturated components. A 1.2 km elevational gradient serves as a natural model system that provides variation in temperatures and associated relative humidity patterns also expected under climate change (up to +5°C by the end of the century; Körner, 2007;IPCC et al., 2018) as temperature declines linearly with elevation (∼0.5°C per 100 m increase in elevation; Körner, 2007). Furthermore, the tree line, i.e., the transition from forest to shrub or grassland, represents a critical threshold, as outside the protection of the tree canopy the various abiotic conditions, to which bumble bees are Frontiers in Ecology and Evolution 03 frontiersin.org exposed, become more extreme (Beals, 1969;Slatyer and Noble, 1992). We first compared the CHC profiles of six bumble bee species that differed in their average elevational distribution and hypothesized that CHC composition and assumed to be climate-associated chemical traits (i.e., mean chain length and proportion of saturated components) correlate with the preferred climate niche of the respective species (Mayr et al., 2021). For simplicity, we refer to these traits climate-associated chemical traits in the following. We also consider phylogenetic constraints, assuming that traits are more similarly expressed in closely related species than in distant ones (Flynn et al., 2011). In a second approach, we investigated intraspecific variation in the CHC profile of Bombus pascuorum (Scopoli, 1763)-a species, which was previously restricted to lower elevation and which has expanded its range to higher elevation during recent climate change (Marshall et al., 2020;Maihoff et al., 2022b). Specifically, we hypothesized that variation in climate-associated chemical traits on the cuticle correlates systematically with environmental change along the elevational gradient and becomes most pronounced above the tree line where abiotic conditions are more extreme.
The third approach entails a field translocation experiment, to test the acclimatization capacity of CHCs in Bombus lucorum (Linnaeus, 1761), a species known to prefer cold and humid forest habitats (Geue and Thomassen, 2020). Field translocation experiments are powerful tools to assess the extent of species' acclimatization and genetic determination, as they elucidate the factors that limit species distributions and have the potential to predict species range shifts in a changing climate (Nooten and Hughes, 2017). We transferred young bumble bee colonies derived from queens, that were collected in the foothills of a comparably cold and wet mountainous region (Berchtesgadener Land, Germany) to (a) higher elevations within the mountainous region and (b) a warmer and drier region in Bavaria's lowlands (Würzburg, Germany). We hypothesized that after 4 weeks differences in CHC profiles and climate-associated chemical traits are indicative for new nesting sites in the way that in the warmer and dryer regions CHC profiles are characterized by a high proportion of saturated components and on average longer mean chain lengths.

Materials and methods
Study area and bumble bee sampling Approach 1: CHC differences between species occupying different elevational niches To test whether differences in CHC profiles between species can be explained by species' elevational niche, bumble bee workers (foraging at flowers) of six species were collected within the National Park Berchtesgaden and its close vicinity (Lat: 47.5477, Lon: 12.9247) for one month (29.07.-31.08.) in two years (2019 and 2020). We selected six coexisting species in the study region to represent differences in elevational distribution. Bombus mucidus Gerstäcker (1869) and Bombus monticola Smith (1849) mainly occur above the tree line at the highest sites of our study region, Bombus soroeensis Fabricius (1777), and Bombus wurflenii Radoszkowski (1859) across the entire studied elevational gradient, and B. pascuorum and B. lucorum at lower elevations below the tree line on average at an elevation of 1,200 m.a.s.l.
[monitored in Maihoff et al., 2022b; see Figure 1 and further information in Supplementary material S1]. The National Park is located within the limestone Alps in southern east Germany-a region characterized by coniferous forest, alpine meadows, and mountain pastures (Konnert and Siegrist, 2000). The tree line in the study area is at an elevation of about 1,500 m.a.s.l. (Köstler and Mayer, 1970). We collected individuals from sites (60 m × 60 m) representative of the species' range and also aimed for different species at the same sites. Therefore, mainly individuals above 1,300 m were collected, except for B. lucorum and B. pascuorum, where we included individuals from lower elevations according to their distribution. In total we collected 65 individuals [B. soroeensis (n = 9), B. mucidus (n = 9), B. lucorum (n = 12), B. monticola (n = 12), B. wurflenii (n = 11), and B. pascuorum (n = 12)] from 18 sites covering a gradient from 641 to 2,114 m.a.s.l. (Figure 1). Note that the species were selected not only according to their preferred distribution, but also according to their relatedness, which does not correlate with the distribution (see Figure 2A for phylogeny).
Approach 2: Intraspecific CHC variation of Bombus pascuorum along an elevational gradient For analyzing intraspecific variation in CHCs, we collected B. pascuorum workers foraging at flowers (n = 58) from 12 sites (60 m × 60 m) across an elevational gradient from 641 to 2,032 m.a.s.l. within the National Park Berchtesgaden and its close vicinity (Lat: 47.5477, Lon: 12.9247). We selected B. pascuorum because of its high abundance in the study region, and its characteristics as a species expanding its range toward higher elevations under recent climate change .

Approach 3: Intraspecific CHC variation in translocated workers of Bombus lucorum
To test whether bumble bees can adjust their CHC profile to new climatic conditions, we translocated young self-reared colonies, originally from cold and wet climate to a warm and dry climate and to distinct elevations. For this, we collected queens, emerging from diapause, of B. lucorum in the alpine study region in spring and reared them in the laboratory under constant climatic conditions [30°C, 60% humidity following Requier et al. (2020)]. Rearing continued until at least 11 workers per queen were hatched. By the end of May (28th of May +/− 3 days), we settled 13 young colonies (queens + workers) in nest boxes and placed them in two climate regions [warm and dry in the Franconian lowlands (264 m.a.s.l.) around Würzburg (n = 5) and the cold and rather wet mountain region around Berchtesgaden, where queens originated from (n = 8); Figure 1B]. The regions differ strongly in their multi-annual mean temperature ( Figure 1B) and precipitation (Supplementary material S2): The warm and dry region has an average annual temperature of 9°C and annual precipitation of about 650 mm, while the alpine region of origin has an average annual temperature of 7°C and > 1,500 mm annual precipitation. Within the warm and dry lowlands colonies were placed at one site. Within the alpine region colonies were placed at three sites [mountain foothill: 752 m.a.s.l. (n = 3); mountain mid: 1,100 m.a.s.l. (n = 3); and mountain high: 1,933 m.a.s.l. (n = 2) for temperature differences per site see Figure 1]. Because temperature increase and precipitation decrease are changing in the same direction within this approach (Supplementary material S2) we consider parallel changing conditions between regions. After at least 4 weeks in the field, we collected 3 individuals per colony (individuals leaving and entering the nest entrance). We cannot prove that the collected individuals really hatched in the field, because bumble bee workers' lifespan averages from 22 to 69 days and depends on the species and environmental conditions (Goldblatt and Fell, 1987;Smeets and Duchateau, 2003;Kelemen et al., 2019). Thus, we here speak of Frontiers in Ecology and Evolution 04 frontiersin.org individuals exposed to different conditions and not of individuals that have developed under different climatic conditions. Collected individuals were flash frozen on dry ice in the field and stored at −20°C for further analyses.

Temperature analyses
In the study region Berchtesgaden, temperature data were modeled at the site level from daily temperature data averaged from neighboring climate stations (18 in total). The statistical details of temperature predictions are provided in Supplementary material S3. Temperature data were used to calculate short term acclimatization temperature (mean temperature one week before sampling = Tacclim) and mean annual temperature (MAT). In the warm and dry study region near Würzburg, temperature data originated from the closest German Meteorological Service (Deutscher Wetterdienst) climate station, which was at the same elevation and 1 km away from the colonies. Within the translocation experiment we calculated mean exposure temperature (i.e., mean temperature experienced by Study design. (A) Study region and sites in Berchtesgaden (Bavaria, Germany). Elevation level is indicated in shades (at 250 m intervals), with darker shades signaling increasing elevation. The National Park border is represented with a green line. Points show study sites. Red sites were used for bumble bee collection in approach 1 and 2, while yellow, green, and blue sites represent sites which were additionally used for the translocation approach (approach 3). (B) Schematic description of the colony translocation (approach 3). Queens collected from the valleys in spring in our study region were reared in the lab and then settled with their first offspring in nest boxes. Numbers indicate new colony locations. In the cold and wet region of origin: mountain foothills (1), mountain mid (2), mountain high (3) (see also Panel A, where mountain sites are shown in respective colors), and in the warm and dry region (4), where the mountain sites are shown according to their color in Panel A. The table lists the temperature regimes to which the colonies were exposed either over the course of the translocation (=Texp) or within the week before sampling along the elevational gradient (=Tacclim). Color code in the map refers to multi annual means of air temperature conditions. Maps were produced in QGIS. Data was obtained from https://www.lfu.bayern.de/umweltdaten, https:// search.earthdata.nasa.gov, and http://www.dwd.de.
Frontiers in Ecology and Evolution 05 frontiersin.org colonies from the establishment in the field to the collection of the workers = Texp) and the short-term acclimatization temperature one week before sampling (=Tacclim). Coordinates for all climate stations are given in the Supplementary material S3.

Chemical analyses
CHC profiles were extracted from legs and wings of all workers, as these body parts were shown to have the highest volume to surface ratio and best reflect climate-associated chemical traits in ants (Sprenger et al., 2021) while being representative for a species and revealing individual differences at the same time (Young et al., 2000;Wang et al., 2016;Mayr et al., 2021;Sprenger et al., 2021). Extraction was performed by immersing pooled body parts in n-hexane for 10 min per individual. The extracts were concentrated under gentle CO 2 stream to approximately 20 μL and transferred to a micro insert. CHC extracts were analyzed with an Agilent 6890 gas chromatograph coupled with an Agilent 5975 Mass Selective Detector (GC-MS, Agilent, Waldbronn, Germany): The GC (split/splitless injector in splitless mode for 1 min, injected volume 1 μL at 300°C) was equipped with a DB-5 Fused Silica capillary column (30 m × 0.25 mm ID, df = 0.25 μm; J&W Scientific, Folsom, United States). Helium served as carrier gas at a constant flow of 1 mL/min. The following temperature program was used: Start temperature 60°C, temperature increase by 5°C per min up to 300°C, isotherm at 300°C for 10 min. The electron ionization mass spectra (EI-MS) were acquired at an ionization voltage of 70 eV (source temperature: 230°C). Chromatograms and mass spectra were recorded and quantified via integrated peak areas with the software HP Enhanced ChemStation G1701AA (version A.03.00; Hewlett Packard). CHC compounds were identified by the compound specific retention indices and their detected diagnostic ions (Carlson et al., 1998). Alongside CHC samples, we run an analytical alkane standard (C 8 -C 20 and C 21 -C 40 ; Sygma Aldrich) for the calculation of the retention index and to check for the sensitivity of the GCMS. Double-bond position in alkenes were identified by DMDSderivatization following Dunkelblum et al. (1985). The GC settings were the same as described before, but oven temperature increased by 5°C per min up to 325°C and isotherm at 325°C for 10 min. Frontiers in Ecology and Evolution 06 frontiersin.org

Statistical analyses
All analyses were performed using the software R version 4.1.3. (R Core Team, 2021).
Approach 1: Differences in CHC composition and climate-associated chemical traits between species occupying different elevational niches In all six bumble bee species, we compared the relative abundances of compounds in the CHC profile of workers. Only CHC compounds which contributed to at least 0.1% of the total abundance of compounds in the CHC profile were analyzed. CHC profiles were assessed by non-metric multidimensional scaling (NMDS), a two-dimensional ordination method to visualize similarity, and agglomerative hierarchical cluster analysis ( Figure 2B). We used the "plot_ordination" function from the "phyloseq" package and, respectively, the "hclust" function from the "vegan" package and "as.dendrogram" function from the "dendextend" package. Dissimilarities between worker profiles were calculated using Bray-Curtis distances. We assessed CHC profile composition differences between bumble bee species by using permutational multivariate analysis of variance (PERMANOVA) in the packages vegan (Oksanen et al., 2020) and pairwiseAdonis (Martinez Arbizu, 2020). Permutations were set on 10,000 or 999, respectively.
Overall differences in CHC profiles between species, however, do not inform per se about the protection capacity of a cuticle against desiccation. We therefore calculated more informative climate-associated chemical traits, i.e., (a) the proportion of saturated components, and (b) the abundance weighted mean chain length. Both values are predicted to increase with evaporation protection capacity (Gibbs and Pomonis, 1995). We used a generalized linear mixed model (function "glmmTMB" Brooks et al., 2017) with mean chain length and the proportion of saturated components, respectively, as response variables and species as the predictive variable. Since individuals were collected in two different years, we included year as a random effect in the model. We run a gaussian-distributed model in each case. Models were determined and checked with the DHARMa package. Model selection (null vs. full model) was performed using the dredge function from the MuMIN package (Barton, 2020) based on Akaike information criterion corrected for small sample sizes (AICc) where the lowest AICc relative to other models (here null model) indicates the preferred model (ΔAICc > 2 indicates statistically relevant differences; Burnham and Anderson, 2004). Post hoc comparison among species were performed with the emmeans function with a Tukey adjustment.
We attempted to elucidate the mechanisms that underly CHC responses to temperature, by testing and comparing CHC responses to two different temperature regimes (MAT vs. Tacclim): We assume that CHC traits that are better explained by MAT than by Tacclim, indicate species adaptation to the respective temperature niche along the elevational gradient. In contrast, we assume that CHC traits that are better explained by Tacclim, suggest a rather spontaneous and likely reversible response of CHCs to the current temperature conditions. We addressed these different temperature effects with separate generalized linear mixed models, as MAT and Tacclim were moderately correlated (corr = 0.76; see Supplementary material S4). In linear mixed models we included temperature as fixed effect and species nested in temperature as random effect (because certain temperatures exist only for certain species). We run a gaussian-distributed model in each case.
To address phylogenetic constraints within CHC profile compositions we performed a Mantel test, which compares phylogenetic distances of species with chemical distances (=calculated Bray-Curtis dissimilarities) between averaged species profiles following Buellesbach et al. (2013) and Martin et al. (2013). We obtained phylogenetic distances from the comprehensive phylogeny provided by Cameron et al. (2007). Furthermore, we analyzed phylogenetic signals (Morans'I) in the climate-associated chemical traits by phylogenetic correlograms (Gittleman and Kot, 1990), which measure the correlation between phylogenetic distances and trait distances, using the R package phylosignal (Keck et al., 2016). Higher correlation indicates that CHC variation is stronger driven by phylogeny. As an initial tree we used again the comprehensive phylogeny of Cameron et al. (2007), but this time we assigned individuals to species (Figure 2A). Hereby we defined the distance between individuals of a given species to 0.02. By doing so, we standardized the unknown genetic distances between individuals to conduct our analyses at the individual level. Our analysis should therefore only be considered as an approximation for the relatedness of individuals.

Approach 2: Intraspecific CHC variation of Bombus pascuorum along an elevational gradient
To test for intraspecific CHC differences in B. pascuorum along the elevational gradient, we first assessed profile similarity by NMDS and performed a PERMANOVA to test for an elevational effect on CHC profile composition. We then conducted a generalized linear mixed model with proportion of saturated components and mean chain length as response variables and elevation as predictor variable. We included site as random effect, because individuals sampled on the same site might per se show a higher CHC similarity than individuals sampled on different sites (e.g., due to relatedness or response to other, non-considered site factors). We run a betadistributed model for the proportion of saturated components and a gaussian-distributed model for the mean chain length, both were determined and checked with the DHARMa package. Due to high correlation between elevation, MAT and Tacclim (Supplementary material S5) we cannot separate the mechanisms underlying temperature responses here. We thus analyzed elevational patterns, which can be interpreted as response to temperature changes. We used a likelihood ratio test (LRT) to test whether the model explains more than a random change. Test was performed with the "anova" function from the "stats" package.
Further, we assessed differences in the coefficient of variation (CV) of the two climate-associated chemical traits along the elevational gradient. The CV, which equals the standard deviation of the respective trait divided by the trait mean, is assumed to decrease with increasing elevation, as increased environmental filtering processes, especially above the tree line may constrain phenotypic variability. We tested for an elevational effect on the coefficient of variation for each site with generalized additive models (GAMs) allowing the flexible detection of linear and non-linear relationships (Wood, 2011). The basis dimension of smoothing term was set to five to have enough degrees of freedom to represent the underlying pattern well, but small enough to maintain reasonable computational efficiency. GAMs were computed with the ´gam´ function from the mgcv package (Wood, 2006).

Bombus lucorum after translocation
To test for a translocation effect on CHC profile composition in B. lucorum, we assessed profile similarity using NMDS and PERMANOVA. We further, tested for differences in climate-associated Frontiers in Ecology and Evolution 07 frontiersin.org chemical traits in translocated workers with generalized linear mixed models. We included colony identity as a random effect in the generalized mixed effect models. See Supplementary material S6 for NMDS displaying CHC-profile similarity from individuals out of the same colonies. We run a gaussian-distributed model in each case. Models were determined and checked with the DHARMa package. To elucidate the prevailing mechanisms behind responses, we again tested two different temperature measures against each other: the average exposed temperature throughout the study period (=Texp), whose explanatory power should be high if workers respond with a change in CHC composition to temperature exposure over several weeks vs. Tacclim, i.e., mean temperature one week before sampling, which suggest a short-time response. Due to high correlation of both Texp and Tacclim with colony location (warm and dry region; mountain foothill, mountain mid, mountain high; see Supplementary material S6) we cannot separate the mechanisms underlying temperature responses here and only present the results of the model with colony location as fixed factor.

CHC difference between species occupying different elevational niches
A high proportion of CHC profile variation was explained by species (F = 56.90, df = 5, p < 0.001, R 2 = 0.83; Figure 3), resulting in speciesspecific CHC profiles (Figure 3; Supplementary material S8, PERMANOVA between all species pairs: p = 0.001). The species-specific compositions stayed consistent within species when individuals are depicted in a cluster dendrogram ( Figure 2B) with minor exceptions in B. lucorum and B. mucidus. The averaged CHC profile distances between species were not correlated with species phylogenetic distance (Mantel: R 2 = 0.08, p = 0.361). See Supplementary material S9 for mean proportion of components for each species.
Interspecific comparison of the mean proportions of saturated components in the CHC profiles revealed that B. pascuorum differed from the five other species, by having a lower mean proportion of saturated components (average: 45% ± 1.7%) than the other species (~56%; Figure 2C, glmm: χ 2 = 36.246, df = 5, p < 0.001). Moran's I correlation was significant at the phylogenetic distance of 0.02, representing the intraspecific level (indicated in red in Figure 2E). This suggest that proportions of saturated components are most similar among individuals within a species. The phylogenetic signal of this trait did not diverge from random expectations at the inter-species level at larger phylogenetic distances. This indicates that the proportion of saturated components show no pattern with phylogenetic distance of the considered species pool. Across all samples, neither MAT nor Tacclim affected the proportion of saturated components (MAT: χ 2 = 0.372, p = 0.541; Tacclim: χ 2 = 0.056, p = 0.815).

Intraspecific CHC variation along an elevational gradient in Bombus pascuorum
In B. pascuorum (n = 58) the CHC composition varied along the elevational gradient, but the degree of variance explained by elevation was very low (PERMANOVA: F = 3.293; p = 0.015, R 2 = 0.06; Figure 4A). Both climate-associated chemical traits, i.e., the proportion of saturated components ( Figure 4B) and the mean chain length ( Figure 4C), did not systematically change with elevation (glmm: χ 2 = 1.46; p = 0.226 and χ 2 = 0.39; p = 0.532). The coefficient of variation (CV) in the proportion of saturated components between individuals from the same site decreased with elevation (gam: ~elevation: edf = 1.94; DE = 85.6%; p < 0.001; Figure 4D)-a pattern driven by the highest sites above the tree line. The CV in mean chain length between individuals from the same site did not change with elevation (gam: ~ elevation: edf = 1; DE = 0.26%; p = 0.87; Figure 4E).

CHC profiles are species-specific
In this study, we explored the role of cuticular hydrocarbons in mediating bumble bees' capacity to cope with different climates in mountainous regions. Bumble bee species preferring different elevational niches revealed distinct and species-specific CHC profiles. Similarities between CHC profiles were not explained by similarities in species elevational niche preference, indicating that CHCs are at least not exclusively shaped by environmental factors correlating with elevation (e.g., mean annual temperature which correlates to 98% with elevation in this study). Rather, species specificity of the profiles suggests a strong genetic component in bumble bee CHCs. Genetic heritability of CHC patterns is known from Drosophila (Ferveur and Jallon, 1996;Holze et al., 2021;Ward and Moehring, 2021) and social insect species like ants and termites (Dronnet et al., 2006;Guillem et al., 2016). Notwithstanding the apparently strong genetic component in the CHCs, the similarity of CHC profiles did not reflect phylogenetic distances among species, suggesting that CHC profiles are under diverging selection within the pool of considered species.
Higher proportion of saturated components within CHCs have been shown to increase the desiccation resistance in some insects . In alpine bumble bees, five out of six investigated species did not differ in this trait, making it unlikely that the proportion of saturated components is of major importance in structuring the occurrence of bumble bee species along the elevational gradient. Only B. pascuorum, which tends to occur in the warmer foothills, expressed a lower proportion of saturated components than all other species. This contradicts the expectation that higher temperatures select for higher proportion of saturated components as improved desiccation barrier. Yet, it is partly in line with findings from other studies analyzing insect desiccation capacity along elevational gradients, which assume higher desiccation stress in higher elevations despite lower mean annual temperatures (Parkash et al., 2008;Mayr et al., 2021). Assuming that the proportion of saturated components really increases desiccation resistance in bumble bees, we suggest that low elevations may not select per se for a higher desiccation barrier, despite higher mean annual temperatures. Instead, less solar radiation, higher oxygen partial pressure, diverse microclimates and sufficient water supply (Körner, 2007) may release species like B. pascuorum from desiccation stress in the foothills. In particular, the canopy cover below the tree line can provide refugia (De Frenne et al., 2019), which allow bumble bees to actively seek places where desiccation stress is low.
Mean chain lengths differed between bumble bee species. Given that each additional carbon atom increases the melting temperature of CHC compounds by 1°C-3°C (Gibbs, 2002), differences in desiccation resistances between, e.g., B. soroeensis (on average 27.9 ± 0.13 carbons) and B. wurflenii (on average 25.5 ± 0.12 carbons) should be striking and require further examination-at best, under consideration of the Similarity of the CHC profiles of six bumble bee species and distribution of these species along the elevational gradient. (A) Similarity of CHC profiles of bumble bees displayed in a two-dimensional graph by non-metric multidimensional scaling (NMDS) based on Bray-Curtis distances. Dots represent CHC profiles of individuals, dot colors the respective species. The closer the dots, the more similar the CHC profile. (B) Bumble bee distribution along the elevational gradient recorded within an intensive monitoring conducted between June and September 2019 (see Supplementary material S1). The green line refers to the tree line in the study region. CHC profiles are clearly determined by species however, species profile similarity does not reflect similarity in the elevational distribution.
Frontiers in Ecology and Evolution 09 frontiersin.org absolute abundances of the different substances. Along the considered elevational gradient, however, differences in mean chain lengths did not predict the preferred elevational niche of bumble bee species. However, habitat preferences in a broader geographical context (Europe; Rasmont et al., 2015) which established 3.5 Ma ago (Hines, 2008) could explain differences in traits mediating desiccation resistance. Interestingly, Frontiers in Ecology and Evolution 10 frontiersin.org species differences in the mean chain lengths were more pronounced in closely related than in less related species ( Figure 2F). CHCs and mean chain lengths are known to be involved in the context of chemical communication in various insects (Nunes et al., 2009;Menzel and Schmitt, 2012;Leonhardt et al., 2016;Maihoff et al., 2022a). The negative correlation between trait similarity and phylogenetic distance suggests that mean chain length could play a role in the recognition of conspecifics, as a character displacement from closely related species might confer an evolutionary advantage in this context of communication. However, it should be noted, that the strength of a phylogenetic signal always depends on the taxonomic level studied and that the power of these analyses depends on the number of considered species (Blomberg et al., 2003;Menzel et al., 2017b), which restricts the interpretability of the phylogenetic signal in our study.

High intraspecific CHC variation and potential environmental filtering at high elevations in Bombus pascuorum
The CHC profile of B. pascuorum showed a high intraspecific variation along the 1,300 m elevational gradient. Profile variation between individuals collected at the same site was comparable to the variation between individuals collected at different elevations ( Figure 4A). We also found no evidence that either of the two climateassociated chemical traits was shaped by factors linked to elevation. We are aware that bumble bee workers of B. pascuorum can perform long vertical flights to exploit resources in different elevational belts (Knight et al., 2005). Such movements could impede the detection of intraspecific acclimatization in the CHCs of this species if the bumble bees are exposed to the changing environmental conditions for a shorter period of time than a change in CHCs can happen. Although the relevant time frame for changes in CHCs in bumble bees is still unknown, laboratory studies with flies indicate, that changes even within hours are conceivable (Stinziano et al., 2015). An experimental approach under controlled climatic conditions, testing short time acclimatization within hours, however, could shed light on this. Meanwhile, our findings suggest that the chemical composition of B. pascuorums' hydrocarbons does not seem to restrict individuals in the performance of such vertical flights in the considered elevational range.
Albeit CHC profiles and climate-associated chemical traits did not change systematically along the elevational gradient in B. pascuorum, we detected a decline in trait variance: Individuals collected in high elevations-especially above the tree line-differed less in the proportion of saturated components than individuals from low-or mid-elevations. A similar reduction of trait variation with elevation, and with the tree line as potential threshold, was detected in other functional traits, including bumble bee tongue length (Sponsler et al., 2022b) and wild bee body size (Classen et al., 2017). Pronounced environmental filtering at habitats with less variable microclimatic conditions (Hoiss et al., 2012) might cause such declines. Here, the phenotype that withstands the potential filtering process, shows a comparable low proportion of saturated components, i.e., on the intraspecific level the desiccation barrier seems rather reduced in individuals that occurred in high elevations. The environmental filtering suggests that the maintenance of phenotypes expressing fewer saturated components leads to an adaptive advantage, even though the mechanisms behind remain elusive. Genetic bottleneck effects at higher elevations (i.e., fewer colonies at high elevation equals less genetic variation) may further constrain trait variations in the highlands (Jackson et al., 2018).

CHC variation in translocated workers of Bombus lucorum
We found no effect of colony location on the CHC profile patterns of translocated B. lucorum workers. The translocation experiment showed that a 4-weeks exposure to different climatic regions, reflecting a 3-fold difference in annual precipitation (Supplementary material S2) and a 25% difference in mean temperature ( Figure 1B), did not trigger a systematic plasticity in the overall CHC composition as a response of warmer and dryer nest location. This is surprising, as precipitation/humidity and temperature differences can affect CHC composition in insects (e.g., ants) earlier-even within weeks-in both field and laboratory experiments Sprenger et al., 2018).
Colony location influenced chain length in translocated workers: workers translocated to warm and dry lowlands tended to have longer chains than those from the cool and wet mountainous foothills. This can be interpreted as a response to differences in regional climate, as longer chains are needed in warmer/drier conditions to prevent water loss through desiccation (Rajpurohit et al., 2017;Menzel et al., 2017a). Similarly, ants subjected to laboratory experiments increased their average chain length in warm conditions , and so did flies (Rajpurohit et al., 2017). However, in our study, mean chain lengths was not found to be shorter in even cooler habitats along the elevational gradient, suggesting again that other factors than temperature predict this trait. Further, chain length increased on average by only 0.3carbons, challenging the question about the physiological relevance of this effect. While chain length can also be influenced by age (Jackson and Bartelt, 1986;Nunes et al., 2009) we found this to be an unlikely effect here, as the colonies for the translocation experiment were reared at the same time, and thus workers were of comparable age. Nevertheless, including life history conditions and other thermal performance related functional traits, e.g., body size, would be informative in future studies as desiccation resistance is a complex trait, likely driven by multiple factors (Davis and Moyle, 2019). Other traits that may influence climate tolerance and hence the distribution pattern of bumble bee species along elevational gradients (Martinet et al., 2021) are genetic regulatory mechanisms like, e.g., the upregulation of heat shock proteins (Jackson et al., 2018;Pimsler et al., 2020) and pile color (Williams, 2007). Bumble bees-unlike ants-are characterized by exhaustive hairiness. The dense cover of hair might add insulation that reduces water loss in general (Kühsel et al., 2017). Thus, bumble bees might be less susceptible to desiccation stress (Parsons, 2019), which might reduce the pressure to alter CHC profile in response to desiccation.

Conclusion
In this study we highlight the high variability of cuticular chemical phenotypes in six species in the genus Bombus. Our results Frontiers in Ecology and Evolution 11 frontiersin.org suggest that the proportion of saturated components and mean chain length are not decisive traits explaining bumble bees' distribution along elevation in our study area. Profiles and chain lengths were rather species-specific, indicating that species recognition or long-established evolutionary divergence may limit trait variation. Yet, differences in profiles and climate-associated chemical traits may explain differences in species response to climate change whereas plastic adjustments within species seem to be unlikely. At present, species range shifts to higher elevations do not appear to be hindered by any particular CHC phenotypehowever, a link to species fitness is urgently needed to clarify this claim. With increasing temperatures or under extreme events, the differences in CHC composition of species might become relevant and determine species' resistance toward warming. Extending this study to more species, including those that already inhabit much warmer habitats, and under the consideration of multiple traits and related fitness consequences may shed light on bumble bee responses to climate change. Our study provides first insights that conserved species-specific signals are maintained under different environmental conditions. This is an important step toward deciphering how potentially conflicting functions in insect CHC profiles can be unified.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions
FM, AC, SN, and TS contributed to conception and design of the study. FM collected the specimens with help of KB, LB, SSc, and AC, performed the statistical analysis with input from SN and AC, and wrote the first draft of the manuscript. FM, SSa, SSc, KK, NS, and AC reared and maintained the bumble bee colonies. All authors contributed to the article and approved the submitted version.

Funding
The GC/MS was granted by the DFG to TS. This project was sponsored by the Bavarian State Ministry of Science and the Arts in the context of the Bavarian Climate Research Network (bayklif).