Symbiosis with microorganisms has been a crucial aspect of evolution, particularly in the divergence of insects, allowing them to colonize all land habitats and occupy all the trophic niches (Cornwallis et al., 2023; Lemoine et al., 2020). Symbiotic relationships between bacteria and their hosts can be classified on various levels, depending either on the stability of host-symbiont association (Perreau and Moran, 2022) or the nature of the relationship (beneficial, neutral, or harmful) (Hammer et al., 2019). Many insect symbioses are long-established, dating back millions of years, and played a key role in insect evolution. On the other hand, we have many examples of relatively new symbiont acquisitions. However, there is no doubt that both ancestral and newly acquired symbionts are still important players in shaping insect biology, as well as their adaptation to changing environmental pressures and that this is an ongoing process (Łukasik and Kolasa, 2024). The implementation of sequencing techniques in symbioses studies enabled a better understanding of many aspects of their functioning.
With the growing availability of advanced molecular techniques the scientific community has begun to ask the question of how microbial dynamics might influence the behavioral variation in social insects (Jeanson and Weidenmüller, 2014). Since then, many eusocial insects' microbiomes have been described in the context of their structure and plasticity (Suenami et al., 2023). Nonetheless, it seems that the best model species for studying the gut-brain axis are ants (Kay et al., 2023), bumblebees (Bombus) and western honeybees (Apis mellifera) due to the relative complexity of their core microbiomes (Liberti et al., 2024) and potential horizontal transfer of symbionts (Kwong et al., 2017). Honeybees have become a preferred model species due to their economic importance, well-known biology and relatively easy breeding. Their gut microbiota has been described as highly stable, consisting of only nine taxa (Engel and Moran, 2013), with most of the composition patterns and functions being performed by only five core taxa: Snodgrassella alvi, Gilliamella apicola, Bombilactobacillus (previously called Lactobacillus Firm-4), Lactobacillus (previously called Lactobacillus Firm-5), and Bifidobacterium asteroides (Motta and Moran, 2024; Zheng et al., 2020). However, variations in workers' microbiomes have been reported depending on the season (Kešnerová et al., 2020; Ludvigsen et al., 2015), diet (Julia C Jones et al., 2018), or geography (Ludvigsen et al., 2017). Recent studies also highlight the role of specific symbiont strains in nutrition provisioning and pathogen resistance (Miller et al., 2021; Parish et al., 2022). However, the honey bee symbionts' role is not limited only to triggering immune responses or coping with environmental stressors.
Recent studies suggest that gut microbiota can influence various aspects of honey bees' behavior. Zhang et al. (2022a) experimentally showed how Lactobacillus modulates learning and memory behavior by regulating tryptophan metabolism. Another study (Zhang et al., 2022b) showed how different bacterial taxa regulate specific modules of metabolites in workers' hemolymph. For example, Bombilactobacillus and Lactobacillus altered amino acid metabolism pathways, leading to upregulations of genes responsible for olfactory functions and division of labour (Zhang et al., 2022a). At the same time, Gilliamella was mainly responsible for circulating metabolites involved in carbohydrate and glycerophospholipid metabolism pathways (Zhang et al., 2022b). Another example is an experiment consisting of monocolonization with Bifidobacterium asteroides, which showed that worker bees colonized with this bacterium had elevated gut concentrations of juvenile hormone III derivatives (Kešnerová et al., 2017). This hormone is responsible for insect growth, development, and reproduction, and in honey bees, is a key factor for transitioning from nurse bees to foraging bees. More recent studies (e.g., by Vernier et al., 2024) showed that a single bacterial taxon can have a significant effect on bee behavior. Namely, they showed that inoculation with Bifidobacterium asteroides significantly increased the frequency of bees' foraging compared with microbiota-depleted controls. Contrary, inoculation with Bombilactobacillus mellis, and Lactobacillus melliventris lowered such behavior. However, the behavioral effects of symbionts are not limited to single individuals. Liberti et al. (2022) compared the social behavior of microbiota-colonized and microbiota-deprived bees. They found that colonized bees had increased levels of brain metabolites (such as ornithine and serine) and increased specialized head-to-head interactions between nestmates. Their complex behavior system and relatively stable microbiota make honey bees an excellent and promising system for studying the gut-brain axis (Liberti and Engel, 2020).
Due to their social nature, honey bees have become a model species in behavioral studies. The most widely known phenomenon is their ability to communicate the location of food sources by dancing (Von Frisch, 1967). However, one phenomenon remains relatively poorly understood - reproductive swarming (hereafter swarming). This complex and multi-step process leads to the division of the colony, with an old queen leaving the hive with approximately 75% of bees (Fefferman and Starks, 2006) and a new queen emerging. Typically, it occurs in mid-spring when the food sources are most abundant. Although scientific literature points out mass queen rearing as the first visible sign of swarming preparation (Grozinger et al., 2014), beekeepers mention other signs that inform them about upcoming swarming, such as the emergence of drones and the cessation of building wax combs (Smith et al., 2017). Swarming preparations require coordinated behavior between thousands of individuals. Studies suggest that there is no nepotistic behavior during swarming (Châline et al., 2005; Rangel et al., 2009). Therefore, an environmental factor is probably responsible for triggering such preparations. Many factors, such as increased colony size, changes in the age proportion in the colony toward younger bees, limitations of available brood cells, and reduced transmission of queen pheromones, have been hypothesized to push the colony toward preparations to swarm (Winston, 1991). However, Grozinger et al. (2014) pointed out that these factors are probably more synergistic, as none alone triggers swarming behavior. Indeed, a modeling approach (Fefferman and Starks, 2006) proved that three factors (colony/population size, brood nest congestion, and skewed worker age distribution) generated swarming patterns similar to those obtained in empirical studies. Nonetheless, each model only resulted in swarming after other variables, defined by their individual models, reached their own thresholds. These results may suggest that a yet undefined critical cue triggers swarming and correlates with the aforementioned factors.
One such factor is royal jelly, known to significantly influence shifts in honey bee biology and behavior. Royal jelly is a secretion of nurse worker bees' hypopharyngeal and mandibular glands (Hassanyar et al., 2023) and is best known for regulating cast determination in honey bees. However, there is still ongoing scientific debate surrounding which compound of royal jelly is responsible for the complex physiological and behavioral changes that distinguish workers from queen bees (Maleszka, 2018). Nonetheless, studies have shown that royal jelly can change metabolic flux (Foret et al., 2012), hormone levels (Hartfelder and Engels, 1998) or launch epigenetic cascades leading to alterations in DNA methylation, histone modifications, and non-protein coding RNAs resulting in global changes in gene regulation (Ashby et al., 2016; Barchuk et al., 2007; Dickman et al., 2013). Interestingly, although royal jelly is nutritionally rich with a high proportion of proteins, carbohydrates, fatty acids, vitamins and minerals (Yu et al., 2023), its effect on honey bees' gut biochemical environment and microbiota remains unexplored. However, studies are showing the impact of the application of royal jelly on mouse health and gut microbiota (Chi et al., 2021). Here, we hypothesise that the overproduction of royal jelly by young nursing bees leads to royal jelly cross-feeding among young bees, influencing the chemical composition of their gut and, as a result, promoting or inhibiting the growth of particular bacterial symbionts. By utilizing the method for simultaneous microbiota metabarcoding and quantification using amplicon sequencing of two hypervariable regions of the bacterial 16S marker gene, we aimed to track changes in young worker bees' microbiomes throughout swarming preparation to propose a bacterial symbiont responsible for shifting honey bee behavior toward swarming.
Specimens used in this study came from three experimental colonies in which queens were artificially inseminated by drones originating from a single queen. In early May 2020, after the first drone brood observation, freshly emerged honey bee workers were marked every three days (Suppl. Table 1, Fig. 1). After ten days, worker bees from each batch (sampling time point) were collected from each experimental hive, immediately placed on dry ice, and put in −60°C until processed. Next, brains were collected for RNA and ATAC-seq analyses, whereas the rest of the body was stored in 96% ethanol and −20°C for further DNA extraction. Finally, 78 specimens from three colonies were chosen for microbiome analysis (Suppl. Table 1).
The experimental workflow of the study
The worker bees' abdomens were homogenized by placing them in 2 ml tubes containing 200 μl of a buffer mixture, which included 195 μl of ‘Vesterinen’ lysis buffer (0.4 M NaCl, 10 mM Tris-HCl, 2 mM EDTA pH 8, 2% SDS) (Aljanabi and Martinez, 1997; Vesterinen et al., 2016) and 5 μl of proteinase K. Additional buffer was added if necessary to ensure the samples were fully submerged. Ceramic beads (2.8 mm and 0.5 mm) were added to each tube, and the samples were homogenized using the Omni Bead Raptor Elite homogenizer for two 30-second cycles at a speed of 5 m/s. The samples were then incubated at 55°C for 2 hours in a thermal block.
After cooling, 40 μl of the homogenate from each tube was transferred to a deep-well plate. Based on previous experiments, 100,000 copies of the quantification spike-in were added to each sample — a linearised plasmid containing an artificial 16S rRNA target Ec5001 (Tourlousse et al., 2016) in 2 μl of TE buffer. DNA was purified using 80 μl of SPRI beads on a magnetic stand, followed by two washes with 80% ethanol. The DNA was then diluted with 20.5 μl of TE buffer, and 20 μl of this solution was transferred to a new 96-well plate for DNA concentration of the subset of samples using the Quant-iT PicoGreen kit.
Amplicon libraries were prepared using a custom two-step PCR protocol (Buczek et al., 2024). In the first step, template-specific primers with Illumina adaptor tails were used: 515F and 806R for the amplification of the V4 region (Parada et al., 2016) and 27F and 338R for the V1-V2 region (Walker et al., 2015) of the bacterial 16S rRNA gene. The PCR solution consisted of 5 μl of QIAGEN Multiplex Master Mix, a mix of primers at concentrations 5 μM (16S V1–V2) and 10 μM (16S V4), 2 μl of DNA template, and 1 μl of water (final volume: 10 μl). The temperature program for the first round of PCR included the initial step of denaturation at 95°C for 15 minutes, followed by 25 cycles of denaturation (30s, 94°C), annealing (90s, 50°C) and extension (90s, 72°C) phases, and the final extension step (10m, 72°C). The products were checked on a 2.5% agarose gel against positive and negative controls and cleaned with SPRI beads.
Illumina adapters and unique index pairs were added to the samples during the second indexing PCR. The temperature program for PCR in this step remained the same, but the number of cycles was reduced to 7. In each laboratory step (DNA extraction, first and indexing PCR), we added a negative control (blanks).
The libraries were pooled approximately equimolarly based on band intensity on agarose gels to ensure a roughly equal representation of each sample in the pool. After the last cleaning step with SPRI beads, the pool was ready for sequencing performed on an Illumina MiSeq v3 lane (2 × 250 bp reads) at the Institute of Environmental Sciences of Jagiellonian University.
The obtained data was analyzed by a set of custom Python scripts developed in the Symbiosis Evolution Research Group at Jagiellonian University (Buczek et al., 2024).
First, R1 and R2 files for each library were split into bins corresponding to the target marker region, and primers were cut off. Due to the unique set of informative indexes used for each sample, the potential cross-talk (Wright and Vetsigian, 2016) between libraries was removed, leaving only the reads characteristic of a particular sample. Next, R1 and R2 reads were joined into high-quality contigs (minimum Phred score 30) using PEAR (Zhang et al., 2014). Subsequently, contigs were de-replicated (Rognes et al., 2016) and denoised (Edgar, 2016) separately for each of the libraries to avoid losing biologically relevant information about rare genotypes that could happen during the denoising of the whole dataset at once (Prodan et al., 2020). Chimeras were recognized using USEARCH, and each sequence was taxonomically classified using the SINTAX algorithm and customized SILVA database (version 138 SSU) (Quast et al., 2013). Afterwards, sequences were grouped based on a 97% similarity threshold into Operational Taxonomic Units (OTUs) with the UPARSE-OTU algorithm implemented in USEARCH. At this stage of the analysis, two tables were produced: ASVs table (Amplicon Sequencing Variant) (also known as zOTUs - zero-radius Operational Taxonomic Units), describing genotypic diversity and OTUs table (Operational Taxonomic Units) – clustering genotypes based on the aforementioned similarity threshold.
Using negative controls generated in each laboratory step, bacterial 16S rRNA gene data were screened for putative DNA extraction and PCR reagent contaminants. We first removed genotypes classified as chloroplasts, mitochondria, Archaea, or chimeras using information about taxonomy classification. Next, we calculated relative abundances and used ratios of each genotype presented in blank and experimental libraries to accurately assign genotypes as putative real insect-associated microbes or PCR or extraction contaminants.
Next, reads identified as quantitative spike-ins were used to reconstruct bacterial absolute abundances in the processed honey bee workers. Specifically, the symbiont-to-extraction spike-in ratio, multiplied by the number of extraction spike-in copies and the proportion of the homogenate, allowed us to estimate amplifiable bacterial 16S rRNA copy numbers in the homogenized specimens (Buczek et al., 2024).
Finally, manual analysis was conducted to remove controls and samples with incorrect indexes or zero abundance of bacteria and create the dataset used in the statistical analysis.
Statistical analysis was performed using RStudio version 2023.03.1+446 (R Core Team, 2023), and Processing 3 software version 3.5.4 (Reas and Fry, 2006) was used for generating a heatmap. Inkscape 1.2.2 (Inkscape Project, 2022) was used to modify generated plots and visualizations. All pictures used for the methodological scheme (Fig. 1) come from https://bioicons.com/ website.
To evaluate the impact of batch effects on the absolute abundances of microbial communities, we used a linear mixed effects model (LME) with log-transformed absolute abundance as a response variable, batch as a fixed effect, and hive as a random effect. The model was fitted using restricted maximum likelihood (REML) estimation. The analysis used the lme4 package (Bates et al., 2015). The model's fit was assessed through the REML criterion at convergence, and the significance of fixed effects was evaluated using Satterthwaite's method for degrees of freedom, as implemented in the lmerTest package (Kuznetsova et al., 2017).
The Bray-Curtis (Sørensen, 1948) dissimilarity matrix was computed from the absolute abundance data. This distance matrix quantifies the compositional dissimilarity between pairs of samples based on their zOTU profiles. To assess the influence of temporal changes (represented by different batches) on microbial community composition, a PERMANOVA (Anderson, 2001) was performed using the adonis function from the vegan (Oksanen et al., 2024) R package. We implemented differential abundance analysis to verify the changes of particular bacterial taxa absolute abundance throughout swarming preparation, which was performed using the DESeq2 R package (Love et al., 2014). The DESeq2 package was used to fit a negative binomial model to the count data and to perform Wald tests (Gourieroux et al., 1982) for differential abundances between each batch.
A Bland-Altman plot (Dewitte et al., 2002) was generated using the ggplot2 R package (Wickham, 2016) to assess the agreement between two quantification methods (based on two gene fragments: 16SV1V2 and 16SV4). The mean and difference between the log-transformed values of the two methods were calculated for each sample. The mean difference, known as bias, was computed, along with the 95% limits of agreement, defined as the mean difference ± 1.96 times the standard deviation of the differences.
Using high-throughput 16S rRNA sequencing targeting both V4 and V1-V2 regions, we successfully reconstructed the core gut microbiota of Apis mellifera worker bees. After quality filtering and decontamination, 76 libraries were retained in each dataset. Contaminants and non-bacterial reads accounted for 5.8%–9.1% of the data and were removed prior to analysis (Suppl. Tables 2.1–2.3, 3.1). After the decontamination, the number of reads in the experimental libraries ranged between 23,788 and 161,303, providing sufficient depth for further analysis (Suppl. Table 2.2).
Both datasets consistently recovered all key members of the core microbiota, including Snodgrassella, Gilliamella, Lactobacillus, Bifidobacterium, Frischella, Bombella, Bartonella, Fructobacillus, and Acetobacteraceae. Among these, Snodgrassella, Gilliamella, and Bifidobacterium were dominant, often exceeding 30–70% relative abundance in individual samples (Fig. 2). While taxonomic composition was broadly similar between the two datasets, the V1-V2 region provided higher genotypic resolution, yielding 869 zOTUs compared to 205 from the V4 region. For example, Gilliamella alone was represented by 166 zOTUs (67 exceeding 1% abundance), while Bifidobacterium included over 250 zOTUs across several OTUs. In both datasets, multiple zOTUs from these core taxa reached ≥1% relative abundance, highlighting substantial intra-species diversity (Suppl. Tables 2.3–2.5, 3.2–3.5). Low-abundance but previously described symbionts, including Bombella, Fructobacillus, Pantoea, Erwiniaceae, Enterobacteriaceae, and Spiroplasma melliferum were also detected, usually below 1% relative abundance.
The bacterial relative abundance across the studied batches. A. Reconstruction based on bacterial 16S rRNA hypervariable region V4. B. Reconstruction based on bacterial 16S rRNA hypervariable region V1-V2
We analysed absolute and relative microbial profiles across experimental batches using the V4 and V1-V2 16S rRNA datasets to investigate temporal changes in the honey bee gut microbiota during the swarming preparation (Fig. 3).
Comparison of the absolute abundance of bacteria between batches and used marker gene regions
Absolute microbial abundance, assessed via spike-in controls and log10-transformed, remained relatively stable over time in both datasets. The V4 dataset exhibited values ranging from 7.06 to 9.19 (mean = 8.22), while the V1-V2 dataset ranged from 7.05 to 9.05 (mean = 8.19) (Suppl. Tables 2.7, 3.7; Fig. 3). Linear mixed model analysis revealed no significant batch-related effect on total microbial load (p = 0.735 for V4; p = 0.281 for V1-V2; Suppl. Tables 2.8, 3.8), indicating that absolute quantities of microbiota remained consistent throughout the sampling period.
In contrast to the stability observed in total microbial abundance, the composition of the microbiota displayed significant temporal variation during the swarming preparation period. PERMANOVA analyses revealed that batch identity significantly affected microbial community structure at both OTU and zOTU levels in both datasets. In the V4 dataset, batch explained 16.6% of the variation in OTU-level composition (F = 1.672, R2 = 0.1664, p = 0.025) and 15.3% at the zOTU level (F = 1.5112, R2 = 0.1529, p = 0.013). Similarly, in the V1–V2 dataset, batch accounted for 18.2% of the variance at the OTU level (F = 1.8411, R2 = 0.1824, p = 0.006) and 13.2% at the zOTU level (F = 1.2528, R2 = 0.1318, p = 0.009) (Suppl. Tables 2.9, 3.9). These findings indicate that, while the overall microbial load remained consistent across time points, the relative abundance and genotypic composition of the community shifted in a temporally structured manner. Taking into consideration the overall stability of honey bees' microbiome, both at the level of host species and caste, this variation indicates a strong dynamic of symbiont composition. Nonetheless, further studies based on a more informative approach (such as metagenomics and transcriptomics) would be needed to fully describe potential functional changes.
Differential abundance analyses further supported these compositional dynamics. At the OTU level, only a small number of taxa exhibited consistent changes across batches: Franconibacter (OTU17) and Tyzzerella (OTU31) in the V4 dataset, and Bombella (OTU29) in the V1-V2 dataset showed significant shifts in relative abundance across at least four consecutive batches (Suppl. Tables 2.11, 3.11). However, at the finer genotypic resolution of zOTUs, more pronounced temporal variability was observed, particularly in the V1-V2 dataset. Out of 193 zOTUs showing significant batch-related changes, 23 exhibited consistent fluctuations across five or more consecutive batches. These dynamic genotypes belonged primarily to core taxa such as Gilliamella (OTUs 1, 16, and 30), Bifidobacterium (OTUs 2, 8, and 13), Snodgrassella (OTU3), Bartonella (OTU7), and Frischella (OTU9). For instance, seven zOTUs from Gilliamella OTU1 varied across five to seven batches, as did multiple genotypes from Bifidobacterium OTU8. The relative abundance of these temporally responsive zOTUs ranged from 0.19% to 34.07% in individual samples, with read counts between 233 and 35,712 across the dataset (Suppl. Table 3.10).
In contrast, the V4 dataset showed minimal genotypic turnover, with only a single zOTU of Bombella (zOTU16) exhibiting significant shifts across multiple batches, specifically between batches three and five (Suppl. Table 2.10). Altogether, these results indicate that although the overall microbial load remains stable, the community composition, particularly at the genotypic level, undergoes subtle, taxon-specific changes throughout the swarming preparation. These shifts, most evident in taxa such as Gilliamella, Bifidobacterium, and Frischella, likely reflect ongoing ecological restructuring or strain-level turnover within the gut microbiota of worker bees.
Insect species differ dramatically in the abundance of microorganisms they host, and these differences often correlate with the microbes' function in insect biology (Hammer et al., 2019). In the case of healthy honey bees, the estimates suggest a stable and constant abundance at a level between 108 and 109 (Motta and Moran, 2024). However, most studies fail to report bacterial absolute abundance with a few notable exceptions (e.g., Schwarz et al., 2016; Kešnerová et al., 2020; Castelli et al., 2022), leaving a gap in understanding the quantitative dynamics of bee-associated microbiota.
Here, we used a new approach to quantify the microbial community as additional information without increasing the costs of multitarget amplicon sequencing (Buczek et al., 2024), which opens up new possibilities to study insect-microbiome interactions on a new higher level. The minimum and maximum absolute abundance estimates were congruent between the V1–V2 and the V4 datasets, ranging between 107 and 109. This consistency supports the accuracy of our method, and it also shifts the described abundance (Motta and Moran, 2024) to lower values. The lack of changes in overall microbiome abundance with simultaneous changes in microbial composition might suggest that physiological factors allow honeybees to sustain bacterial quantities on a particular stable level during the swarming preparation. We can not also exclude the possibility that the stability in microbiome quantities among the analyzed bees is related to their age. Specifically, since we used 10-day-old worker bees throughout swarming preparation, the overall homogeneity of the abundances might be the outcome of the relatively short time from their emergence from pupae, as the adult honey bees emerge sterile and acquire symbionts through trophallaxis with other nestmates (Moran, 2015). As the detected minimal absolute abundance (107) was one order of magnitude lower than previously reported (108), it would be possible that the quantities of relatively young bee workers do not reach their optimal abundances yet and may increase in the following days of their life. On the other hand, Kešnerová et al. (2020) demonstrated the seasonal variability in bee microbiota abundance, revealing statistically significant differences in total microbiome load between winter bees and foragers. Notably, our results correspond well with their findings, showing that foragers sampled in May and June (the same month we sampled our bees) tend to have lower abundances (107). Thus, it is possible that the microbial quantities may be influenced not only by the bee's age (behavioral tasks) but also by season-dependent environmental factors such as pollen source and quality.
The microbial community of honey bees is relatively consistent in terms of composition; however, beyond the core microbiota, these insects may host additional bacteria whose presence may be host-specific and related to various factors. In this study, we analyzed the honey bee microbiome using the V1-V2 and V4 regions of the bacterial 16S rRNA gene, revealing slight differences in microbiome composition depending on the region used. Our results confirmed previous findings showing that relying on a single hypervariable region of 16S rRNA is not sufficient for accurate reconstruction of the honey bee microbiome, with the V1-V2 region complementing the widely used V4 region (Romero et al., 2019). Nonetheless, both datasets consistently identified five core clusters characteristic of the microbiome of honey bees: Snodgrassella, Bifidobacterium, Lactobacillus (currently distinguished into Lactobacillus and Bombilactobacillus), and Gilamella (Luo et al., 2024; Motta and Moran, 2024). Additionally, we captured honey bee symbionts that are not included in the core microbiome: Frischella, Bartonella, Bombella, Fructobacillus, and bacterium belonging to the Acetobacteraceae group (Motta and Moran, 2024). However, we observed a clear bias in terms of OTU diversity. The V4 primers distinguished more OTUs within the Lactobacillus group, whereas V1-V2 primers were better at detecting Snodgrassella, Gilamella, Bifidobacterium, Bartonella, and Bombella. Surprisingly, we did not find Commensalibacter, a taxon typically present in honey bee microbiota (Motta and Moran, 2024). However, this might be explained by results indicating that this taxon is much more prevalent in winter bees than spring/summer bees (Kešnerová et al., 2020). Other bacteria often described as parts of the honey bee microbiome, including Pantoea, Serratia (only in the V4 dataset), and Enterobacteriaceae (Kešnerová et al., 2020; Motta and Moran, 2024) were also found in our datasets.
Previous studies have shown that microbiome changes in honey bees are not dependent on the location of the apiary but rather on the time point at which bees were collected (Almeida et al., 2023). Our results support that claim, showing a significant effect of the batch on microbial composition with no effect of the hive factor. Results of our analyses indicate that, depending on the dataset (V1-V2 vs V4) and clustering level (zOTUs vs OTU), different bacterial taxa or genotypes can be distinguished with their absolute abundance shifting through the swarming preparation. In the V4 dataset, we identified two OTUs with statistically relevant changes: Franconibacter and Tyzzerella. The former one was described as a member of a saccharification agent used to initiate fermentation in the production of Chinese liquor and vinegar (Gao et al., 2017), whereas the second was recently proposed as a factor playing a potential role in larval biology (Maigoro et al., 2024). Hence, although they might be members of the honey bees' microbiota, considering their low abundance and the limitations of amplicon-based bacterial metabarcoding, we would be far from assessing their significant role in triggering the swarming preparation or playing an essential role in the process. On the other hand, the changes in the prevalence of Bombella zOTU might indicate its function in this process. Bombella has been recognized as a diverse taxon with different strains contributing to various aspects of honey bees' larval and adult biology, ranging from protection against fungal pathogens, to larvae supplementation in lysine and buffering nutritional stress (Miller, 2023; Parish et al., 2022). Lysine is directly involved in nitric oxide synthesis, a neurotransmitter affecting honey bees' brain functions (Gage et al., 2020). Notably, Bombella is primarily found in the royal jelly (Corby-Harris et al., 2014), further supporting our hypothesis that the swarming behavior might be directly or indirectly triggered by an overconsumption of this compound by young bees. Hence, further studies involving the mix of omics techniques and proper experiments should be conducted to fully describe the potential role of Bombella in the swarming process.
The V1-V2 dataset analysis showed much more variation on both zOTU and OTU levels. In the zOTU dataset, we found 23 genotypes with significant changes in at least four consecutive batches. Those genotypes belonged to Gilliamella, Bifidobacterium, Snodgrassella, and Bartonella groups. Considering that the V4 dataset showed a stable diversity pattern of the first three species (belonging to the core microbiome), we tend to interpret those changes as an outcome of a bias caused by V1-V2 primers to overrepresent the genetic diversity of those taxa, as shown in other studies comparing microbiomes based on those two gene regions (Åhlén Mulio et al., 2024). Naturally, the observed differences in zOTU composition could be explained by biological factors, such as changes in the abundance of particular bacterial strains during swarming preparation. However, we still lack evidence of whether such mechanisms occur in individual bees (Ellegaard and Engel, 2019). Bartonella was a bacterial symbiont captured in the V1-V2 dataset but absent in the V4 dataset. This bacterium has been described as especially abundant in workers involved in hive tasks such as food processing (Jones et al., 2018) and requiring food rich in pollen to grow in abundance successfully (Kešnerová et al., 2020). Its change over the sampling period might be explained by a higher pollen supply available from the environment between May and June. The genomic-based reconstruction of Bartonella metabolism showed that it is capable of tryptophan and phenylalanine secretion (Li et al., 2022). Tryptophan regulation by Lactobacillus has been shown to influence the memory behavior and learning of worker bees (Zhang et al., 2022a), suggesting its potential role in triggering the changes in behavior during swarming preparations. Interestingly, the only zOTU showing statistically significant changes in at least four consecutive batches was Bombella, which overlaps with results obtained by analyzing the OTU V4 dataset.
Although our results show only subtle changes in microbiome composition during the swarming preparation, the detected shifts may still have important functional consequences for bee physiology and behavior. The process of swarming involves complex social coordination, increased responsiveness to queen pheromones, and altered nutritional and metabolic demands (Seeley, 2010). Recent studies have demonstrated that gut bacteria can influence host neurochemistry and signaling molecules, including those derived from amino acid metabolism such as tryptophan, phenylalanine, and lysine (Zhang et al., 2022; Gage et al., 2020; Li et al., 2022). In this context, the observed increase in Bombella abundance could modulate lysine availability and downstream nitric oxide synthesis, potentially affecting neuronal activity and decision-making processes associated with swarming initiation. Likewise, changes in Bartonella, a taxon capable of producing aromatic amino acids, may impact the regulation of serotonin and dopamine pathways, thereby influencing motivation and task allocation within the colony.
Beyond neurotransmitter precursors, altered microbiome composition may also affect nutritional processing and immunity, indirectly shaping the bees' readiness to swarm. Bombella and Bartonella have both been linked to enhanced resilience against nutritional stress and pathogens (Corby-Harris et al., 2014; Kešnerová et al., 2020; Hamilton et al., 2021; Miller et al., 2023), traits that may be beneficial during the energetically demanding pre-swarming phase. Therefore, even minor compositional adjustments in these bacteria might represent adaptive fine-tuning of the colony's physiological state, facilitating the transition from nursing to swarming. Future studies combining metagenomics, metabolomics, and behavioral assays will be crucial to determine whether these microbial shifts contribute causally to the behavioral changes characteristic of swarm preparation.
Honey bees' taxonomically stable microbiota, with its potential strain diversity combined with their complex behavior, makes them the perfect model species for studying gut-brain interactions. We tested the hypothesis that the gut-brain axis can trigger swarming preparation by implementing an approach for simultaneous microbiota metabarcoding and its quantification based on artificial 16S plasmids. Our results indicate that it is plausible that some taxa may play a role in changing the behavior of young honey bees. Especially Bombella holds significant promise, as it was shown to vary over time and is known to be involved in lysine supplementation, which is a neurotransmitter affecting bee brain functions. The primary source of lysine is royal jelly, which is more common in colonies preparing to swarm. Therefore, it is possible that Bombella is involved in an important taxon in the process of swarming preparation. Although further studies involving gene expression analysis and genome reconstruction of a symbiont strain would be required, similar to an approach used in ants showing the correlation between microbiome composition and brain gene expression (Kay et al., 2023).