Introduction

Winter presents a significant challenge for honey bee colonies in temperate regions. Over the last two decades, beekeepers have reported high colony losses during winters. These losses have been associated with low temperature stress, geographic relocation due to climate change, genotype, diseases, pesticides, and poor nutrition1,2,3,4. In preparation for overwintering, a honey bee colony undergoes notable alterations in both the behavior and physiology of the worker bees. These changes include decreased individual activity, alterations in behavior regulated by endocrine systems, increased nutrient reserves, extended longevity, and thermoregulating clustering at the colony level1. Previous research has shown that honey bees live longer and grow larger when exposed to low temperatures5. However, other studies indicate that low temperatures negatively impact the honey bee immune system6and survival7through higher larval mortality and reduced life expectancy in adulthood8. To mitigate these deleterious effects, honey bees have evolved life history strategies like cessation of brood rearing during the winters, increased immunity, and increased antioxidant expression9. Additionally, honey bees form a dense winter cluster to withstand harsh winter conditions. The thermoregulation of this winter cluster relies heavily on the insulation provided by the tightly packed outer layer of bees, known as the mantle bees, which is crucial for survival at low temperatures10. Moreover, to reduce the cold temperature stress on the bee colonies, and to improve overwintering survival, beekeepers often store the hives indoors. Storing hives under controlled cold climate conditions of 5 °C—7 °C and 25% relative humidity has been shown to improve health and survival11,12,13.

The symbiotic gut microbiota of honey bees is vulnerable to temperature stress14. Winter losses of honey bee workers are primarily caused by harsh weather conditions, resulting in starvation and increased vulnerability to diseases15. Additionally, the honey bees experiencing cold stress have an increased risk of diseases and infections, which in turn raises the likelihood of colony losses. A recent transcriptomic study provides insight into the molecular mechanisms of cold resistance in Apis cerana. This study identified various temperature-sensitive proteins, such as heat shock proteins and zinc finger proteins, which were upregulated when bees were exposed to 0 °C and are considered candidate genes for honey bees’ tolerance to cold stress16. Apart from host physiology, the role played by the honey bee gut microbiome has also emerged as one of the most important physiological aspects to investigate in relation to honey bee cold tolerance. The gut bacterial species within honey bees undergo alterations due to various environmental and developmental stresses, with temperature being a crucial one17.

The honey bee gut microbiome contains five core host-specific bacterial species that are highly conserved and comprise 95% of the total microbes18,19. The hindgut of every adult worker contains Snodgrassella alvi, Gilliamella apicola, Lactobacillus Firm-5, Bombilactobacillus (formerly Lactobacillus Firm-4), and Bifidobacteriumspecies17. Along with these core bacteria, there are non-core species including Bartonella apis, Commensalibacter spp., and some other identified and several unidentified species that have been collected from the surrounding environment like hives and plants20,21. These bacterial species promote honey bee growth and physiology22, facilitate breakdown of toxic dietary compounds23, and modulate immune functions24. Dynamic changes to the microbiota composition have been observed throughout the seasons25,26,27. The gut microbiota has been reported to differ between winter and summer honey bees, with reduced α-diversity and higher levels of Bartonella and Commensalibacterduring winter28. Moreover, temperature during different months affects the composition of the gut microbiota by gradually changing the bacterial diversity between the seasons. In subtropical conditions, precipitation affects the composition of the honey bee gut microbiota27. Numerous studies have been conducted to understand the seasonal dynamics of honey bee gut microbiota. However, very few studies have explored the effects of prolonged exposure to harsh climatic conditions on honey bee gut microbiome. For instance, direct exposure to cold temperatures has been found to influence the gut microbiota in both honey bees and bumble bees29. Additionally, elevated temperature has shown to affect nectar microbes where the abundance of bacteria increased in the warmer temperature, influencing bumble bee forager preference30.

Beyond seasonality and temperature, vast differences in gut microbiota at the strain level has also been observed in two closely related honey bee species, Apis mellifera and Apis cerana31,32. Additionally, significant variation in both composition and function among diverse Asian honey bee populations has been reported, showing distinct patterns of their gut microbiota33. A previous study showed genetic divergence and functional convergence of gut bacteria in the eastern honey bee, Apis cerana and the western honey bee, Apis mellifera34 and between Apis cerana and Apis florea35. Differences in gut microbes between strains of Apis mellifera may influence the performance of those honey bee strains under different climate conditions.

Given that differences in gut microbes between honey bee strains may affect their performance under varying climate conditions, beekeepers will need to adopt different strategies to manage overwintering colonies effectively and prevent colony losses. Our study is important for beekeepers in two key ways. First, winter bees feed on pollen that has been stored in the hive for several weeks or months. Previous research has shown that consuming aged pollen affects the gut microbiota composition of nurse bees, leading to impaired development, increased mortality, and a higher risk of Nosema infection36. Second, the winter bees may retain feces in their gut for extended periods, which can likely affect nutrient availability and physicochemical conditions in the gut. Without defecation, bacteria may accumulate over time in the guts of winter bees, whereas more frequent defecation in nurses or foragers may result in faster microbiota turnover28. Hence, the goal of this study was to assess whether overwintering storage conditions disrupt the stability and diversity of the gut microbes in overwintering honey bees. To understand long term storage effects on gut microbiota, we compared the whole gut microbiota of two commercial strains of Apis mellifera i.e., claimed cold hardy bees Bolton Bees, MN (https://boltonbees.com/pages/mn-hardy-hives) that are bred and reared in Minnesota and Italian bees (Mann lake Bees, MN) of unknown rearing. Hives were stored either at a constant 6 °C indoors or kept outside in natural conditions during winter. We predicted that the hives stored outside in natural fluctuating temperatures will have different gut bacterial communities due to cold stress compared to hives stored at a constant temperatures. We also predicted that bacterial communities would change during the three-month sampling period of falling temperatures, viz. October, November, and December.

Material and methods

Bee samples

Honey bee hives were managed outside in the field (Fargo, North Dakota) from the first week of May to the last week of October and then 3 hives from each bee strain were shifted to refrigerated cargo containers for overwintering storage (6 °C) in November while 3 hives from each bee strain were left outdoors in the field. The first sampling was done on 10/17/2022 (base sampling of outdoor bees), followed by moving some hives inside storage on 10/21/2022. The second sampling was done on 11/15/2022 from hives both inside and outside. The third sampling was done 12/14/2022, again for both inside and outside hives (Fig. 1, Table 1). Hives were routinely inspected for signs of disease and mite infestations. While no diseases were detected, some hives with mite infestations were treated on 09/28/2022 using Mite Away Quickstrips (Apistan Anti-Varroa Mite Strips, Mann Lake Bee & Ag Supply). During the dearth period (when hives were moved to inside storage) in winter, the hives were supplemented with winter patties and Pro-Sweet liquid feed (Bee-Pro Patties, Mann Lake Bee & Ag Supply). Moreover, air temperatures were recordedfrom the North Dakota Agricultural Weather Network (NDAWN) weather station in north Fargo (latitude: 46.9003°, longitude: −96.8323°), near the filed site (Fig. 2). In addition, the temperature-managed cargo containers were frequently checked to ensure setpoint consistancy.

Fig. 1
figure 1

Illustration of the management of honey bee hives from October to December. In October, hives were managed outside in the field. In November, three hives from each bee strain were shifted to refrigerated cargo containers (6 °C) for overwintering storage, while another three hives from each strain remained outdoors in the field. The same management was repeated in December. Two genetic strains of bees were used in the experiment, Bolton Bees (B) and Mann Lake bees (M).

Full size image
Table 1 Representation of honey bee sample collected from various storage environments across different months.
Full size table
Fig. 2
figure 2

Hourly air temperature at the hive location.. All hives experienced outside air temperatures in October. Hives that remained outside and were not moved to indoor storage experienced these temperatures through the end of December. Dotted lines represent sampling dates, while the solid line marks the date on which half of the hives were moved indoors.Grey dots represents the average temperature for each day.

Full size image

DNA extraction

Bolton bees and Mann Lake bees were kept separate at two different locations in the field. The hives were rocked gently to prompt the bees to exit. The exiting bees were then captured at the entrance for sampling. Three bees from each treatment were used for whole gut dissection after surface sterilizing bees using 1% sodium hypochlorite followed by three washes using sterilized water in sterile conditions. During dissection, we ensured that the guts were intact and had not been perturbed or defecated prior to dissection. This served as a visual marker to maintain consistency in the sample collection process. Whole gut samples were prepared by bead-beating the samples on a Qiagen Tissuelyser for 6 min at 30 Hz to disrupt recalcitrant cells. Samples for cell lysis were prepared by combining two 3-mm chromium steel beads and approximately 50 μl of 0.1 & 0.5 mm ZR BashingBeads™ inside lysis tube (Biospec, Bartlesville, OK) with 750 μl of ZymoBIOMICS™ lysis solution and 20 μl of proteinase K. A subsequent round of bead beating involved rotating the samples for 6 min at 30 Hz, followed by an incubation period at 56 °C for one hour. DNA was extracted using ZymoBIOMICS™ DNA microprep kit collection (Zymo research, Irvine, California) including 2 blank extractions as a no template control for further downstream analysis.

PCR amplification and illumina sequencing

16S rRNA gene libraries for paired-end reads were generated following a previously published protocol37,38. The V5-V6 region of the 16S rRNA gene was amplified using 16S rRNA gene primers (799F mod3: CMGGATTAGATACCCKGG and 1115R: AGGGTTGCGCTCGTTG), each incorporating a unique barcode sequence. PCR1 reactions were conducted with 2 μl of DNA, 10 μl of 2 × Pfusion High-Fidelity DNA polymerase (New England Biolabs, Ipswich, MA), 10 μl of ultrapure water, and 0.5 μl of 10 μM 799F-mod3 and 0.5 μl of 10 μM 1115R primers. The reaction condition for PCR1 was 94 °C for 3 min, 94 °C for 45 s, 52 °C for 1 min, 70 °C for 1:30 min, repeated step2 29X, and 72 °C for 10 min to amplify this region. We ran gel electrophoresis to confirmed if 16S rRNA primers were successfully attached to our samples. To complete the Illumina adapter sequence, we initiated the process by cleaning the PCR product with exonuclease and shrimp alkaline phosphatase. Exonuclease was applied to remove excess primers, while shrimp alkaline phosphatase was used to eliminate residual deoxynucleoside triphosphates (dNTPs). Four µL of the 1 × ExoSAP was mixed with each 7µL PCR1 sample. This reaction was incubated at 37˚C for 30 min and then 95˚C for 5 min. Following this initial cleanup, the purified PCR products were employed as templates for a second PCR. This subsequent PCR (PCR2) utilized 1 μl of the cleaned PCR product as a template, using the same primers (PCR2F and PCR2R), and was conducted under conditions identical to the initial PCR. For PCR2, a 1 μl aliquot of Exo-Sap (Thermo Fisher Scientific) PCR1 product was utilized to conduct the second step of Illumina library preparation. In PCR2, the Exo-Sap PCR1 products were further amplified, incorporating linker poly-A primers for recognition on the sequencing platform. The PCR2 reaction conditions included an initial denaturation at 94 °C for 3 min, followed by cycles of 94 °C for 45 s, 58 °C for 1 min, 72 °C for 1 min and 30 s, with step 2 repeated 14 times, and a final extension at 72 °C for 10 min, with final volume of 25 μl. To ensure uniform DNA input quantities, Invitrogen DynaMag TM-96 Side Skirted kit was employed for the normalization step across all samples prior to submitting a final volume of 10 µl for sequencing. The 16S rRNA gene sequencing was performed on NextSeq2000 P2 600 cycle kit (2 × 300xN/A) that produced total of 382.53 M reads at Q30 of 88.96%.

Bioinformatics

QIIME 2–2019 was used for the visualization and trimming of low-quality ends in reads from raw 16S rRNA sequence libraries. Subsequently, DADA2 was utilized for the assignment of sequences to amplicon sequence variants (ASVs), which represent 16S rRNA gene sequences with 100% matches. Taxonomy was assigned to the ASVs through the sklearn classifier, trained specifically for the 799–1,115 region of the 16S rRNA gene, utilizing the SILVA 138.2 database (accessed November 2023). Additionally, local BLASTn searches against the NCBI 16S microbial database (accessed October 2023) were conducted. Features were filtered from the resultant ASV table, corresponding to contaminants identified in the blanks using the R package ‘decontam’ (version 1.10.0) at a conservative threshold of 0.5. This process aimed to identify contaminants while also eliminating ASVs identified as chloroplast and mitochondria. To standardize the number of sequences per library, alpha rarefaction in QIIME2 was employed, and 4000 reads per sample were selected. This approach ensured the retention of 107 out of 108 samples while capturing most of the diversity.

Generalized Linear Mixed-effects Models (GLMMs) were used to model the Shannon diversity index. All GLMMs were fitted using the lme4 package39and lmer Test package40in R41, using R studio version 4.2.2. Bee strain, storage treatment (exposure to natural vs. storage at 6 °C temperature), and sampling month were incorporated as fixed effects. Hive ID was incorporated as a random effect to account for the repeated measuring of each hive42. Hive ID corresponds to the hives, while the sample ID refers to the samples taken from those hives. A stepwise backward selection process was employed for model selection. Initially, a full model was run with all fixed effects and random effects. Then all possible interactions among fixed effects were examined, followed by the systematic removal of one fixed effect at a time. The process continued until no further model improvement (e.g., lower AIC) could be observed. A difference of ≤ 2 in the AIC indicated that models were similar, in which case the simpler model (i.e., with fewer parameters) was selected as the best-fit model43. The best fit model was diagnosed for model assumptions through residual diagnostic plots using the ‘DHARMa’ package44. Finally, package ‘ggeffects’ version 0.16.045 was used to extract means, standard deviations, and 95% confidence intervals for the predicted values of Shannon diversity index for different sampling times. Tukey HSD tests for multiple comparisons between treatments was employed with estimated marginal means comparisons (EMMs) using the ‘emmeans’ function46. Other matrices like observed species richness, shannon diversity, and simpson diversity are provided in supplementary data (Data available; alpha amp.csv, Supplementary Fig. 1).

For beta diversity analysis, Adonis Bray–Curtis distance dissimilarities and nonmetric multidimensional scaling (NMDS) ordination were conducted in R v4.3.0 with ‘vegan’ package. The ‘phyloseq’ (ver 1.34.0) R package was used for beta diversity analysis. Permutation analysis of variance (PERMANOVA) was applied using the ‘adonis’ function on distance matrices with 999 permutations. The differences between bee gut microbiomes including distance as a dependent variable and months (October, November, and December), bee strain (Bolton bees and Mann lake bees) and Storage Status (Inside and Outside) as independent variable were assessed using separate Adonis models. In all models, bee hives were accounted as a block- strata, where permutations were allowed among levels of all the variable, but within each level of bee hives (random variable), no permutations occurred.

To understand the effect of different variables on the abundance of candidate bacterial species, Generalized Linear Mixed-effects Models (GLMMs) was used to model the abundance of particular bacterial species. Bee strain, storage treatment (exposure to natural vs. storage at 6 °C), and sampling month were incorporated as fixed effects. Hive ID was incorporated as a random effect to account for the repeated measures of each hive42. The fitted GLMM and LM was diagnosed for issues with heterogeneity, overdispersion or missing co-variates using the residual plots. To evaluate the significance of fixed effect, analysis of variance tables for the GLMM was calculated using Wald chi-square tests. Tukey HSD tests was employed for multiple comparisons between treatments with estimated marginal means comparisons (EMMs) using the ‘emmeans’ function46.

Results

Alpha diversity

To investigate alpha diversity patterns, generalized linear mixed models (LMM) that included Shannon diversity as a response variable was applied (Table 2). The most parsimonious, best-fit model included month as a fixed effect, and random effects for hive ID (different hives in each treatment) (Table 3). In contrast to the predictions, alternative models incorporating bee strain and storage treatment did not significantly improve the best-fit model. The random effects structure in the best fit model showed relatively little variation among hives (SD = 0.1604), when compared to the residual error of the model (SD = 0.3560). This indicates variation between hives and the replications do not account for the majority of the variation in Shannon diversity. Furthermore, the Shannon diversity index in October differed significantly from November (p value = 0.0482), while December showed no significant difference (p = 0.1891) (Data available; alpha diversty R file). Month explained only 8% of the variation in Shannon diversity, whereas the full model explained 23% of the variation in Shannon diversity. Overall, results indicate a temporal variation in Shannon diversity, but the strain of bee and whether the hives were stored inside or outside had no effect (Fig. 3). Non modelled matrices like observed species richness, Shannon diversity, and Simpson diversity were not effected by any of the predictor variable (p > 0.05).

Table 2 Comparisons of Akaike Information Criterion (AIC) value for different Generalized Linear Mixed Models (GLMMs) with relevant predictor and random variables.
Full size table
Table 3 Results of best fit Generalized Linear Mixed Model from Shannon diversity index.
Full size table
Fig. 3
figure 3

Shannon diversity index by month, bee strain and storage status. Shannon diversity was significantly different in November (p = 0.0482) as compared to October. October and December were not significantly different (p = 0.1891). Neither bee strain nor storage condition explained Shannon Diversity (see Table 2). Model predictions for effect of sampling time on Shannon diversity index (see Table 3). In A and B, the line inside the box represents the median, while the whiskers represent the lowest and highest values within the 1.5 interquartile range. Individual sample values are shown as points. In C, grey shading represents the 95% confidence interval.

Full size image

Beta diversity

To investigate beta diversity patterns, we conducted PERMANOVA using the adonis2 package in R, testing various models with Bray–Curtis distance matrices (phyloseq) as a response variable. The simplest and best-fitting model included storage status and month in interaction with bee strain as a fixed effect, and block (strata) as bee hive ID; Distance ~ Storage status + Month * Bee strain, and strata = hive ID (Table 4). The overall composition of bee gut microbiomes differed by bee type (adonis F = 3.19, df = 1, R2 = 0.0272, p value = 0.001), month (adonis F = 5.09, df = 2, R2 = 0.0872, p value = 0.001), but not storage status (adonis F = 0.57, df = 1, R2 = 0.0048, p value = 0.887). The interaction between month and bee strain was significant (adonis F = 1.57, df = 2, R2 = 0.0269, p value = 0.021). Month explained only 8.7% of the variation, whereas the bee type, and interaction between month and bee strain explained 2.7%, and 2.6% of the variation, respectively. Overall, these findings suggest that the microbial dynamics in overwintering honey bee workers are influenced by bee strain over the months, with stability observed in storage conditions (Fig. 4). Pairwise comparison showed significant differences in gut microbiota when compared October – November (p value = 0.002), November – December (p value = 0.001), and October—December (p value = 0.001) where a significant decline in beta diversity in November was observed.

Table 4 Best fit permutation analysis of variance (PERMANOVA) model applied using the adonis function with distance matrices with 999 permutations and bee hives accounted as a blocks – strata.
Full size table
Fig. 4
figure 4

Non-metric multidimensional scaling (NMDS) performed with Bray–Curtis dissimilarity. Month (p = 0.001) and bee strain (p = 0.001) significantly affected beta diversity, and their interaction was also significant (p = 0.021). Each point on the plot represents a honey bee gut sample and samples that share greater similarity are ordinated closer together. Ellipses on the plot indicate 95% confidence intervals.

Full size image

Taxonomical abundance

Following the removal of barcodes and primer sequences through QIIME2 (Demultiplexing), we obtained a collective sum of 34,879,494 paired-end 16S rRNA gene amplicon reads. The average read count per sample was 317,086, with an overall average quality score of 34. Across the 108 samples analyzed, the range of reads per sample varied from a minimum of 56,533 to a maximum of 1,314,819. Analyzing the taxonomic profile revealed a predominance of Firmicutes (70.96%), Proteobacteria (20.18%), and Actinobacteriota (5.83%) at the phylum level. We investigated the presence of the core gut microbiota and identified only a limited number of reads for Snodgrassella.

In contrast to Snodgrassella, we identified various species within the genus Lactobacillus, including unclassified Lactobacillus, Bombilactobacillus mellis (formerly Firm-4 group), Lactobacillus apis (Firm-5 group), Apilactobacillus micheneri, and Lactobacillus uncultured Firmicutes, which dominated across all bee samples. In this study we focused on examining the relative abundance of all the species under Lactobacillus genus and two non-core (Bartonella and Commensalibacter) and two core (Gilliamella and Bifidobacterium) microbes, along with the endosymbiont Wolbachia. Moreover, considering the importance of genus Bartonella and Commensalibacterin a previous study28, we analyzed the abundance of two Bartonella species: Bartonella apis and another denoted as uncultured Bartonella and two Commensalibacter species: Acetobacteraceae_bacterium and another denoted as Unclassified Commensalibacter. The bacterial species-specific analysis for relative abundance was performed using mixed model with formula: Abundance ~ Storage status + Month + Bee strain + (1 | hive ID).

Our analysis revealed that the month significantly influenced the abundance of unclassified Lactobacillus (df = 99.29, F-value = 31.8402, p-value < 0.005). Specifically, October showed a significantly higher abundance compared to December (SE = 3.56, df = 112.7, p value = 0.0261), while October had a significantly lower abundance compared to November (SE = 3.56, df = 112.7, p value = 0.0002). Additionally, December had a significantly lower abundance compared to November (SE = 3.10, df = 97.5, p value < 0.0001) (Fig. 5). Neither storage status (df = 28.17, F-value = 0.1432, p value = 0.7079) nor bee strain (df = 12.15, F-value = 1.0835, p value = 0.3182) significantly affected the abundance of unclassified Lactobacillus. The abundance of Lactobacillus apis did not vary by storage status (df = 107, F-value = 0.885, p value = 0.3489) or bee strain (df = 107, F-value = 1.4531, p value = 0.2306), but month had a significant effect (df = 107, F-value = 4.7334, p value = 0.0107). We found Lactobacillus apis had significantly higher abundance in October compared to December (SE = 1.50, df = 110, p value = 0.0125), whereas no significant difference was observed between October and November (SE = 1.50, df = 110, p value = 0.0673), nor between November and December (SE = 1.34, df = 102, p value = 0.7473). None of the predictors (storage status, month, or bee strain) significantly affected the abundance of Bombilactobacillus mellis (storage status: df = 107, F-value = 0.2015, p = 0.654; month: df = 107, F-value = 2.530, p value = 0.0844; bee strain: df = 107, F-value = 0.6398, p value = 0.4256), Apilactobacillus micheneri (storage status: df = 107, F-value = 2.0463, p value = 0.1555; month: df = 107, F-value = 0.5612, p value = 0.5722; bee strain: df = 107, F-value = 2.1700, p value = 0.1437), and Lactobacillus uncultured Firmicutes (storage status: df = 107, F-value = 1.2865, p value = 0.2592; month: df = 107, F-value = 1.680, p = 0.1896; bee strain: df = 107, F-value = 3.987, p value = 0.05).

Fig. 5
figure 5

Relative abundances of Unclassified Lactobacillus present in honey bee gut across different months. Each column represents an individual bee. The relative abundance, represented in percentages, is shown on the y-axis.

Full size image

Bartonella apis abundance did not differ by month (df = 99.11, F value = 2.27, p value = 0.1085), bee strain (df = 38.26, F value = 0.1824, p value = 0.6717), or storage status (df = 38.26, F value = 0.1824, p value = 0.6717). However, the Bartonella uncultured ASV significantly differed by month (df = 98.72, F value = 4.74, p value = 0.0107) while storage status (df = 34.58, F value = 0.0528, p value = 0.81954) and bee strain (df = 11.559, F value = 0.1322, p value = 0.72274) were not significant. We found a significant difference in abundance between December and November (SE = 1.03, df = 97, p value = 0.0085), but not between October and December (SE = 1.21, df = 112, p value = 0.358) or October and November (SE = 1.21, df = 112, p value = 0.456). For Commensalibacter (Commensalibacter Acetobacteraceae_bacterium), bee strain (SE = 1.48, df = 102, p value < 0.0001) was the only factor significantly influencing abundance with Mann Lake bees having a higher abundance than Bolton bees. The effect of storage status (df = 107, F value = 0.261, p value = 0.6097) and month (df = 107, F value = 1.6663, p value = 0.1938) was not statistically significant. Moreover, none of the predictors significantly influenced the abundance of unclassified Commensalibacter (storage status: df = 107, F-value = 0.4946, p value = 0.4834; month: df = 107, F-value = 0.1538, p value = 0.8577; bee strain: df = 107, F-value = 1.0855, p value = 0.2998).

Bee strain was found to be a significant predictor of Bifidobacterium (Bifidobacterium uncultured_Bifidobacterium) abundance where Bifidobacterium showed higher abundance in Bolton bees compared to Mann Lake bees (SE = 0.077, df = 102, p value = 0.0014). The storage status (df = 107, F value = 0.0682, p value = 0.7944) and month (df = 107, F value = 0.8069, p value = 0.4489) variables did not contribute significantly to the variation in this model. The storage status (df = 37.34, F value = 6.67, p value = 0.013) and month (df = 99.15, F value = 14.18, p value = 3.808 × 10−6) significantly affected the abundance of Gilliamella (Unclassified Gilliamella), while Bee strain did not show a significant effect (df = 12.10, F value = 0.675, p value = 0.426). Pairwise comparison showed significant differences in abundance during October and December (SE = 0.197, df = 112, p value < 0.0001), November and December (SE = 0.166, df = 97.5, p value = 0.002), but no significant difference during November and October (SE = 0.197, df = 112, p value = 0.5361). A significant effect of the storage status on the abundance of Gilliamella was seen between the Inside and Outside bees (SE = 0.215, df = 45, p value = 0.0272). storage status (df = 48.66, F value = 0.0377, p value = 0.846) and Month (df = 99.42, F value = 0.97, p value = 0.38) is not a significant predictor of abundance in Wolbachia (Unclassified Wolbachia). Wolbachia abundance was significantly different between the bee strains (df = 11.78, F value = 5.352 p value = 0.03959). This results suggest a potential difference in gut microbial abundance between Bolton bees and Mann lake bees, but the significance was marginal (SE = 0.243, df = 14, p value = 0.0570). Overall, this taxonomic data emphasizes the marginal differences in gut microbiome composition between bee types, revealing a temporal changes in microbial composition that are not influenced by storage conditions in overwintering worker honey bee (Fig. 6A, 6B, 6C).

Fig. 6
figure 6figure 6figure 6

Relative abundances of the top fifteen bacterial species present in honey bee gut in different storage conditions in October (A), November (B), and December (C). Less than 10% is a category of low abundance species that made up less than 10% of the median number of reads in each gut samples. Each column represents an individual bee. The relative abundance, represented in percentages, is shown on the y-axis.

Full size image

Discussion

We used a 16S rRNA gene amplicon sequencing approach to study/characterize the gut microbiota of overwintering adult worker bees between (a) two commercial bees (Bolton Bees and Mann Lake bees), (b) two storage conditions and, (c) three months during a season. In this study, we hypothesized that temperature would be one of the most important environmental factors affecting the gut microbiota in overwintering honey bees. Surprisingly, we found no changes in the diversity, composition, or abundance of the microbiota when the bees were stored at constant 6 °C or outside during the winter. Interestingly, we found that the overwintering bees’ gut microbiome was dominated by the genus Lactobacillus, contributing to an increased relative bacterial abundance across all samples, regardless of storage conditions. An increase in Lactobacillushad been previously detected at relatively lower abundances in earlier studies on winter bees47,48.

The dominant Lactobacilli in our study included unclassified Lactobacillus, Bombilactobacillus mellis, Lactobacillus apis, Lactobacillus uncultured Firmicutes, and Apilactobacillus micheneri. Recently, Zheng et al., 2020 reclassified the genus Lactobacillusinto 25 distinct genera, 23 of which are newly established based on phenotypic and genotypic characteristics49. Three previously described species, along with Bombilactobacillus mellis (formerly Lactobacillus mellis50), have been assigned to the newly defined genus Bombilactobacillus, which belongs to the Firm-4 group51,52,53,54. Another major group within the reclassified bacterial generus is the former Lactobacillus Firm-5, which includes Lactobacillus apisas a member55. Additionally, Lactobacillus kosoi has been recently reclassified as a heterotypic synonym for Apilactobacillus micheneri56, a species found abundantly in solitary bees. Each of these species has been studied for their distinct functions. For example, Lactobacillus apis(Firm-5) has been shown to convert tryptophan to an indolic aryl hydrocarbon receptor agonist, improving learning and memory behaviors in honey bees55. Honey bees rely on the bacterial metabolite succinate, primarily derived from Lactobacillus apis, to maintain energy metabolic homeostasis57. Both Firm-4 and Firm-5 are closely related Firmicutes58, with Lactobacillus apis having well-documented distinct functions, while Bombilactobacillus mellis remains relatively understudied. Apilactobacillus, the most abundant microbe in the solitary bee microbiome, is environmentally acquired59and exhibits antibiotic resistance60. An unbalanced microbial community can regulate the abundance of Apilactobacillus micheneri, potentially influencing solitary bee biology and compromising their growth and survival59. Moreover, many studies have found genus Lactobacillusserves as the most important genus in the guts of honey bees61,62that plays an important role in the health of bees by showing probiotic characteristics63,64,65. Lactobacillushas shown numerous function in regulating honey bee health like stimulating immune system66, nutrient metabolism64, ability to digest flavonoids and other compounds present in the pollen wall67, and to inhibit pathogens62,68,69. Moreover Lactobacillihave been shown to release short chain fatty acids and vitamins which are utilized by midgut cells as an energy source70,71and play a role in regulating metabolism67. This suggests that the higher abundance of Lactobacillus during the overwintering phase may regulate host energy metabolism and support heat production by worker bees, thereby stabilizing colony temperature. Despite the physiological stress induced by cold temperatures, this study demonstrates that honey bees have the ability to acclimatize their gut microbiota in cold conditions. One possible explanation for these findings is that honey bees rely on specific bacterial species to maintain their health, utilizing the functions performed by the gut microbiota. Thus, bees may retain these bacterial species, particularly those belonging to Lactobacillus, to extract sufficient resources at low temperatures, regardless of bee strain or storage conditions.

We observed differences in the abundance of several bacterial species between the two bee types. Bolton bees exhibited significantly higher abundance of Bartonella (in some instances), Bifidobacterium, and Wolbachia, whereas Mann lake bees showed significantly higher abundance of Commensalibacter Acetobacteraceae_bacterium when compared to Bolton bees. In particular, the bees obtained from Bolton Bees represent specific genetic lines known for their resilience in Minnesota’s harsh climate. This genetic line, named after the location, where the parent Queen overwintered (Bolton, Minnesota), is derived from the base stock called MN-Hardy. The queens from the MN-Hardy line, used as grafts for queen production, have successfully endured the challenges of a prolonged and frigid Minnesota winter. Typically, such severe winters could affect honey stores, cleansing flights, and brood-laying, but these MN-Hardy Queens are claimed to have demonstrated ability to thrive under these cold winter conditions (https://boltonbees.com/pages/mn-hardy-hives). On the other hand, “OHB Italian” Mann Lake bees, sourced from Olivarez Honey Bees (OHB) in Northern California, are accompanied by claims for their hive performance, disease resistance, and overall robust health. The Italian Queens from OHB have reportedly undergone extensive breeding and careful genetic selection. The incorporation of traits such as “Minnesota Hygienic” and “VSH” (Varroa Sensitive Hygiene) further strengthens their innate ability to resist diseases and combat mite infestations. These two distinct types of bees have been specifically employed for their resilience to withstand the challenging winters typical of the Midwest, making them interesting subjects to study the gut microbiome during the overwintering period. In a previous study on honey bees, the colonization of specific sets of bacteria, such as Lactobacillus Firm5 and BifidobacteriumBifido-1.2, in genetically varied hosts strongly suggests that the genotype significantly influences the microbiota structure32. Another study on honey bees has identified marked differences in the core gut microbial community when comparing different lineages, which include Maltese honey bees (lineage A) to the Italian honey bees (lineage C). Notably, Maltese honey bees exhibited an inverse proportion of Lactobacillaceae and Bartonellaceaewhen compared to Italian honey bees72. The functionality of these differences between bee strains have not been studied. This finding underscores the extent of strain-level diversity within the bacterial communities. However, our study showed extensive overlap of the gut microbial diversity between bee strains and with significant differences in only a few bacterial species specifically, Bartonella exhibited higher abundance levels in Bolton bees, while Mann Lake bees showed increased abundance of Commensalibacter.

Several studies have demonstrated differences in gut bacterial communities between bee adults based on geographical locations at the same point in time33,73,74. However, very few studies have explored differences between bee genotypes based on their geographical origins. Given that most colony loss occurs during the winter months and that honey bee gut microbial communities fluctuate seasonally17, significant differences exist between the gut microbiota of summer and winter bees. Winter bees in previous studies have shown a dominance of Bartonella and Commensalibacter28. Additonally, the variation in Lactobacillusspecies may be influenced by winter conditions in different geographical regions. For instance, the study conducted at the University of Lausanne, Switzerland28 reported less extreme cold temperatures compared to Fargo, ND, where our study was conducted (average temperature in Oct- 13.33 °C, Nov- 8.8 °C and Dec- 3.3 °C, data collected from weatherspark.com). The geographical variation could account for the higher abundance of Lactobacillusobserved in the Fargo, ND study compared to the Lausanne, Switzerland study. Lausanne (46.5197° N, 6.6323° E, elevation 1,624 ft), located on the northern shore of Lake Geneva, features a hilly landscape, whereas Fargo (46.8772° N, 96.7898° W, elevation 906 ft) is situated in a flat region predominantly covered by agricultural land. These differences in topography and land use may influence the microbial communities. Therefore, research on the gut bacteria of winter bees is still in its early stages, with only a few studies available47,75,76. Additionally, a recent study found that Lactobacilluswas one of the most abundant bacterial genera in hives that either survived or failed to survive the winter of 2022 in northern Virginia, USA, hence supporting its role in winter survival. Previous studies suggested that the broad environment77, hive location78, artificial pollen supplements79have subtle influence on the relative abundance of the honey bee microbiota. These findings suggest that the higher abundance of Lactobacilli may be attributed to the significant impact of environmental acquisition routes like solitary bees, just like in solitary bees59. Moreover, in-depth analysis of this study revealed that Commensalibacter and Snodgrassellawere positively associated with winter hive survival, underscoring their importance for hive resilience76 Similarly, in our study, we observed that Commensalibacter exhibited differential relative abundance between bee strains, with Lactobacillus dominating in all winter bees. However, we detected little to no Snodgrassella in our samples. Therefore, this study observed that the winter honey bee gut microbiota differs from that in other studies conducted across different months, particularly under natural winter conditions in the midwestern United States.

Overall, the findings in this study reveal varying degrees of presence of bacterial species in overwintering honey bees across different months. There were significant differences in the richness and evenness (Shannon Diversity Index) between October and November, but no differences were noted in the bees sampled in November when compared to December samples. Additionally, beta diversity shows difference in species composition and relative abundance in different months. Workers showed significant differences in species composition across months, with a significant variation in beta diversity, consistent with previous studies47,80. Several studies have examined the gut microbiota composition in various types of worker bees across different seasons, consistently demonstrating that the microbial community structure remains relatively stable over time79,81,82. However, recent study have noted that the gut microbiota differs between winter and summer honey bees, with the long-lived winter bees exhibiting a stable microbiota with reduced α-diversity and higher levels of Bartonella and Commensalibacter28. This lower community alpha diversity with Bartonella and Commensalibacter as dominant bacterial species may confer certain physiological benefits. Hence, these studies also collectively reveal minor variations in gut microbial communities in temperate honey bee colonies during winter, indicating a shift in dominance, with the non-core bacterium Bartonellasurpassing the core bacterial species28,80,83. Another recent study also highlighted the significance of Bartonella, highlighting its expanded capability to convert metabolic wastes such as lactate and ethanol into pyruvate, which potentially provides energy for the host as well as other symbionts75. Considering the importance of Bartonella in previous studies, we analyzed two different species of Bartonella for their taxonomic abundance. However, we found the Lactobacillus genus as dominant, present in higher abundance, contrasting with the findings of a previous study where the Bartonella genus showed dominance in winter bees.

The role of Lactobacillus in overwintering honey bee workers in this study is limited by the nature of 16S rRNA sequencing. Additionally, 16S rRNA sequencing leaves a gap in understanding the absolute abundances of bacterial species as it only reliably provides relative abundances. To improve bacterial species detection resolution and functional profiling at the gene level, future studies could utilize qPCR or shotgun sequencing instead of 16S amplicon sequencing. Additionally, our study was limited by its timeframe, as sampling was conducted only in the month of October, November, and December. We observed a significant decline in beta diversity in November, which may vary if bees are sampled in later months (e.g., January, February, and March) as temperatures continue to drop during these colder months.

Conclusion

In conclusion, our study shows that the gut microbiota of overwintering honey bees is dominated by one specific bacterial genus, Lactobacillus, which likely contribute to the bees’ ability to withstand environmental stressors61,62,67,68,69,70. The results of this study raise two key questions. (1) Why was the Lactobacillus genus predominant in winter bees and why was this pattern not observed in previous studies? Could this genus contribute to honey bee survival under extreme cold conditions, such as those found in North Dakota? (2) While the overall gut microbiota remained relatively stable, there were slight shifts in the composition of specific host-associated bacterial communities across genetically distinct honey bee colonies. Do these subtle compositional differences have functional implications that influence colony health or survival? Overall, our study adds to the growing body of research on honey bee gut microbiota and its role in bee health during winter. Future research should focus on specific Lactobacillus species to better understand the microbial dynamics that support their functional roles, and to explore the mechanisms underlying the resilience of honey bee gut bacteria and its implications for overall bee health. Lastly, but importantly, only a few studies have reported the presence of Wolbachia in honey bees, including findings detailed in84,85. Interestingly, Wolbachia, which is exclusively maternally transmitted and known to alter host reproduction through various mechanisms, has been detected in the gut of these winter bees. Although its abundance is notably low, this observation is surprising and intriguing. This finding opens new avenues for research on the prevalence of Wolbachia in the honey bee gut and highlights the need for future studies to explore the potential relationship between Wolbachia and honey bee workers.