Composition and diversity of bacterial communities associated with honey bee foragers from two contrasting environments

Abstract

The honey bee is associated with a diverse community of microbes (viruses, bacteria, fungi, and protists), commonly known as the microbiome. Here, we present data on honey bee microbiota from two localities having different surrounding landscapes – mountain (the Rhodope Mountains) and lowland (the Danube plain). The bacterial communities of abdomen of adult bees were studied using amplicon sequencing of the 16S rRNA gene. The composition and dominance structure and their variability within and between localities, alpha and beta diversity, and core and differential taxa were compared at different hierarchical levels (operational taxonomic units to phylum). Seven genera (Lactobacillus, Gilliamella, Bifidobacterium, Commensalibacter, Bartonella, Snodgrassella, and Frischella), known to include core gut-associated phylotypes or species clusters, dominated (92–100%) the bacterial assemblages. Significant variations were found in taxa distribution across both geographical regions and within each apiary. Lactobacillus (Firmicutes) prevailed significantly in the mountain locality followed by Gilliamella and Bartonella (Proteobacteria). Bacteria of four genera, core (Bartonella and Lactobacillus) and non-core (Pseudomonas and Morganella), dominated the bee-associated assemblages of the Danube plain locality. Several ubiquitous bacterial genera (e.g., Klebsiella, Serratia, and Providencia), some species known also as potential and opportunistic bee pathogens, had been found in the lowland locality. Beta diversity analyses confirmed the observed differences in the bacterial communities from both localities. The occurrence of non-core taxa contributes substantially to higher microbial richness and diversity in bees from the Danube plain locality. We assume that the observed differences in the microbiota of honey bees from both apiaries are due to a combination of factors specific for each region. The surrounding landscape features of both localities and related vegetation, anthropogenic impact and land use intensity, the beekeeping management practices, and bee health status might all contribute to observed differences in bee microbiota traits.

Introduction

The western honey bee (Apis mellifera) is a key pollinator of many agricultural crops and wild plants (Hung et al., 2018). The ongoing trend of honey bee decline throughout the world has brought the attention to the great loss of ecological services provided by this species, its health status, and numerous threats it faces (Paudel et al., 2015; Boncristiani et al., 2020). Bees are associated with a diverse community of microbes (viruses, bacteria, fungi, and protists), commonly known as the microbiome. As a superorganism, the honey bee has complex interactions with its own microorganisms, the microbial communities of food stores, hive surfaces, and surrounding environments; therefore, four major groups of bee-associated microbiomes have been distinguished (Smutin et al., 2022). These include gut and body surface microbiota, microbiota of the in-hive environment (honeycombs, brood combs, propolis, royal jelly, etc.), and microbiota of the surrounding environment. Recent advances in DNA technologies and biomolecular analyses enabled new insights into the understanding of microbiota composition, dynamics, and impact on host biology and health. Except for a few bacterial pathogens that are widespread in apiaries and causing high mortality (Forsgren, 2010; Genersch, 2010), the functional role of most members of bee-associated microbiota on host and colony health is not completely understood (Engel et al., 2016; Khan et al., 2020). The gut microbiome is the most extensively studied part of the bee-associated microbiota (Romero et al., 2019) used as a model for gut microbiota research (Zheng et al., 2018). It is relatively simple and specialized, consisting of 8–10 recurring phylotypes or species clusters called Alpha-1 (Bartonella apis), Alpha-2.1 (unclassified species of Acetobacteriaceae family), Alpha-2.2 (Saccharibacter floricola), Beta (Snodgrassella alvi), Bifido (Bifidobacterium asteroides), Gamma-1 (Gilliamella apicola), Gamma-2 (Frischella perrara), Firm-4 (Lactobacillus mellis and L. mellifer), Firm-5 (other Lactobacillus spp. including (Apilactobacillus) kunkeei) (Sabree et al., 2012; Kwong and Moran, 2016). Despite the conserved bacterial phylotype composition, it was found that each honey bee harbors a distinct community at the functional level (Ellegaard and Engel, 2019).

It is well known that S. alvi, G. apicola, two species of Lactobacillus (L. apis and L. helsingborgensis), and Bifidobacterium spp. are ubiquitous and can be found in every adult honey bee worker worldwide, often considered as part of the core gut microbiome composition (Erban et al., 2017; Nowak et al., 2021). Other microbial species, e.g., B. apis, Apibacter adventoris, F. perrara, and species of the family Acetobacteraceae, have been sometimes found in the gut of many workers, but exceptions exist (Jones et al., 2018; Raymann and Moran, 2018). The gut bacteria showed specific spatial distribution related to different sections of the honey bee digestive system (Martinson et al., 2012; Callegari et al., 2021). Snodgrassella alvi and G. apicola dominate in the ileum region of the hindgut, F. perrara occupy the pylorus near the beginning of the ileum, while the Lactobacillus strains and Bifidobacterium spp. are most abundant in the rectum region of the hindgut (Raymann and Moran, 2018).

The functional role of the gut microbiome was mainly related to host nutrition and metabolism, regulation of the immune system, and pathogen resistance (Emery et al., 2017; Kwong et al., 2017; Bonilla-Rosso and Engel, 2018; Daisley et al., 2020; Castelli et al., 2021), and alterations in gut microbiota (such as dysbiosis) may affect host metabolic homeostasis and immunity (Yuan et al., 2019; Zhu et al., 2020). Multiple factors, such as bee age, caste and social activity, presence of pathogens and parasites, exposure to agrochemicals and anthropogenic toxins, the microbiota of in-hive environment and stored food, as well as the landscape heterogeneity and related local environmental conditions (including vegetation diversity, land-use type, season), are known to affect honey bee gut microbiome composition and dynamics (Anderson et al., 2013; Graystock et al., 2015; Engel et al., 2016; Kwong and Moran, 2016; Jones et al., 2018; Bosmans et al., 2018; Ribière et al., 2019; Kešnerová et al., 2020; Muñoz-Colmenero et al., 2020). The versatility and potential utility of the western honey bee, including its associated microbiota, as biomonitors of environmental pollutants, pesticides use, presence of pathogens, climate change, and ecosystem health, has been recognized and recently reviewed (Cunningham et al., 2022).

Traditionally, the microorganisms associated with honey bees and their food are determined by examination, and subsequent isolation and identification, using appropriate tests and taxonomic keys (see Snowdon and Cliver, 1996; Gilliam, 1997 for reviews). Very often taxonomic identification is difficult to implement in species with similar morphological features. This requires an application of high-sensitivity methods such as DNA sequencing including DNA barcoding (Porter and Hajibabaei, 2018; Eberle et al., 2020). Recently, a novel method (environmental DNA metabarcoding) comprising DNA isolation from mixed or bulk environmental samples, amplification, high-throughput sequencing (HTS), and bioinformatics was widely applied for studying and monitoring biological diversity including microbiome diversity (Deiner et al., 2017; Piper et al., 2019; Kennedy et al., 2020). HTS of 16S and internal transcribed spacer ribosomal RNA (rRNA) markers is now one of the most effective and extensively used approaches in microbiome studies that also provides quantitative estimates of the microbial community composition (Bukin et al., 2019).

In this study, we investigated the bacterial communities associated with honey bees from two regions in Bulgaria (mountain and lowland) by applying DNA metabarcoding approach based on 16S rRNA gene. The objectives were: (i) to detect ‘core’ bee-associated bacteria and differential taxa specific to each locality and (ii) to compare the bacterial community structure and diversity within and between localities.

Materials and methods

Sites description and sampling

Bee specimens were collected during honey production season (June 2018) from two different localities: the ‘Experimental Apiary of the Research Center of Stockbreeding and Agriculture’, located on the outskirts of Smolyan city (41°35′40″N, 24°40′39″E), the Rhodope Mountains (RM), Southern Bulgaria and from a private apiary in a lowland village near to Byala city (43°25′6″N, 25°49′14″E), Danube plain (DP) region, Northeastern Bulgaria. Both sites differed considerably by their topography and landscape characteristics, altitude, dominant land use type, and vegetation. The locality in the Rhodope Mtn is mountainous with an average altitude of about 1160 m, a very heterogeneous environment, steep slopes, more diverse landscape structures, numerous grasslands and pastures, and coniferous tree vegetation. In contrast, the landscape of the lowland locality in the Danube plain region is slightly hilly (altitude ranges 100–300 m a.s.l.), surrounded by agroecosystems with annual and perennial crops, broadleaved vegetation, small forest patches.

To identify bee-associated microbiota, 15 hives from each apiary were randomly selected. From each hive, one bee forager returning to the hive with pollen on their legs was collected with sterile tweezers from each hive near to hive entrance. Each bee specimen was placed in a sterile 50 μl falcon and immediately frozen on dry ice in the field. The samples were brought to the laboratory of the Institute of Biodiversity and Ecosystem Research and stored at −20°C before processing.

DNA extraction, PCR amplification, and sequencing

Prior to DNA extraction, the posterior body segment of bee specimens (known also as metasoma or abdomen) was cut-off with sterile scissors. For pragmatic reasons, we will use the term abdomen in this text. After that, the bees were mechanically homogenized by sterile pestle in a total of 350 μl Lyse T buffer (GeneMATRIX Tissue DNA Purification kit, Cat.no. E3550, EURx Sp. z o.o.; Gdansk, Poland) to complete the destruction of the abdomen. After removal of all tissue fragments, total DNA isolation was performed according to the manufacturer's instructions. DNA concentration and quality were measured using an Analytikjena spectrophotometer (Analytik Jena GmbH, Germany). The extracted DNA was then sent to the AIM laboratory (https://www.aimethods-lab.com/) for DNA metabarcoding using NGS (sequencing and bioinformatics processing).

From each sample, 5 μl of extracted genomic DNA was used, along with Plant MyTAQ (Bioline, Luckenwalde, Germany). The V3–V4 hypervariable regions of the bacterial 16S rRNA gene were amplified with the adapted primers 16s-barcode primers (341f – TACACGACGCTCTTCCGATCTTCATCCTACGGGNGGCWGCAG – and 785r – CAGACGTGTGCTCTTCCGATCCGCTCAGACTACHVGGGTATCTAATCC) (Morinière et al., 2016; Thijs et al., 2017; Morinière et al., 2019). Amplification success and fragment length were assessed using gel electrophoresis. Amplified DNA was cleaned up and resuspended in 50 μl molecular water for each sample before proceeding. Illumina Nextera XT (Illumina Inc., San Diego, USA) indices were ligated to the samples in a second PCR reaction applying the same annealing temperature as for the first PCR reaction but with only seven cycles, and ligation success was confirmed by gel electrophoresis. DNA concentrations were measured using a Qubit fluorometer (Life Technologies, Carlsbad, USA), and samples were combined into 40 μl pools containing equimolar concentrations of 100 ng each. Pools were purified using MagSi-NGSprep Plus (Steinbrenner Laborsysteme GmbH) beads. A final elution volume of 20 μl was used. HTS was performed on an Illumina MiSeq using v3 (2 × 300 bp, 600 cycles, maximum of 25 mio paired-end reads) chemistry.

Bioinformatic processing

Paired-end merging was done using the -fastq_mergepairs utility of the USEARCH suite v11.0.667_i86linux32 (Edgar et al., 2011). Adapter sequences were removed using CUTADAPT (Martin, 2011) (default parameters). The remaining pre-processing steps, namely quality filtering, dereplication, chimera filtering, and clustering were carried out using the VSEARCH suite v2.9.1 (Rognes et al., 2016). Sequences were trimmed to a length of 400 bp and were only included in the analyses that matched perfectly to the barcode. 16S rRNA genes were classified using the SILVA database (Quast et al., 2012) and NCBI database (O'Leary et al., 2016) both with a minimum similarity threshold of 97%. The resulting file could then be used to map the reads to the operational taxonomic units (OTUs) to create the OTU table. To reduce likely false positives, a cleaning step was employed which excluded read counts in the OTU table that constituted less than 0.01% of the total number of reads in the sample. OTUs were blasted against a custom database downloaded from GENBANK using Geneious (v.10.2.5 – Biomatters, Auckland – New Zealand), and following methods described in Morinière et al. (2016). OTUs were additionally removed from the results based on negative control samples, i.e., if the combined number of reads in the negative controls constituted more than 20% of the total number of reads in the OTU. OTUs were also annotated with the taxonomic information from NCBI (downloaded from https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/), followed by a creation of a taxonomic consensus between BOLD, NCBI, and RDP.

A total of 1,667,218 raw reads of the V3–V4 region of the bacterial 16S rRNA were obtained and assembled in 833,609 clean sequences with an average sequence length of 411 bp, and with a maximum of one expected error. These sequences were obtained from the 20 assayed samples (13 from the Rhodope Mountains and 7 from the Danube plain). Before chimera detection, 419 bacterial OTUs were identified at the 97% sequence similarity cut-off. After chimera filtering and removing data with low-quality scores for three bee specimens from Danube plain locality and OTUs identified as organelle DNA (chloroplast and mitochondria), 66 OTUs and 17 samples remained for further analyses (Supplementary table S1).

Data analyses

The microbiome data analyses were performed using the web-based platform MicrobiomeAnalyst (Dhariwal et al., 2017; Chong et al., 2020) based on the OTUs abundance table for 17 samples containing raw counts of 66 OTUs. Taxonomy data in text table format and the metadata file were uploaded and further analyzed in the Marker-gene Data Profiling (MDP) module. Taxa having less than a 10% prevalence with a read number ≤2 were filtered out. The 39 most frequent and abundant OTUs were further involved in the analyses after data filtering. The OTU abundances were rarefied to the minimum library size (Supplementary fig. S1a) (Kandlikar et al., 2018) and normalized. Rarefaction curves per each sample are presented (Supplementary fig. S1b). Core bacterial taxa were distinguished and visualized by setting up a 20% prevalence value and 0.01 for the relative abundance (RA) value. Three alpha-diversity measures (Shannon diversity, Simpson diversity, and Chao1 indices) were calculated at OTUs and generic levels and visualized (McMurdie and Holmes, 2013). Chao1 index estimates taxa richness by accounting for taxa (e.g., OTUs, species, or genera) that are undetected because of low abundance. Shannon and Simpson's indices are based on taxa richness and their abundances (or evenness). Pairwise comparisons and the significant differences between both groups were calculated using the non-parametric Mann–Whitney U test (P ≤ 0.05). To observe the average dissimilarity in the composition of microbial communities in both regions of the country, beta diversity was estimated based on Bray–Curtis similarity distances and visualized with non-metric multidimensional scaling (NMDS) ordination (McMurdie and Holmes, 2013). The statistical significance of the clustering pattern in ordination plots was evaluated using the permutational multivariate analysis of variance (PERMANOVA). Differences in taxonomic classification between both geographical regions were illustrated using heat tree analysis (Foster et al., 2017). It leverages the hierarchical structure of taxonomic classifications to quantitatively (using the median abundance) and statistically (using the non-parametric Wilcoxon rank-sum test, P < 0.05) depict taxonomic differences between microbial communities (view mode: comparison) or abundance profiles (view mode: abundance). Linear discriminant analysis of effect size (LEfSe) was applied to identify bacterial taxa specific to each locality using a significance level of P < 0.05 (Segata et al., 2011).

Results

Bacterial community structure, composition, and variability

The taxonomic assignment of the 39 most common and abundant OTUs revealed five phyla, six classes, 12 orders, 15 families, and 22 genera (fig. 1, Supplementary fig. S2). Among bacterial phyla, Proteobacteria, Firmicutes, and Actinobacteria dominated the bee-associated communities at both locations (99–100%, fig. 1a). Proteobacteria prevailed significantly in all bee specimens from the lowland locality (DP, 68–87%) and in five bee specimens from the Rhodope Mtn (RM, 52–79%), followed by Firmicutes (RM: 17–63%, DP: 10–27%) and Actinobacteria (0.4–14%). OTUs assigned to Bacteroidetes and Tenericutes were present in a very low RA (0.3–0.9%) and occurrence in both localities. A clear variability within and between localities was also observed in the RAs for bacterial families (Supplementary fig. S2a) and genera (fig. 1b). In the bee specimens from the Rhodope Mtn, on average 43% of reads belonged to the family Lactobacillaceae (mainly Lactobacillus, 43%); Orbaceae make up 14% of Bacteria (Gilliamella, 12% and Frischella, 2%), followed by 10% RA of Rhizobiaceae (Bartonella), 9% – Enterobacteriaceae (mainly Morganella, 8%), 8% – Neisseriaceae (Snodgrassella), and about 7% of Bifidobacteriaceae (Bifidobacterium) and Acetobacteraceae (Commensalibacter). Other families, such as Leuconostocaceae (Fructobacillus), Pseudomonadaceae (Pseudomonas), Spiroplasmataceae (Spiroplasma), and Sphingomonadaceae (Sphingomonas), were found in less than 1% of the total OTU abundance (Supplementary fig. S2a).

Figure 1. Bar charts illustrating the distribution and relative abundances of bacterial phyla (a) and genera (b). Data are represented for each honey bee specimen and locality. RM, Rhodope Mountains; DP, Danube plain.

In bee specimens from the Danube plain locality, on average 30% of OTUs belonged to the family Rhizobiaceae (Bartonella), 21% – Lactobacillaceae (mainly Lactobacillus – 20%), 18% – Enterobacteriaceae (Morganella – 10%, Klebsiella – 6%, Serratia – 1%), 14% – Pseudomonadaceae (Pseudomonas spp.), and 7% to Orbaceae (Gilliamella – 6% and Frischella – 1%). The average RA of bacteria belonging to families Neisseriaceae (Snodgrassella), Bifidobacteriaceae (Bifidobacterium), and Acetobacteraceae (Commensalibacter, Bombella) was lower than 3% (2–4%). Bacteria from five families (Carnobacteriaceae, Flavobacteriaceae, Oxalobacteraceae, Sphingobacteriaceae, and Sphingomonadaceae) were found in low frequency and RAs (>0.3%) (Supplementary fig. S2a). Heatmaps visualizing bacterial taxa composition and their RAs at various taxonomical levels (species, genera, families, and orders) for each locality and individual bee specimens are presented in the Supplementary fig. S3a–d.

Bacterial diversity, core, and differential taxa

The total number of OTUs found in all bees from both localities was 66 (Supplementary table 1) with one-third of them (or 27 OTUs) revealed in a very low number of reads (≤10 per bee specimen) mainly from the Danube plain locality. Lactobacillus and Gilliamella showed the highest OTU richness. Because of the limited number of successfully sequenced specimens from this locality, we could not be sure if these are indeed rare bacterial phylotypes. Despite this bias in bee specimens' quantity, the OTU richness of microbial communities in bees from the Danube plain was higher as compared to those from the Rhodope Mtn (29 vs. 35, in total) and per bee specimen (fig. 2). Twenty-five OTUs and 15 genera were common for microbial communities of both localities, of them ten OTUs and five genera were unique for the Danube plain area (Supplementary fig. S4).

Figure 2. Alpha diversity measure showing the observed variation of OTU numbers (richness) in the bacterial communities. Data are represented as a box plot. RM, Rhodope Mountains; DP, Danube plain; Mann–Whitney U test statistic: 0.05, P value: 0.004.

Alpha diversity assessed by Shannon, Simpson, and Chao1 indices revealed different trends when calculated at OTUs and genus levels (Supplementary fig. S5a–f). The bacterial communities of bees from the Danube plain were significantly more diverse than those of the Rhodope Mtn when genus data were compared (Mann–Whitney U test, Shannon (P = 0.04), Simpson (P = 0.059), and Chao1 (P = 0.01); Supplementary fig. S5d–f). Except for Chao1's diversity index (P = 0.004; Supplementary fig. S5c), no significant differences in OTU diversity between both localities were revealed (P > 0.05; Supplementary fig. S5a, b). Beta diversity analyses clearly distinguished the bacterial communities from both localities (P = 0.002) (NMDS; fig. 3). The microbial communities of bee specimens from the Rhodope Mtn were more similar to each other than those of the Danube plain region.

Figure 3. Non-metric multidimensional scaling (NMDS) ordination of honey bee communities at OTU level using Bray–Curtis as a dissimilarity measure. Results from the PERMANOVA test (F, stress, and P values) are represented for each taxonomic level. Individuals from the same colony and explained variances are represented by the same colors.

Several core bacterial taxa occurring with high frequency and abundance specific to each locality were revealed (fig. 4). Bartonella (B. apis), Bifidobacterium sp., Snodgrassella (S. alvi), Frischella (F. perrara), and phylotypes of Gilliamella, Lactobacillus, and Commensalibacter, known as honey bee gut-associated bacteria, dominated in almost all bee individuals (92–100%). However, the percentage of these specific core gut bacterial genera was less abundant (32–73%) in a few bee specimens. OTUs assigned to other four genera (Pseudomonas, Serratia, Morganella, and Klebsiella) were present with higher RA in bee specimens from the Danube plain and occasionally – in bees from the Rhodope Mtn (Supplementary fig. S3a, b).

Figure 4. Heatmap representing the core bacterial genera dominating the honey bee communities from the Rhodope Mountains (a) and Danube plain (b). The y-axis represents the prevalence of core genera across the detection threshold range (relative abundance) on the x-axis. The prevalence of each core species is indicated by a color gradient from dark blue (low) to dark red (high).

The geographical differences in the bacterial microbiome composition and structure were also supported by the results of two additional analyses – the heat tree analysis (fig. 5) and the LEfSe (Supplementary fig. S6, table S3). Both analyses showed differences in the distribution of seven bacterial genera, i.e., Lactobacillus (LEfSe, P = 0.02), Gilliamella (LEfSe, P = 0.03), Snodgrassella (LEfSe, P = 0.03), Bartonella (LEfSe, P = 0.02), Pseudomonas (LEfSe, P < 0.00), Serratia (LEfSe, P < 0.00), and Apibacter (LEfSe, P = 0.01) with the first three genera being significantly more abundant in the Rhodope Mtn (Supplementary fig. S7). Bacteroidetes was found in honey bees from the Danube plain only (Wilcoxon's P = 0.001) with Apibacter sp. occurring in three out of four specimens (LEfSe, P = 0.01, Supplementary table S3). The order Pseudomonadales, composed of four Pseudomonas spp., prevailed significantly in bee specimens of the Danube plain region (Wilcoxon's P = 0.001) (fig. S6). Other differential taxa distinguishing bee-associated bacterial communities of both localities are presented in the Supplementary materials (fig. S6 and table S3).

Figure 5. Heat tree graph representing the geographical differences in the bacterial communities associated with honey bees. Significantly different taxonomic ranks are shown on the tree as a light red and dark blue nuance for communities of the Rhodope Mountains and the Danube plain, respectively. Node diameter is proportional to the number of OTUs classified as that taxon.

Discussion

In our study, we focused on the composition and diversity of bacterial communities associated with honey bees from two apiaries located in the Rhodope Mtn and Danube plain regions having different landscape topography – mountain and lowland, respectively. We characterized the bee-associated bacteria at different hierarchical levels (OTUs to phylum) by applying various community-based analyses such as composition and dominance structure, and their variability within and between localities, α and β diversity, analyses for revealing core and differential taxa.

Bacterial OTUs assigned to seven genera (Bartonella (B. apis), Bifidobacterium sp., Snodgrassella (S. alvi), Frischella (F. perrara), and Gilliamella, Lactobacillus, and Commensalibacter), known as specific gut-associated bacteria of A. mellifera, dominated (92–100%) the bee-associated assemblages in most specimens. Phylotypes or species clusters from these seven genera are known as specific gut microbiota that is well-adapted to honey bees worldwide, mainly transmitted through social interactions between individuals (Kwong and Moran, 2016). They occupy different sections of the bee digestive system and have evolved specific host–microbe interactions and functions (Martinson et al., 2012). Two ubiquitous Gram-negative genera Snodgrassella (S. alvi) and Gilliamella (members of the Proteobacteria phylum) and several Gram-positive Lactobacillus bacteria (Firmicutes phylum) prevailed significantly in the Rhodope Mtn locality. One Proteobacteria species (B. apis), usually found with lower abundance and frequency, dominated the bacterial communities of bees from the Danube plain locality. These primary bacterial genera are part of the core gut microbiota of honey bees (Kwong and Moran, 2016). The F. perrara usually localized in a restricted niche of the bee gut (Kwong and Moran, 2016), and one unidentified Bifidobacterium species were present in all specimens from both localities with low and relatively high RA, respectively. Except for F. perrara (Engel et al., 2015a, 2015b; Maes et al., 2016), these bee-specific bacterial taxa are considered beneficial gut symbionts of social bees that have protective and nutritional functions (Kwong and Moran, 2016). For example, Lactobacillus spp. has essential roles in carbohydrate metabolism and in protecting hosts by producing antimicrobial metabolites (such as organic acids, diacetyl, benzoate, and bacteriocins) and sugars (Zhang et al., 2022). Bifidobacterium spp. are involved in the biosynthesis and assimilation of trehalose that is used for energy storage in insects and its presence and abundance are considered bacterial indicators for assessing the age of bee workers (Dong et al., 2020).

The honey bee biology and colony behavior, as well as the microbiome composition of the in-hive environment, strongly affect the routes of acquisition and the initial composition of bacteria associated with worker bees (Powell et al., 2014; Smutin et al., 2022). Subsequently, the microbiomes of various in-hive structures (surfaces), stored foods (e.g., honeycombs, bee wax, propolis, pollen, royal jelly, honey, etc.), and out-hive environments (floral nectar and pollen) can contribute to the acquisition of non-core bee-associated microbiota (Snowdon and Cliver, 1996; Kačániová et al., 2009; Anderson et al., 2011; Rering et al., 2018; Cullen et al., 2021). Species of the genus Lactobacillus, the most recognized and diverse honey bee symbionts, were commonly found in the in-hive environment including bee bread, royal jelly, bee pollen, and nectar (Anderson et al., 2013, 2014). For example, the fructophilic lactic acid bacteria L. (Apilactobacillus) kunkeei isolated from all beehive substrates and crop sections of the bee alimentary tract (Corby-Harris et al., 2014) have been found worldwide on different plants, especially associated with floral nectar and pollen (McFrederick et al., 2012; Neveling et al., 2012; Anderson et al., 2013; Rering et al., 2018; Casalone et al., 2020; de Oliveira Scoaris et al., 2021; Rering et al., unpublished).

In some specimens, a substantial deviation from the core gut-associated taxa was revealed. Several genera (Pseudomonas, Morganella, Klebsiella, Serratia, and Providencia), belonging to families Pseudomonadaceae and Enterobacteriaceae, were present with high RAs mainly in bees from the Danube plain area and occasionally in a few individuals from the Rhodope Mtn. We relate the presence of these bacteria to body surface-associated microbiota. These genera have widespread distribution, found in various environments (water, soil, plants, insects, and other animals including humans). However, some species are also known as potential and opportunistic bee pathogens that can disrupt gut microbiota and increase host susceptibility to infection and reduced lifespan (Kwong and Moran, 2016; Raymann et al., 2017; Lakhman et al., 2021; Al Naggar et al., 2022). When foraging on flowers, honey bees can transfer numerous pollen-associated bacteria that differ significantly in plant species and pollination type (Manirajan et al., 2018), including some plant pathogenic species (e.g., Pseudomonas), thereby facilitating species-specific crop diseases (McArt et al., 2014; Pattemore et al., 2014).

Alpha diversity indices and beta diversity analyses confirmed the observed differences in honey bee-associated microbiota from both localities. One of the possible explanations for the observed community shifts in bacterial composition and higher diversity in bees from the Danube plain locality might be related to the higher infectious rate with Nosema microsporidian parasites found in the same apiary (Radoslavov et al., 2017). Alterations in the community composition and diversity of bee gut microbiome were found to be positively correlated with the presence of eukaryotic parasites and pathogens (Cariveau et al., 2014; Hubert et al., 2017; Felden et al., 2021). Moreover, a higher microbial diversity was associated with a lower bumblebee resistance against infection by an intestinal parasite (Näpflin and Schmid-Hempel, 2018).

The presence and higher diversity of non-specific gut bacteria in bees from the lowland area may have been related to the landscape features of the location and the level of anthropogenic impact. This is a rural area with higher land use intensity. Increased exposure to agricultural chemicals due to the prevalence of agroecosystems may have had long-term consequences on bee health and its microbiome structure and function. As a main pollinator of various crops, often located in areas with high land use intensity, A. mellifera has diverse interactions with the local environment that can significantly affect bee health and colony performance. In general, each location is a complex of interacting local, global, and regional environmental factors, some of which acting as bee stressors. Land use intensity, agrochemicals, and nutrition quality have been found to contribute to alterations in the honey bee-associated microbiome (Kakumanu et al., 2016; Jones et al., 2018; Motta et al., 2018; Muñoz Colmenero et al., 2020; Cullen et al., 2023). In contrast to the lowland area, the topography of the Rhodope Mtn creates a heterogeneous environment with more diverse landscape structures and ecosystems owing to the variations in temperature, vegetation, lack of homogeneous monocultures, and lower pollution levels. This leads to a larger diversity of non-managed habitats, such as natural grasslands, low-intensity pastures, forests, patches of shrubs, and single trees that provide higher plant diversity and nutritional resources for honey bees. Both sites differed in many aspects, including elevation, landscape heterogeneity, land use type and vegetation, access to floral resources, application of agrochemicals that may affect alone or in combination. However, the lower number of sampling locations, smaller sampling sizes, and the combined effect of multiple abiotic and biotic factors hinder the clear interpretation of our findings. Therefore, further studies with higher sampling rate, involving pathogen and fungal community analyses, would be required to disentangle the influence of local and regional factors.