Introduction

Bees are vital to global ecosystems, functioning as primary pollinators and enhancing biodiversity, which in turn supports agricultural productivity and food security worldwide1,2. The honey bee (Apis mellifera), as a key pollinator, directly boosts crop yields and underpins biodiversity across diverse landscapes, from managed farmland to wild habitats.

The gut microbiota of bees plays a critical role in their health and productivity, fulfilling essential nutritional, immune, and physiological functions. This microbiota is largely maintained through both horizontal and vertical microbial transmission, which contributes to its resilience and adaptability across bee generations3,4,5,6,7. Additionally, the vertical inheritance of microbial communities has shaped the gut microbiota composition of corbiculate bees over evolutionary timescales6, resulting in a relatively conserved core microbiota that comprises fewer than ten major bacterial taxa. These include Snodgrassella alvi, Gilliamella apicola, Lactobacillus spp., Bifidobacterium spp., Frischella perrara, Bartonella apis, Parasaccharibacter apium, and Commensalibacter spp8. These microorganisms support key digestive and immune functions, promote weight gain, and enhance bees’ resistance to pathogens, highlighting the evolutionary importance of this symbiosis9,10,11,12.

Studies have demonstrated the influence of different environmental factors, such as urban and rural landscapes, on bee gut microbiota composition5,13,14,15,16,17,18,19,20. However, there is a lack of information on how the gut microbial structure of bees is shaped across different biomes in Brazil. Understanding how landscapes affect the gut microbiota of bees can help improve conservation efforts and sustainable management of these ecosystems. Despite significant progress in understanding bee gut microbiota globally, there remains a knowledge gap concerning the microbiota diversity within the specific biomes of Brazil. This study addresses this gap by conducting a comparative analysis of the gut microbiota composition of Apis mellifera from two ecologically contrasting Brazilian biomes: the Atlantic Forest and the Caatinga.

Materials and methods

Biome characterization

The study was performed in two distinct biomes in Paraíba State, Northeastern Brazil: The Atlantic Forest and the Caatinga. These biomes are remarkably different in terms of soil composition, climate characteristics, and vegetation21,22. The Atlantic Forest isa type of tropical montane forest thatextends from the state of Rio Grande do Norte to Rio Grande do Sul. It encompasses a diversity of forest ecosystems with markedly differentiated floristic structures and compositions shaped by the specific climatic characteristics of each region23. In contrast, the Caatinga, located in Northeastern Brazil, is characterized by a semi-arid climate, vegetation adapted to aridity, and a remarkable diversity of species24. Within the state of Paraíba State, the Atlantic Forest covers approximately 7.31% of the land area, while the Caatinga covers 92.69%25.

Study design and samplings

Five pooled samples, each consisting of an average of 35 Africanized Apis mellifera nurse bees, were collected from five different hives in apiaries located within the two contrasting biomes. The five samples from the Atlantic Forest region were collected from the Beekeeping Sector of the Federal University of Paraíba (Latitude: 6°58’19.77"S, Longitude: 35°43’16.88"W, Altitude 522 m). The other five samples were collected from nested hives in the apiary of São João do Cariri, belonging to the Experimental Station of UFPB (7° 12’ 43’’ S, 36° 37’ 46’’ W; altitude 510 m above sea level) in the Caatinga region (Fig. 1).

Fig. 1
figure 1

Geographical locations of apiaries sampled in the Atlantic Forest and Caatinga biomes of Brazil. Satellite imagery obtained from Google Earth Pro Google Earth Pro (version 7.3.6.9796, https://earth.google.com) and subsequently edited in Microsoft PowerPoint (Office 2013). S - São João do Cariri – PB, Latitude: 7°23’0.77"S, Longitude: 36°31’43.07"W, Altitude 428 m, Precipitation: 381.4 mm; A - Areia - PB, Latitude: 6°58’19.77"S, Longitude: 35°43’16.88"W, Altitude 522 m, Precipitation: 1,400 mm.

The samplings were carried out in the dry period of the year, by the end of November, which is characterized by increased temperatures and the absence of rainfall in both the Atlantic Forest and the Caatinga regions. Appropriate personal protective equipment necessary for the management of Apis mellifera species bees was used during sample collection. The bees were placed into sterile tubes containing 70% ethanol and stored at -20ºC until dissection.

Sample preparation and DNA extraction

The bee intestines were dissected in a sterile environment. Before extraction, the bee specimens were left on a sterile filter paper for 10 min to thaw and allow the ethanol to evaporate. In the dissection process, a cross-section was made in the last abdominal segment of each bee, and the gut contents were then transferred to sterile microtubes. Each sample represented pooled gut specimens from 35 individual bees.

DNA was extracted using a commercial kit (PowerSoil, Qiagen, Hilden, Germany), according to the manufacturer’s instructions. DNA integrity was evaluated by agarose gel electrophoresis, and concentration was determined using a microvolume spectrophotometer (Colibri LB 915, Titertek-Berthold, Germany).

Library preparation and sequencing

The libraries for 16 S rRNA metabarcoding were prepared as previously described following a reference protocol (Illumina, 2013).

The V3-V4 hypervariable region of the 16 S rRNA gene was amplified by PCR. The 25-µL master mix consisted of 2.5 µL of DNA template (5 ng/µL), 5 µL of each primer (341 F: 5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG CCT ACG GGN GGC WGC A-3′ and 805R: 5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GGA CTA CHV GGG TAT CTA ATC C-3′), and 12 µL 2× KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, MA, USA). The thermal cycling conditions consisted of an initial denaturation at 95 °C for 3 min, followed by 25 cycles at 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min.

Amplification products were visualized on 1.5% agarose gel before purification using magnetic beads (AMPure XP, Beckman Coulter, USA). The dual indices and Illumina sequencing adapters were attached using the Nextera XT Index Kit (Illumina). A second clean-up step was performed using magnetic beads. Purified PCR products were quantified by fluorometry (Qubit 2.0, Life Invitrogen, USA). Pooled libraries were denatured with NaOH, diluted with hybridization buffer, and heat-denatured. Paired-end sequencing was performed on an Illumina MiSeq platform with a V2 kit (2 × 250 cycles), with 15% PhiX DNA (Illumina) as a sequencing control.

Bioinformatics and statistical analysis

The raw demultiplexed paired-end sequences were processed using QIIME 2-2020.226. Reads were filtered, denoised, truncated to a length of 248 base pairs, and parsed for non-chimeric sequences using DADA2, generating Amplicon Sequence Variants (ASV)27. Sequences were aligned using the “qiime fragment-insertion sepp” command for phylogenetic analysis28. The taxonomic composition of the samples was determined using a pretrained naive Bayes classifier with a 99% sequence similarity threshold for V3-V4 reference sequences (SILVA-132-99-nb-classifier.qza) and the “Qiime feature-classifier classify-sklearn” command. Negative control samples were examined for potential contaminants, and no taxa overlapped with the true samples.

Microbial diversity was quantified using Pielou’s index for evenness and the Shannon index for richness and abundance. ANOVAs were used to compare diversity between groups in R version 4.1.029 after testing for normality using the Shapiro-Wilk test. For alpha diversity assessment, we used phylogenetic diversity indices: Shannon (abundance and richness), Faith (phylogenetic diversity), Pielou (uniformity of abundance distribution), and number of observed ASVs in QIIME 2-2020.226. Differences in alpha diversity between samples from the two biomes were assessed in R version 4.1.029 after normality testing using Shapiro-Wilk test. When data were normally distributed, we used ANOVA followed by Tukey’s post hoc test, and when data did not meet normality assumptions, we applied the Kruskal-Wallis test followed by Wilcoxon rank-sum tests at a 5% significance level (α = 0.05) to compare microbial diversity between honey bee biomes, using the laercio package version 1.0–030.

For qualitative beta diversity analysis, we used the Jaccard coefficient to evaluate the dissimilarity across samples regarding gut microbial structure based on presence and absence of taxa, and Unweighted UniFrac as a phylogenetic-based metric for dissimilarity assessment. For quantitative beta diversity, we used the Bray-Curtis coefficient and Weighted UniFrac metrics. The Bray-Curtis coefficient was chosen because it combines both abundance and presence/absence data for the compositional dissimilarity assessment, while Weighted UniFrac is a phylogenetic-oriented metric for taxa. Beta diversity analyses were performed using the QIIME 2-2020.2 platform26. PERMANOVA was used for each metric to test putative differences in taxa composition between biomes31.

The core taxa (defined as taxa present in all samples) were determined using the “qiime feature-table core-features” command. Analysis of Composition of Microbiomes (ANCOM) was applied to assess differences in relative abundances between biomes. The relative abundances of ASVs observed in core taxa and identified as significant by ANCOM were assessed. Data were tested for normality using the Shapiro-Wilk test, and an analysis of variance (ANOVA) was applied to identify significant differences at P < 0.05. The Emperor 2020.2.0 plugin in QIIME 2-2020.226 was used to visualize beta diversity outputs. We used R version 4.1.029 to generate alpha diversity figures.

Results

16 S rRNA metabarcoding outcomes

A total of 623,229 reads were obtained across all 10 sampled pools, with 590,517 reads (94.75%) passing initial quality filtering and 585,109 reads (93.88%) retained after denoising. Sequence processing yielded 614 unique Amplicon Sequence Variants (ASVs), with an overall read count of 341,529. Taxonomic assignment was successful for 98% of the sequences at the genus level.

Alpha diversity analysis of gut microbiota)

Alpha diversity indices, including the Shannon index, Faith’s phylogenetic diversity (PD), Pielou’s evenness, and observed features, showed no statistically significant differences in microbial diversity between honey bees from the Atlantic Forest and Caatinga biomes (ANOVA: Shannon, p = 0.614; Faith’s PD, p = 0.311; Pielou’s Evenness, p = 0.471; Observed Features, p = 1.000) (Fig. 2A–D).

Fig. 2
figure 2

Box plots showing alpha diversity indices (the Shannon, Faith’s phylogenetic diversity (PD), Pielou’s evenness, and observed ASVs) for gut microbiota in honey bees from the Caatinga and Atlantic Forest biomes. Box plots illustrate outliers, first and third quartiles, and median values for samples from the Caatinga (Caa) and Atlantic Forest (Mata) groups, with no significant differences observed (Shannon index, p = 0.614; observed ASV, p = 1; Faith’s PD, p = 0.311; Pielou’s evenness, p = 0.471).

Beta diversity analysis of microbial composition

Beta diversity analyses revealed significant differentiation in microbial composition between the Atlantic Forest and Caatinga honey bee populations, as assessed by the Bray-Curtis, Jaccard, and Unweighted UniFrac indices (PERMANOVA: Bray-Curtis, p = 0.015; Jaccard, p = 0.007; Unweighted UniFrac, p = 0.010). The Weighted UniFrac analysis did not reveal significant differences between the groups (p = 0.271) (Fig. 3A–D). The Bray-Curtis index indicated compositional differences, suggesting that microbial community structure is influenced by biome-specific environmental factors. Similarly, the Jaccard analysis underscored significant qualitative divergence in gut microbiota composition between biomes, emphasizing the effect of presence-absence patterns in microbial taxa.

Fig. 3
figure 3

Beta diversity of honey bee gut microbiota in the Caatinga and Atlantic Forest. Box plots display differences in microbial community composition (PERMANOVA: Bray-Curtis, p = 0.015; Jaccard, p = 0.007; Unweighted UniFrac, p = 0.010; Weighted UniFrac, p = 0.271).

Core Microbiota and differential abundance analysis (ANCOM)

Using ANCOM, we identified Apibacter as a differentially abundant genus between biomes (ANOVA: Apibacter, p = 0.024, W = 7) (Fig. 4A). Core microbiota analysis identified seven genera consistently present in all samples: Lactobacillus, Commensalibacter, Rhizobiaceae, Snodgrassella, Gilliamella, Orbaceae, and Bifidobacterium, collectively accounting for 63% of the total genera. No significant differences in relative abundance were observed among these core genera (ANOVA: Lactobacillus, p = 0.822; Commensalibacter, p = 0.299; Rhizobiaceae, p = 0.161; Snodgrassella, p = 0.415; Gilliamella, p = 0.717; Orbaceae, p = 0.619; Bifidobacterium, p = 0.248) (Fig. 4B–H; Supplementary Figs. 1 and 2).

Fig. 4
figure 4

Relative abundances of core bacterial genera in the gut microbiota of Africanized honey bees from the Atlantic Forest and Caatinga biomes. Box plots show the median, first and third quartiles, and outliers for each biome. Apibacter was the only genus showing significant differences between biomes (p = 0.024). No significant differences in abundance were detected for other core genera.

Discussion

The identification of 614 ASVs (614) in the gut microbiota of honey bees from both biomes indicates considerable microbial richness, consistent with previous studies reporting similar diversity levels in the gut microbiota of honeybees6,11,32. Although microbial composition differed significantly between biomes, alpha diversity metrics showed that richness and evenness of the microbial communities within each biome were consistent. This finding suggests that environmental differences between apiary locations may not be the primary factor shaping alpha diversity in these bees´ gut microbiota. This observation aligns with a previous study6, showing that the microbiome of corbiculate bees tends to remain stable across different settings. Higher-resolution methods, such as strain-level identification through shotgun metagenomics, may reveal finer-scale landscape effects on microbial composition not detected in this study33.

Beta diversity analysis revealed significant differences in microbial composition between bees from the Caatinga and Atlantic Forest biomes. The gut microbiota composition differed between biomes, except when analyzed using the weighted UniFrac index. This index measures the evolutionary distances between microbial communities while accounting for bacterial taxa abundances in each group.

Contact among eusocial organisms, particularly in bee hive communities, maintains gut microbiota through microorganism transfer within populations4,34, which could explain the fact that microbiome composition correlates more strongly with host genetic and seasonality than with geography6,35.

Honey bees from both biomes exhibited significant differences in their gut microbial communities, with dissimilarities observed between Caatinga and the Atlantic Forest samples, as well as hives within each biome. The Jaccard similarity index, which measures qualitative dissimilarity, showed significant differences in bee microbiota beta diversity between biomes, with greater dissimilarity among Caatinga hives compared to Atlantic Forest hives. These patterns suggest that biome-specific conditions influence gut microbiota structure. Environmental factors, including floral resource availability, and climatic conditions play key roles in selecting bee gut microbiota. Microbial symbioses as crucial for insect socialization36. Consequently, both ecological parameters and behavioral traits drive persistent symbiosis across the solitary-social insect spectrum spectrum37. The hive structure is further influenced by environmental factors, niche construction, nest substrate, and food storage practices38,39,40.

Our findings on the core microbiome (Lactobacillus, Commensalibacter, Snodgrassella, Gilliamella, Bifidobacterium, Rhizobiaceae and Orbaceae) in both biomes align with previous reports for Apis mellifera3,11,41,42,43,44. Species of the genus Apibacter encode a type VI secretion system (T6SS)45, which confers resistance to pathogenic bacteria through antibacterial protein synthesis. In our study, samples were collected during a period of food scarcity, particularly in the Caatinga region. This resource scarcity may have contributed to the decreased relative abundance of Apibacter. Resource limitations can significantly affect bee microbial populations through various mechanisms, impacting gut microbiota and overall colony health34,46,47,48. Further research should examine the relationship between Apibacter and plant resource availability, as declining populations of these bacteria can compromise hive health and increase suscepribility to pathogens,

The dominant microbial community composition varied with environmental features across landscapes. Similar patterns emerged across different metrics (Bray-Curtis, Jaccard, weighted and unweighted UniFrac), suggesting that environmental conditions, including site specific differences, influence gut microbial community composition and relative abundance of specific taxa49. Gut bacterial taxa abundance has been shown to vary with insect habitat, correlating with environmental oxygen availability50.

Honey bees from different biomes shared common bacterial families and genera, particularly Lactobacillus, Gilliamella, and Snodgrassella. The prevalence of these genera reflects the characteristic gut microbiota of Africanized bees, which comprises nine bacterial phylotypes predominantly found in the large intestine51.

Five species – Snodgrassella alvi, Gilliamella apicola, two Lactobacillus species, and one Bifidobacterium species – occur ubiquitously in corbiculate bees worldwide6; constituting the core gut microbiome. Other taxa, such as Bartonella apis, Apibacter adventoris, Frischella perrara and Acetobacteraceae, are variably present in adult bee guts. While environmental species primarily occur in the foregut52. , core microbiota species show distinct spatial distribution within the large intestine. S. alvi and G. apicola dominate the ileum region, with S. alvi forming a continuous layer along the longitudinal folds, and G. apicola occurring above it53.

The Atlantic Forest bees showed higher abundance of the genus Apibacter, Commensalibacter and Rhizobiaceae. Gilliamella, Lactobacillus, Snodgrassella, and Enterobacteriaceae exhibited low abundance in both biomes, while Frischella, Bifidobacterium, and Proteobacterium were relatively abundant.

We acknowledge several limitations in this study that should be considered for a correct and transparent interpretation of our findings. Firstly, the current study employed five pooled samples per biome only, which presents inherent constraints on statistical power and generalizability. This sample size limitation potentially restricts our ability to comprehensively capture the full microbiota variability within and between the Atlantic Forest and Caatinga biomes. We recognize that a larger, more diverse sample set would provide more robust and representative insights into the gut microbiome dynamics of Apis mellifera across these distinct ecological regions. Moreover, while our research identified significant differences in gut microbiota composition between the two biomes, we explicitly acknowledge the complexity of environmental interactions. The current study did not establish direct causal relationships between specific environmental variables (including floral diversity, pesticide exposure, and climatic conditions) and the observed microbiota variations. These environmental factors likely operate through intricate, multilayered interactions that profoundly influence microbial community structures.

In conclusion, this study provides important insights into Apis mellifera gut microbiota composition and variability across biomes. While environmental conditions influenced microbial composition, microbial diversity remained consistent across populations. These findings highlight the need for further research on how microbiota variations affect bee health and performance in different landscapes. The identification of a common core microbiota emphasizes the importance of specific microbial groups in maintaining bee health. This work stablishes a foundation for future studies examining the ecological roles of these microbes and their contribution to honey bee adaptation across environments. Future research should focus on expanding sample sizes, implementing comprehensive environmental data collection, and employing advanced system biology tools. Integrative approaches such as multi-omics analyses, network-based modeling, and machine learning algorithms could provide more nuanced insights into the complex interactions between environmental factors and Apis mellifera gut microbiota composition. By leveraging these computational and analytical strategies, we could develop a more holistic understanding of microbiome variations across different biomes, ultimately enhancing the ecological and biological comprehension of honey bee microbiome dynamics.