Global distribution of honeybee gut microbiome and pesticide-driven adaptations in opportunistic microbial species

Global distribution of honeybee gut microbiome and pesticide-driven adaptations in opportunistic microbial species | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Global distribution of honeybee gut microbiome and pesticide-driven adaptations in opportunistic microbial species View ORCID Profile Nelli Vardazaryan , View ORCID Profile Lusine Adunts , View ORCID Profile Inga Bazukyan , View ORCID Profile Magdalina Zakharyan , View ORCID Profile Honglian Liu , View ORCID Profile Chrats Melkonian , View ORCID Profile Lilit Nersisyan doi: https://doi.org/10.1101/2025.02.25.640077 Nelli Vardazaryan 1 Armenian Bioinformatics Institute , Yerevan, Armenia 2 Institute of Molecular Biology, National Academy of Sciences of Armenia , Yerevan, Armenia Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Nelli Vardazaryan For correspondence: lilit.nersisyan{at}abi.am nelli.vardazaryan{at}abi.am Lusine Adunts 1 Armenian Bioinformatics Institute , Yerevan, Armenia 2 Institute of Molecular Biology, National Academy of Sciences of Armenia , Yerevan, Armenia Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Lusine Adunts Inga Bazukyan 3 Yerevan State University, Institute of Biology , Yerevan, Armenia Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Inga Bazukyan Magdalina Zakharyan 2 Institute of Molecular Biology, National Academy of Sciences of Armenia , Yerevan, Armenia Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Magdalina Zakharyan Honglian Liu 4 SciLifeLab, Department of Microbiology, Tumor and Cell Biology, Karolinska Institutet , Solna, Sweden Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Honglian Liu Chrats Melkonian 5 Theoretical Biology and Bioinformatics, Department of Biology, Faculty of Science, Utrecht University , 3584 CH Utrecht, The Netherlands 6 Bioinformatics Group, Wageningen University and Research , Wageningen, The Netherlands Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Chrats Melkonian Lilit Nersisyan 1 Armenian Bioinformatics Institute , Yerevan, Armenia 2 Institute of Molecular Biology, National Academy of Sciences of Armenia , Yerevan, Armenia Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Lilit Nersisyan For correspondence: lilit.nersisyan{at}abi.am nelli.vardazaryan{at}abi.am Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Honeybees ( Apis mellifera ) rely on a specialized gut microbiome shaped by climate, flora, agrochemicals, and dietary supplements. Yet, how these factors alter microbiome composition and function remains unclear. We integrated 16S rRNA and shotgun metagenomics data from eight published studies across six regions, alongside newly generated 16S data from Armenia, to assess how environmental and agrochemical factors influence the honeybee gut microbiome. We also introduced a novel co-abundance and functional analysis pipeline to identify treatment-affected bacterial networks and associated pathways. We observed a stable core of six to twelve phylotypes in amplicon and metagenomic datasets respectively. However, we note significant geographic variation in relative abundances, likely reflecting differences in diet, climate, and local flora. Armenian data revealed distinct seasonal shifts, particularly elevated Commensalibacter in autumn and minor urban-rural differences. Pesticide treatments elicited varying responses: oxalic acid drove pronounced beta-diversity shifts; neonicotinoids had subtler effects, both primarily impacting opportunistic pathogens; and glyphosate disrupted core taxa with stronger effects in newly emerged bees under prolonged exposure. Co-abundance network analysis highlighted that the pesticide-associated community was enriched in adaptive pathways, including potential glyphosate degradation by Pseudomonas , biofilm formation, and aromatic amino acid synthesis. These findings reaffirm the stability of the honeybee core microbiome yet underscore that environmental and anthropogenic stressors induce distinct compositional and functional shifts. We emphasize the need for longitudinal metagenomic approaches that enable high-resolution functional profiling and co-abundance network analysis to clarify how these microbiome shifts impact bee health and colony sustainability. Introduction The Western honey bee ( Apis mellifera ) is a key pollinator, renowned for its broad geographic distribution and effective foraging [ 1 , 2 ]. In recent years, honeybee populations have declined globally, with annual losses reaching 30-45% [ 3 – 6 ]. Several factors have been linked to this decline, including climate change [ 6 ], heavy use of pesticides [ 7 ], use of antibiotics [ 8 ], and pathogen infection [ 9 ], but no single factor has yet been implicated as the cause. The gut microbiome plays a vital role in bee health by regulating metabolism, pollen digestion, immunity, development, and behavior [ 10 ]. The core microbiome of adult workers comprises highly specialized bacteria. Five consistently reported phylotypes include Proteobacteria - Gilliamella ( G. apicola ) and Snodgrassella ( S. alvi ), which mostly reside in the ileum and form a biofilm; and Bacillota - Bombilactobacillus ( a.k.a Lactobacillus Firm-4 ), Lactobacillus ( a.k.a. Lactobacillus Firm-5), and Bifidobacterium , predominantly found in the rectum [ 11 – 14 ]. Other genera, such as Frischella, Bartonella, Commensalibacter, and Bombella , are also common [ 11 , 15 , 16 ], while Apilactobacillus kunkeei, Apibacter, Serratia marcescens, and some Enterobacteriaceae, are sometimes reported at low abundance [ 17 ]. Environmental factors, health status, and diet strongly affect the honeybee gut microbiome, often leading to increased mortality [ 4 , 18 ]. In particular, agricultural pesticides can disrupt the microbiota, altering behavior and immune function [ 19 – 21 ]. Neonicotinoid insecticides, such as thiacloprid, clothianidin, and imidacloprid, not only directly impair the nervous system affecting the memory and locomotion of honeybees, but may also disrupt the gut microbiome [ 22 – 25 ]. Oxalic acid, an acaricide insecticide commonly used against Varroa parasites and considered relatively safe for bees [ 26 ], can nonetheless lower the gut pH and alter the microbiome composition at high doses [ 19 ]. Glyphosate, the most widely used herbicide globally, negatively affects honeybee behavior and physiology [ 27 ]. In plants, it targets the EPSPS enzyme in the shikimate pathway, crucial for the synthesis of essential aromatic amino acids [ 14 ]. Since many bacteria possess glyphosate-sensitive class I EPSPS enzymes, concerns present about its adverse impact on the bee gut microbiome [ 14 , 28 – 30 ]. Among the glyphosate-affected species, S. alvi is most consistently reported. However, the overall impact of glyphosate on the gut microbiome remains inconclusive [ 4 , 29 , 30 ]. Although the honeybee gut microbiome is widely recognized as essential for bee health, conflicting evidence remains regarding how pesticides affect its taxonomic composition and functional potential. For example, studies on neonicotinoids yield contradictory results, some report negligible impacts under short-term or laboratory conditions, while others document significant shifts after long-term field exposure [ 25 , 31 ]. Similarly, while in vitro experiments show that oxalic acid increases the abundance of Gilliamella, Bifidobacteria, Lactobacillus , and Snodgrassella , these effects are not consistently observed in vivo [ 19 ]. The impact of glyphosate also varies considerably with dose, exposure duration, and the state of gut colonization [ 29 , 32 , 33 ]. Importantly, these studies do not fully explain the functional consequences or reveal the adaptation mechanisms underlying these taxonomic shifts. To address these gaps, we performed a meta-analysis by aggregating honeybee gut microbiome datasets from eight published studies across six geographic regions and incorporated newly generated 16S data from Armenia, an understudied region in apicultural research. We employed a unified bioinformatics workflow to reprocess the data, thereby reducing technical variability and enabling robust comparisons of community structure and predicted functional capacity across diverse environments and pesticide treatments. Furthermore, we introduced a novel co-abundance network analysis to identify treatment-affected bacterial communities and elucidate their putative functional roles. Materials and methods Dataset search and inclusion criteria We conducted a systematic search for datasets in the NCBI SRA database using the terms: (microbiome AND “gut” AND (“Apis mellifera” OR “honeybee”)) . We included both 16S rRNA amplicon and whole-genome sequencing (WGS) datasets. Relevant datasets were additionally gathered from recent literature. The inclusion criteria for studies were as follows: (i) paired-end sequencing; (ii) study subjects being worker bees; (iii) inclusion of the entire gut (including crop, midgut and hindgut); (iv) availability of raw data and corresponding metadata; (v) for 16s amplicon sequencing, libraries targeting the V4 region (e.g., V4 only, as well as V3-V4) to minimize batch effects related to variable amplicon region choice and enable cross-study comparison [ 34 , 35 ]. The data was downloaded from SRA using the fastq-dump tool from the NCBI SRA Toolkit, or downloaded directly through personal communication with authors. Detailed information about the datasets used is available in Table S1. DNA extraction and sequencing of the Armenian samples Honeybees from different colonies were collected from the Gegharkunik, Kotayk, Syunik, Vayots Dzor, Lori, Ararat, Artashat, Jermuk, Hanqavan, and Aghavnadzor provinces of Armenia. Samples were retrieved from a single colony in each province to limit the effect of genetic variation. Guts were dissected using flame-sterilized forceps under aseptic conditions. The head was separated from the thorax to detach the esophagus, and the entire intestinal tract was carefully extracted through the anterior end by grasping the last abdominal segment with forceps and gently pulling. Additional segments of the abdominal exoskeleton were carefully detached to ensure the rectum was removed intact, as it was notably enlarged in most bees. DNA was extracted using the cetyltrimethylammonium bromide (CTAB) bead-beating method, followed by phenol-chloroform purification and ethanol precipitation, as described in [ 36 ]. The extracted DNA was sent to Macrogen Inc. (Seoul, Korea) for paired-end 16S sequencing on the Illumina MiSeq platform. The V3-V4 hypervariable regions of the 16S rRNA gene were amplified using the primers Bakt_341F (5′-CCTACGGGNGGCWGCAG-3′) and Bakt_805R (5′-GACTACHVGGGTATCTAATCC-3′). 16S sequencing data processing Raw paired-end reads were filtered out if they contained any of the following: (i) ambiguous bases, (ii) homopolymers longer than eight bp (using fastq_qual_trimmer v1.0), or (iii) quality scores less than 28 (using fastq-mcf v.1.0). The trimmed sequence reads were processed using the dada2 package for R (v.1.22.0) to infer exact amplicon sequence variants (ASVs) [ 37 ]. Briefly, reads were truncated with the filterAndTrim function, specifying the truncLen parameter to trim bases with a quality score below 30; reads containing more than two errors (forward) and three errors (reverse) were removed with the maxEE parameter; sequence reads were dereplicated, denoised, and merged using dada2 default parameters, with pooled sample inference applied for each dataset. The datasets were merged into a single ASV table using the mergeSequenceTables function of the phyloseq package (v.1.48) [ 38 ]. To reduce batch effects, ASVs were truncated to the common V4 region of the 16S rRNA gene through the following steps: (i) sequences were aligned using the msa R package v.1.30.1 [ 39 ]; (ii) aligned ASVs were truncated from the 3′ end, based on the C3 conserved sequence region (GTGCCAGCAGCCGCGGTAAA), allowing up to five possible mismatches; (iii) reads were further truncated at the 5′ end, to the most abundant length of 60 bp; (iv) alignment gaps were removed. Taxonomy was assigned to the ASVs using the SILVA database [ 40 ] via the assignTaxonomy function in dada2 . The taxonomy-assigned data was converted into a phyloseq object. Chimeric ASVs and those observed in only one sample were filtered out. ASVs unclassified at the genus level, assigned to mitochondria, chloroplasts, or eukaryotes, or not assigned to bacterial or viral lineages were excluded. To improve taxonomic accuracy, we retained only taxa meeting the following criteria: (i) at least 10 ASVs in one or more samples, and (ii) present in more than 70% of samples in each dataset or treatment group within a study. Rarefaction curves were constructed using the rarefy function from the vegan package (v.2.6.8) to assess the number of observed ASVs relative to library size. Each sample was rarified to 12,000 reads using the rarefy_even_depth function in phyloseq ( Fig. S1 ). WGS data processing Taxonomic classification of the WSG reads was performed with Kraken 2 (v.2.1.2) [ 41 ] with default settings and the NCBI Standard database ( https://benlangmead.github.io/aws-indexes/k2 , accessed on 17 May 2022), without removal of host sequences. 1-9% of reads were assigned to prokaryotes. Species-level counts were re-estimated using Bracken (v.2.9). Downstream analysis of Bracken results was conducted using the phyloseq package [ 42 ]. The species abundance table was rarified to 800 000 reads (Fig. S1). Diversity analysis and statistics Shannon alpha diversity indices were calculated using phyloseq. For community composition (multivariate) analyses, pairwise Bray-Curtis dissimilarity between samples was estimated, and principal coordinates were computed using the dist_calc and ord_calc functions in the microViz package v.0.12.4 [ 43 ]. PERMANOVA analysis was performed using the adonis2 function in the vegan package v.2.6.8. The envfit function in vegan was used to fit environmental variables onto the ordination axes using multiple regression. The statistical significance of these relationships was determined using a permutation test. The fitted environmental variables were visualized as vectors, projected onto the ordination diagram using the geom_segment function in ggplot2 (v.2 3.5.1). Differential abundance analysis To identify species or genera associated with treatment, we performed differential abundance analysis with the ANCOM-BC2 v.2.6.0 [ 44 ]. Unlike other methods, ANCOM-BC2 assumes that, under the null hypothesis, the proportions of taxa in a sample remain constant across conditions, making it particularly effective in addressing the biases and covariances common in microbiome datasets. It also adjusts for biases derived from varying sequencing and extraction methods. While running ANCOM-BC2, we encountered an error related to object length mismatches that we corrected in our pipeline using the steps suggested by the community ( https://github.com/FrederickHuangLin/ANCOM-BC/issues/231 ). Co-abundance network analysis We developed a new methodology and an R package, micronet , for the construction and functional analysis of microbial taxon-taxon co-abundance networks. We first selected taxa with at least 0.05% and a prevalence of 20% or more in at least one treatment group (e.g. control, treatment, and, if applicable, experiment group). If some taxa within a genus were identified at the species level, the rest of the unidentified ASVs were aggregated and marked with the genus name and “other” prefix. For genera with no ASV identified at the species level, we aggregated those and marked them with the genus name (Table S10). We calculated Spearman correlation coefficients between the relative abundances of taxa within each treatment group in each study. To ensure robustness, we performed bootstrapping using the co_occurrence function of the phylosmith R package: in 200 iterations per treatment group, a random subset of samples was selected, and only taxa pairs with an absolute correlation coefficient >= 0.6 and FDR-adjusted p-value < 0.1 were retained. We kept taxon pairs that appeared in at least 80% of the iterations with a consistent correlation sign. These consistently correlated pairs were then merged across groups into a single co-abundance network. Nodes represented species or genera, and edges represented significant co-abundances, with edge weights defined as the median correlation coefficient across the groups where the correlation was observed. This pipeline is implemented in the generate_network function of the micronet package. Next, we aimed to detect treatment-affected subnetworks. We computed log-fold changes (logFC) of taxon relative abundances between treatment and control groups with ANCOM-BC2 . The absolute logFC values were discretized into node weights as follows: |logFC| > 1.5: weight 3; |logFC| between 1 and 1.5: weight 2; |logFC| between 0.5 and 1: weight 1; |logFC| < 0.5: weight -2; missing values: weight 0. We then extracted treatment-associated subnetworks using the maximum weight connected subgraph (MWCS) algorithm implemented in the mwcsr R package (v.0.1.9). MWCS identifies a connected group of nodes that maximizes the sum of their weights, representing the most relevant set of connected taxa. We have implemented this pipeline into the subgraph_analysis function of micronet . Finally, the networks were visualized in Cytoscape v.3.9.0 using an edge-weighted spring-embedded layout, in which nodes attract or repel each other based on their correlation coefficients [ 45 ]. Genome assembly and taxonomic characterization Metagenome-assembled genomes (MAGs) were reconstructed from WGS datasets. The resulting MAGs and functional annotations were later used as proxies to enhance functional calls in 16S datasets (see below). Raw library sizes ranged from 2.1×10 6 to 1.0×10 7 reads. Reads were trimmed and quality-controlled with fastp (v.0.23.2); those 10% low-quality bases (phred quality 5 ambiguous bases, or with low complexity (with repetitive or homogeneous patterns) were discarded, and Poly-G/X tails were trimmed. Host DNA was eliminated by aligning reads to the Apis mellifera genome (GCF_003254395.2), using Bowtie2 (v.1.3.1). Post-processing yielded 7.8×10 4 to 2.2×10 6 paired-end reads per sample. Both WGS datasets were co-assembled with MEGAHIT (v.1.2.9), and contigs were binned using MetaBAT, CONCOCT, and MaxBin. The three bin sets were merged via MetaWRAP’s bin_refinement module (v.1.3 [ 46 ]). Bin quality was assessed with CheckM (v.1.0.18); only bins with > 80% completeness and <5% contamination were retained, yielding 15 bins. Finally, MAGs were taxonomically annotated using GTDB-Tk (v.2.4.0) with default parameters [ 47 ] (Table S2). Gene function annotation We predicted microbial gene functions using PICRUSt2 (v.2.4.1) [ 48 , 49 ]. ASVs with species/genus level taxonomic assignments were used to annotate genes with KEGG Orthology (KO). We refined KO predictions using functional profiles from corresponding WGS-derived MAGs to overcome the limitations of high-level ASV taxonomic annotations. MAGs were annotated with eggNOG-mapper (v.2.1.12) [ 49 ]. For species with both MAG and ASV annotations, only KOs from MAGs were retained (Table S2). Pathway enrichment analysis We analyzed the refined KO table to identify enriched KOs in treatment-affected microbial communities. For each subnetwork, we separated nodes into up- or down-regulated groups based on positive or negative logFC values. KO-related gene counts in these groups were compared against the entire network using a hypergeometric test in R’s stats package (v.4.0.2). We computed the total probability of observing the given count (with dhyper ) or any higher count (with phyper , lower.tail = FALSE). KOs with p-values < 0.1 were deemed significantly enriched. Next, we performed KEGG pathway overrepresentation analysis with enrichKEGG in the ClusterProfiler package (v.4.12.6) [ 50 ], using the full network’s KO distribution as a reference. Pathways with FDR-adjusted p-values ≤ 0.2 were considered enriched. Only pathways enriched among up-regulated nodes and not enriched in down-regulated nodes were retained. Results Building a cross-study dataset with batch effect correction We collected, reprocessed, and reanalyzed raw data from eight published studies. We systematically searched the NCBI SRA database for 16S and WGS datasets from adult workers. These datasets comprise samples collected from hives or sterile lab conditions, with or without pesticide or sugar supplement treatments, and span diverse geographic regions. Additionally, we generated an in-house 16S rRNA dataset from the entire guts of worker honeybees collected across Armenia ( Fig. 1 ). Download figure Open in new tab Figure 1 Geographic distribution and sample composition of studies included in the meta-analysis. The map shows study locations and the number of samples per treatment (red) and control (green) group. Colors indicate the sampling location: (AM: Armenia, DN: Denmark, NZ: New Zealand, US:ST: Southern US, CH: China, UK: United Kingdom). Shapes indicate whether the samples are from 16S rRNA or WGS sequencing. In total, we analyzed 55 WGS and 461 16S samples. All WGS samples were untreated, whereas the 16S dataset comprised 190 untreated controls and 242 treated samples ( Table 1 ). Treatments included herbicides (glyphosate), insecticides (oxalic acid and neonicotinoids), and sugar supplements. Mean library sizes varied widely, from 10K to 160K reads per library for 16S and WGS datasets, respectively. 57.2% of 16S samples targeted the V4 and 42.8% - the V3–V4 regions ( Table 1 , Fig. S1,and Table S1). View this table: View inline View popup Download powerpoint Table 1. Summary of dataset features included in this study. To address potential batch effects in 16S datasets due to sequencing method heterogeneity [ 34 , 35 ], we performed Principal Coordinates Analysis (PCoA) based on Bray-Curtis dissimilarity. Clustering was primarily driven by variations in amplicon regions (PERMANOVA: R² = 0.135, p < 0.001) and library sizes (R² = 0.149, p < 0.001) (Fig. S3A). These effects were minimized after trimming libraries to a common V4 region and standardizing library sizes via rarefaction and filtering based on prevalence and abundance thresholds (see Fig. S3B and Methods). Although beta diversity differences related to study comparisons and experiment types (adult vs. newly emerged bees) persisted, we conclude that batch effects from technical variations, especially amplicon region and library size, are negligible after trimming and normalization and can be disregarded when comparing gut microbiome compositions across geographies. Honeybee gut microbiota is composed of a stable core with region-depend abundance patterns We first examined the most abundant taxa in untreated (field control) honeybee gut microbiomes using 16S datasets from Armenia, Denmark, New Zealand, and the southern US, along with WGS datasets from China and the UK. Consistent with prior research, the ten most abundant genera across all samples (both 16S and WGS) included: Bacillota : Lactobacillus ( formerly Lactobacillus Firm-5) and Bombilactobacillus ( formerly Firm-4) ; Pseudomonadota : Gilliamella, Snodgrassella, Frischella, Bartonella, Bombella, Commensalibacter, and Klebsiella; Actinobacteriota : Bifidobacterium [ 51 , 52 ]. Although Klebsiella, an opportunistic pathogen, was infrequently detected, it reached an average abundance of 3.2% in some samples ( Fig. 2A ). Download figure Open in new tab Figure 2. Comparison of honeybee gut microbiome composition across studies. (A) Relative abundances of the 10 most abundant bacterial genera in field control samples, with median relative abundances indicated for each study (AM: Armenia, DN: Denmark, NZ: New Zealand, US:ST: Southern US, CH: China, UK: United Kingdom). ( B) Prevalence and median relative abundance of core genera (present in >80% of field control samples) across studies. (C) Prevalence of species-resolved core taxa across all datasets. ( D) Alpha diversity of the honeybee gut microbiome in field control samples, measured using Shannon index; asterisks indicate significant pairwise differences (Kruskal-Wallis post-hoc test, p < 0.05). The central line indicates the median, the box represents the interquartile range (25th–75th percentiles), and whiskers extend to the 10th and 90th. (E) Principal coordinate analysis (PCoA) based on Bray-Curtis distances, with samples colored by study; between-group differences were assessed using PERMANOVA. Arrows denote the genera discriminating community composition (multiple regression test with permutations, p < 0.05). The microbiome composition was relatively stable within and between colonies (ANOSIM, R = 0.183, p = 0.024), yet varied significantly between states within studies (ANOSIM, R = 0.246, p = 0.002) and among studies (ANOSIM, R = 0.366, p = 0.001). In 16S datasets, Denmark and the southern US exhibited higher abundances of Gilliamella (26% and 24% median relative abundance) and Snodgrassella (19% and 16%), while Bombilactobacillus (21% and 22%) dominated in Armenia and New Zealand. Armenia also showed elevated Bifidobacterium levels. In contrast, WGS datasets were dominated by Bartonella (17%-31%), whereas 16S datasets were dominated by Lactobacillus (median 22%-42%). Bartonella was less prevalent in 16S datasets (0.4%–9%) and absent in the southern US ( Fig. 2B ). Next, we defined core taxa as those present in at least 80% of field control samples per study, analyzing 16S and WGS datasets separately ( Fig. 2B, C ). In the 16S datasets, six core taxa were identified: Lactobacillus Bombilactobacillus, Gillimella , Snodgrassella, Frischella, and Bifidobacterium . WGS datasets revealed six additional core genera: Entomomonas, Moraxella, Vibrio (Gammaproteobacteria) ; Bartonella, Apilactobacillus and Commensalibacter (Alphaproteobacteria). Notably, five core genera, Lactobacillus, Bombilactobacillus, Gillimella , Snodgrassella, and Bifidobacterium, are well documented [ 11 , 53 , 54 ]. In contrast, Bartonella, Commensalibacter, Frischella, and Apilactobacillus were reported as core only in specific studies [ 2 , 19 , 55 ]. In WGS datasets from the UK and China, we also found three low-abundance genera not previously reported as core: Etomomonas, a novel genus from Asian honeybee gut [ 56 ], Moraxella, previously identified in Apis mellifera [ 57 ] and linked to nicotinoid exposure [ 24 ], and Vibrio, enriched in pollen [ 58 ] (Table S3). At the species level, in the 16S datasets, Lactobacillus resolved to L. apis and L. melliventris , whereas WGS datasets revealed a broader range of core species, including L. apis, L. kullabergensis , L. helsingborgensis, and L. sp. wkB8 . Notably, L. melliventris was not identified by Kraken 2 in the China and UK WGS datasets despite its presence among reads, suggesting database annotation issues (see Supplementary Note 1, Fig. S4). Core Gilliamella species were resolved only in WGS samples, dominated by G. apicola along with G. sp. ESL0441, ESL0405, and ESL0443 . In Bifidobacterium, B. asteroides was core in both 16S and WGS, with B. coryneforme and B. indicum additionally core in WGS. In Snodgrassella, S. alvi was the only core species consistently identified in 16S and WGS. Additional WGS-only core species were Bombilactobacillus bombi, Frischella perrara, Bartonella apis , Commensalibacter sp. ESL0284 and AMU001, Apilactobacillus kunkeei, Entomomonas moraniae , Moraxella osloensis, and Vibrio anguillarum ( Fig. 2C ). We then compared microbiome diversity across studies. Alpha-diversity (Shannon index) was significantly higher in New Zealand compared to Southern US and Denmark (Kruskal-Wallis, post-hoc p < 0.001). Armenia also exhibited higher diversity than the Southern US ( p < 0.001) ( Fig. 2D , Table S4). PCoA using Bray-Curtis dissimilarity revealed that New Zealand and Armenia samples clustered distinctly from other studies, indicating greater within-group homogeneity ( Fig. 2E , Table S5). The features influencing discrimination of Armenia and New Zealand samples were increased Lactobacillus and Bombilactobacillus , and decreased Gilliamella ( Fig. 2E ). Notably, a subgroup of Southern US samples from Texas differed from those in Florida and Tennessee due to higher Pantoea abundance ( e.g. P. sp. SO10, P. vagans, P. agglomerans, P. ananatis ). These species may act as plant pathogens or contribute to pathogen defense in both plants and bees, with some linked to urban areas [ 59 – 63 ]. Overall, our analysis reveals a conserved core microbiota across diverse geographic regions, although their relative abundances vary significantly by location. We next investigated how seasonal and environmental factors shape the gut microbiota of honeybees in Armenia. Seasonal and environmental factors influence gut microbiota in Armenia, with increased Commensalibacter abundance in autumn We next analyzed original data from Armenia, where the gut microbiomes of worker honeybees were collected from nine regions during 2017, 2020, and 2023, in spring, summer, and autumn. In total, 21 samples, each comprising 5-6 pooled bee guts were examined. The microbiome was largely dominated by Lacotibacillus, Bombilactobacillus , and Bifidobacterium ( Fig. 3A ). Notably, some samples exhibited increased levels of potential pathogens such as Escherichia , and Shigella . Download figure Open in new tab Figure 3. Impacts of seasonality and urbanization on honeybee gut microbiome composition in Armenia. (A) Relative abundance of major taxa in Armenian honeybee gut samples, stratified by season, collection date, and anthropization level. (B) Alpha diversity of the gut microbiome (Shannon index) analyzed according to seasonal and anthropization factors. The central line indicates the median, the box represents the interquartile range (25th– 75th percentiles), and whiskers extend to the 10th and 90th. (C) Principal coordinate analysis (PCoA) based on Bray–Curtis distances; arrows indicate the most discriminant genera (multiple regression with permutations, p < 0.05). Samples are colored by collection season. The adjacent scatterplot highlights seasonal differences in Commensalibacter abundance, with group differences determined using ANCOM-BC2. (D) PCoA is colored by anthropization level, with an accompanying scatterplot illustrating differences in Bombilactobacillus abundance between urban and natural habitats. Seasonal comparisons among spring, summer, and autumn samples revealed significant differences between autumn and summer ( Fig 3C ; PERMANOVA, R 2 = 0.155, p = 0.03, Table S6). Multiple regression analysis indicated that these differences were primarily driven by higher abundances of Commensalibacter and Snodgrassella in autumn, and increased Bombilactobacillus in summer. Differential abundance analysis using ANCOM-BC2 confirmed a statistically significant rise in Commensalibacter during autumn (summer-autumn, q < 0.001, LFC = 3.94; spring-autumn, q < 0.016, LFC = 4.56) ( Fig. 3C , Table S7). Similar seasonal shifts have been observed elsewhere, with Swiss studies reporting increased Commensalibacter in winter bees [ 64 ], although other European [ 65 ] and Chinese studies have noted alternative patterns, e.g. increase in Bartonella [ 66 ] or a decrease in Snodgrasella [ 67 ]. We further assessed anthropogenic impacts by comparing bees from urban/near-urban apiaries (Ararat and Kotayk, near Yerevan) with those from rural or natural regions. PCoA revealed distinct microbiota profiles between urban and natural sites ( Fig. 3D ; PERMANOVA, R² = 0.119, p = 0.04), primarily due to an increased abundance of Bombilactobacillus in the urban samples. However, differential abundance analysis with ANCOM-BC2 did not yield statistically significant differences (q < 1, LFC = 0.32), likely due to the small sample size ( Fig. 3D , Table S7). To date, no specific taxon has been consistently linked to anthropogenic factors in urban regions. However, a decrease of Lactobacillus and Snodgrasella, with a simultaneous increase in Bifidobacterium and Enterobacteriacea, seen in near-urban regions, although not significant, was also reported previously [ 68 ]. Although some studies have reported increased microbial diversity in urban environments [ 61 , 69 , 70 ], we did not observe such an effect ( Fig. 3B ). Differential Impact of Pesticides and Diet on Honeybee Gut Microbiome We investigated the impact of pesticides, herbicides, and dietary supplements on the diversity and composition of the honeybee gut microbiome. To ensure comparability, all datasets were processed using a uniform preprocessing and analysis pipeline. In the studies considered, treatments were administered via feeding, with controls provided either as sugar water (sugar controls [ 19 ]) or as field controls, depending on the study. Comparisons were made between treated samples and their corresponding controls. Oxalic acid treatment produced the most pronounced shift in beta diversity (PERMANOVA, R² = 0.197, p = 0.001), while alpha diversity did not differ significantly from sugar controls ( Fig. 4A ). Taxonomic analysis using ANCOM-BC2 revealed that oxalic acid was associated with significant decreases in Bombella intestini (logFC < -3.8, adjusted p < 0.001) , Bartonella (logFC < -1.4, adjusted p = 0.003) , Apilactobacillus (logFC 1.2, adjusted p = 0.013) ( Fig. 4E , Table S8). These findings are consistent with the original study [ 19 ], which employed the ALDEx2 method [ 71 ] and detected similar taxa, although additional species were reported. Download figure Open in new tab Figure 4. Effects of treatments on alpha diversity, community composition, and differential abundance of specific taxa. (A–D) Changes in the honeybee gut microbiome following various treatments. Alpha diversity (Shannon index) is presented as boxplots, where the central line indicates the median, the box represents the interquartile range, and whiskers extend to the 10th and 90th percentiles. Group differences in alpha diversity were assessed using the Kruskal-Wallis test (p-values indicated). Beta diversity was evaluated using Bray-Curtis distances and visualized via PCoA, with ellipses representing 95% confidence intervals; PERMANOVA (R² and p-values) was used to assess differences between groups. The panels depict: (A) Oxalic acid (vs. sugar controls); (B) Glyphosate treatment in adult honeybees from Uruguay and Texas, USA, as well as newly emerged bees in Texas (vs. sugar controls); (C) Manuka honey (vs. field controls); and (D) Inverted sugar (vs. field controls). (E) Differentially abundant taxa identified by ANCOM-BC2 following treatment. Log₂ fold changes (logFC) are plotted on the x-axis, and only taxa with statistically significant differences (adjusted p-value < 0.05) are shown. Taxa names in bold indicate those that passed the sensitivity analysis for pseudo-count addition. Squared taxa belong to the core, round taxa are non-core. In contrast, treatments with other insecticides, including acetamiprid, thiacloprid, and imidacloprid did not lead to significant shifts in overall microbiome composition (Fig. S6), consistent with previous studies [ 4 , 19 , 25 ]. Nonetheless, differential abundance analysis identified specific ASVs that were altered by these insecticides, including increased abundances of Bombella apis, Hafnia-Obesumbacterium, Klebsiella , Lonsdalea, and Serratia , and a decrease in Fructobacillus and Apilactobacilluse ( Fig. 4E , Table S8). We did not detect significant composition differences for other insecticides evaluated in the same study [ 19 ] (Fig. S6). For glyphosate, we analyzed data from studies conducted in Uruguay [ 72 ] and Texas, USA [ 29 ]. In the Texas study, two experimental setups were performed, one treating adult bees and another treating newly emerged bees (NEW samples), while in the Uruguay study, only NEW samples were considered. Beta diversity analysis revealed a significant shift in Texas NEW samples (PERMANOVA, R 2 = 0.238, p = 0.001). In contrast, no significant changes were detected in adult bees from Texas or NEWs from Uruguay, in agreement with the original reports [ 4 , 29 ]. Although the Uruguay study reported compositional changes in NEW bees, these findings were based on comparisons involving multiple treatments rather than glyphosate alone (Table S9). Notably, ANCOM-BC2 analysis showed consistent reductions in Snodgrasella alvi in both the Texas and the Uruguay NEW samples ( Fig. 4E , Table S8) [ 4 , 29 ]. In adult bees from Texas, glyphosate treatment also resulted in a decrease in Fructobacillus fructosus and an increase in Lonsdalea , whereas the Uruguay NEW revealed reductions in Commensalibacter and S. alvi coupled with an increase in Staphylococcus . In contrast to pesticide treatments, dietary supplements such as manuka honey and inverted sugar significantly decreased alpha diversity relative to field controls, consistent with previous findings from New Zealand [ 18 ]. Both supplements induced shifts in microbiome composition as determined by beta diversity analysis, though the effect sizes were smaller than those observed with oxalic acid and glyphosate in Texas NEW samples (PERMANOVA R 2 = 0.091, p = 0.003 for manuka honey; R 2 = 0.095, p = 0.003 for inverted sugar) ( Fig. 4C, D ). At the taxonomic level, manuka honey supplementation was associated with a decrease in Apilactobacillus , a finding not reported in the original study ( Fig. 4E , Table S8, S9). Additionally, we did not identify any differentially abundant taxa with inverted sugar treatment, in contrast to the original study, which utilized ANOVA for statistical analysis [ 18 ]. Collectively, these results highlight that different treatments distinctly modulate the honeybee gut microbiome, with some (e.g., oxalic acid) eliciting robust compositional shifts and others (e.g., certain insecticides) causing more subtle, taxon-specific changes. However, given that focusing solely on individual taxa may not fully capture the complexity of the microbial community’s adaptations, we next employ co-abundance network analysis to elucidate collective effects and functional responses to these treatments. Co-Abundance Networks Reveal Adaptive Pesticide Responses and Potential Glyphosate Degradation in Honeybee Gut Microbiome To capture community-level responses to pesticide and sugar supplement treatments, we used the microbiome abundance data of all the studies and constructed a co-abundance network where each taxon is represented as a node and edges denote significant positive or negative correlations in their relative abundances. We then identified the subnetwork in which the sum of log fold change differential abundance weights (compared to controls) was greatest, thereby pinpointing the group of taxa that collectively respond most strongly to the treatment. Finally, we performed functional enrichment analysis on this subnetwork to determine which pathways and functional profiles may enable these bacteria to adapt and increase in abundance ( Fig. 5 and Fig. S6; see Methods for details). Download figure Open in new tab Figure 5. Co-abundance networks of the honeybee gut microbiome under treatments and corresponding KEGG pathway overrepresentation analysis. Panels depict treatment-associated subnetworks for: (A) glyphosate in Uruguay hive samples, (B) glyphosate in US Texas NEW samples, (C) oxalic acid, (D) acetamiprid, and (E) thiacloprid. For each treatment, the upper panel shows the co-abundance network and the lower panel displays the corresponding functional enrichment results. In the networks, nodes represent taxa and edges denote significant Spearman correlations (ρ > 0.6, adjusted p < 0.2); solid lines indicate positive correlations, and dashed lines negative correlations. Core taxa are represented by rectangular nodes, and non-core taxa by circles. Node size reflects the absolute log fold change (LFC) between treatment and control groups, with a green-to-blue color scale indicating up- and down-regulation, respectively. Nodes for which LFC values were not calculated by ANCOM-BC2 due to sample size limitations are colored gray. Nodes and edges highlighted with red borders indicate the treatment-specific subnetwork identified using the maximum weight connected subgraph algorithm (see Methods). Only KEGG pathways significantly enriched in the treatment groups (FDR < 0.1) are shown. The glyphosate-associated subgraph identified in the Uruguay and Texas NEW samples included interconnected taxa, with decreases in S. alvi in Uruguay and Texas, and decrease in Serratia in Uruguay. In the Uruguay NEW samples, this subnetwork expanded to encompass additional taxa: Frischella , Commensalibacter , Lactobacillus , Bifidobacterium asteroides , and Bombilactobacillus (all decreased), along with increases in Gilliamella , Fructobacillus (other), Lactobacillus apis and Pseudomonas . Pathway enrichment within this glyphosate-associated subnetwork in Uruguay identified several amino acid metabolic pathways, such as tryptophan and tyrosine metabolism, suggesting a direct influence of glyphosate exposure ( Fig. 5A, B , Table S12). Furthermore, our analysis revealed an enrichment of pathways for xenobiotic and aromatic compound degradation, nitrogen cycling, and two-component system, suggesting potential capacity for glyphosate degradation. Pseudomonas is known to degrade glyphosate either by converting it to sarcosine via the C-P lyase pathway or by transforming it into aminomethylphosphonic acid (AMPA) and methylamine [ 73 – 76 ]. In our samples, we observed enrichment of the phosphonate and phosphinate metabolism pathway, including the phnJ gene, which encodes C-P lyase; a key enzyme in both degradation routes. This gene is part of the phn operon, which also encodes putative glyphosate transporters ( phnC , phnD , and phnE ) [ 77 ], indicating glyphosate uptake capacity. Additionally, enriched pathways for glycine metabolism, nitrogen cycling, and aromatic compound degradation may further facilitate the breakdown of glyphosate degradation intermediates. Although the precise degradation route remains unclear, the increased abundance of Pseudomonas and these functional enrichments support the hypothesis that glyphosate degradation is an adaptive response to glyphosate exposure. In Uruguay samples, glyphosate exposure altered bacterial chemotaxis and biofilm formation pathways, likely mediated by both Pseudomonas and Gilliamella . This response may enhance overall bacterial stress adaptation [ 10 ]. Furthermore, a disconnected subnetwork in these samples showed upregulation of Staphylococcus and Streptococcus , aligning with previous reports of impaired immune function following glyphosate treatment [ 78 ]. In contrast, treatments with oxalic acid and neonicotinoids (acetamiprid and thiacloprid) affected a different subnetwork that included Klebsiella , Hafnia-Obesumbacterium , and Serratia . Klebsiella and Serratia increased in both cases, while Hafnia-Obesumbacterium increased under neonicotinoid exposure, and decreased under oxalic acid treatment. ( Fig. 5C ). For neonicotinoids, this subnetwork further expanded to include Fructobacillus , which decreased in abundance ( Fig. 5D,E ). Oxalic acid’s impact also extended to include decreased abundances of Snodgrassella alvi, and Commensalibacter , highlighting the documented sensitivity of these taxa to pesticides [ 19 ]. Unlike glyphosate, we did not detect specific oxalic acid degradation pathways in the samples. Instead, the oxalic acid–affected subnetwork showed enrichment for stress response pathways. ABC transporters could enhance the regulation of internal pH after its reduction by oxalic acid application [ 79 ]. Upregulation of siderophore biosynthesis suggests a compensatory response to iron limitation caused by oxalic acid sequestration [ 80 ]. Similarly, neonicotinoid exposure commonly enriched pathways related to stress responses, such as the two-component system, which is typically involved in insecticide and xenobiotic detoxification [ 81 , 82 ]. Notably, the pathway for cationic antimicrobial peptide (CAMP) resistance was enriched in neonicotinoid-associated subnetworks, potentially reflecting increased antimicrobial peptide production by the host, as previously observed in bumblebees exposed to these insecticides [ 72 ] ( Fig. 5D-E , Table S12). Discussion A major strength of our study is the standardized re-analysis of diverse datasets. By using both 16S and WGS data and a unified pipeline, we minimize methodological discrepancies that may have contributed to conflicting findings in previous studies [ 2 , 34 , 35 ]. Our novel co-abundance network approach links compositional shifts with functional potential. We also provide a fully reproducible pipeline to facilitate future research across diverse biomes of interest, along with an R package, micronet, for co-abundance network analysis. Nevertheless, we acknowledge several limitations. Despite efforts to normalize library size and limit primer biases, residual heterogeneity remained across studies, especially due to differing extraction methods, primers, and sequencing approaches. This, however, is an inherent issue of any meta-analysis. Additionally, variations in laboratory and field pesticide exposures, cross-sectional sampling designs, and the inherent limitations of ASV-based annotations (which miss strain-level differences) constrain our interpretations. Notably, the WGS data were available only for untreated samples, limiting strain-level functional inference in treated groups. Our results support the consensus that honeybees maintain a defined core gut microbiome, typically including six phylotypes: Lactobacillus, Bombilactobacillus, Gilliamella, Snodgrassella, Frischella , and Bifidobacterium . Notably, incorporating low-abundance taxa via WGS datasets suggests an expanded core of up to twelve taxa, including Bartonella, Commensalibacter, Frischella, Apilactobacillus, Moraxella, and Vibrio . Studies show that honeybee microbiota are shaped by diet, climate, and local flora [ 67 , 68 , 83 – 85 ]. We also observed significant geographic variation. However, the factors driving this variation are challenging to disentangle. Seasonal factors and urbanization likely play a role, as evidenced by our newly generated data from Armenia. Other biological factors unaccounted for in our study, such as bee age, temperature, floral resources, and colony management, may also contribute to these differences. This may explain why different studies often conflict; for instance, our finding of elevated Commensalibacter in autumn agrees with a Swiss study but contrasts with reports from other European and Asian sites [ 64 , 66 ]. Moreover, variations in data analysis pipelines have contributed to inconsistent results, which underscores the value of our unified reprocessing of both 16S and WGS datasets from multiple studies. A pressing question in bee research is how agrochemicals affect colony health. Recent studies link insecticides, and herbicides to shifts in gut community composition, with conflicting or inconclusive findings regarding affected taxa and functions [ 4 , 19 , 23 , 25 ]. Our meta-analysis supports that gut responses vary by chemical class, dose, and bee life stage. Glyphosate affected the core taxa. We consistently observed depletion of Snodgrassella alvi confirming its reported sensitivity. In Texas, newly emerged bees showed marked beta diversity shifts, and in Uruguayan newly emerged bees exhibited drastic network changes suggesting higher vulnerability in younger bees. Adult bees in Texas displayed milder shifts, hinting at a more resilient microbiome in older foragers [ 4 , 29 ]. Co-abundance network analysis revealed potential adaptation mechanisms. In Uruguayan newly emerged bees, a glyphosate-affected subnetwork involving increased Pseudomonas, Gilliamella, and enrichment of several pathways, suggesting a potential for C-P lyase-mediated glyphosate degradation. While Pseudomonas species have demonstrated glyphosate utilization in other environments [ 76 , 77 ], this, to the best of our knowledge, is the first indication of such an adaptation mechanism in honeybees. Enrichment of tryptophan and tyrosine metabolism pathways further suggests compensation for glyphosate-induced shikimate pathway inhibition. In contrast, the Texas newly emerged bees exhibited a smaller affected subnetwork. Strain-level differences or variations in experimental setup could affect these differences. The Texas bees were exposed to a higher glyphosate dose (169 mg/l versus 10 mg/l in Uruguay), however at a shorter period (2 days versus 7 days in Uruguay). This also suggests that prolonged exposures have a stronger impact on microbiome shifts and adaptations than higher doses. Insecticides, on the other hand, did not affect the core taxa but induced an increased abundance of opportunistic pathogens that exist as part of the normal microbiota. Neonicotinoids (imidacloprid, thiacloprid, acetamiprid) induced subtler microbiome alterations [ 19 ]. Although overall community composition remained relatively stable, low-abundance opportunistic pathogens, such as Serratia , Klebsiella , and Hafnia-Obesumbacterium thrived under treatment. These taxa were enriched for pathways of xenobiotic detoxification and stress response (e.g., two-component system, ABC transporters), and cationic antimicrobial peptide (CAMP) resistance. These adaptations could confer a selective advantage if the host produces more antimicrobial peptides under pesticide-induced stress [ 72 ]. In contrast, the acaricide oxalic acid caused the most pronounced beta-diversity shift, and also increased Serratia , and Klebsiella , similar to neonicotinoids. Additionally, it led to the depletion of Bartonella, Bombella, Apilactobacillus , and Snodgrassella [ 19 , 80 ]. Although considered relatively “bee-safe,” oxalic acid lowers gut pH, imposing selective pressure that appears to upregulate general pathways of response to low pH and iron limitation (e.g., ABC transporters, siderophore biosynthesis), rather than dedicated oxalate-degradation mechanisms [ 80 , 86 ]. All in all, our meta-analysis reaffirms a stable honeybee core microbiota while demonstrating that environmental factors (geography, season) and anthropogenic influences (pesticides, dietary additives) induce taxon-specific shifts. Pesticides tend to favor bacteria with enhanced stress-response capabilities or opportunistic degraders. Longitudinal metagenomics and mechanistic studies remain essential to clarify how these microbiome shifts impact colony health, ultimately guiding sustainable strategies for apiculture and pollinator conservation. Data and software availability Raw 16S rRNA gene sequencing and WGS data sets evaluated in this study are available from the NCBI Sequence Read Archive (PRJNA719169, PRJNA732842, PRJNA775827, PRJNA432211, PRJNA432210, PRJNA483763, PRJNA531038, PRJNA685398, PRJNA787435 and PRJNA494922). The 16s amplicon datasets from the Armenian samples generated in this study are deposited in the NCBI SRA database (PRJNA1220608). The micronet package for co-abundance network analysis is available for download and use at Github: https://github.com/abi-am/micronet . The rest of the scripts used in this study are available at https://github.com/abi-am/bee-microbiome-metaanalysis . Contributions LN and NV conceived the study. NV performed the analyses, produced the figures, and wrote the manuscript. IB generated the Armenian data. LA developed the micronet package and conducted the co-abundance network analysis. MZ and HL carried out additional experiments. LN supervised the study and prepared the manuscript. CM supervised computational analysis and revised the manuscript. All authors read and approved the final version of the manuscript. Acknowledgements The authors acknowledge funding from the Higher Education and Science Committee (HESC) MESCS RA PhD Support program awarded to NV (24AA-1F065), the HESC RA Prospective Directions grant awarded to LN (24FP-2I061), and Yerevan State University Inner Grant 2022 to IB. References 1. ↵ Hung KLJ , Kingston JM , Albrecht M et al. The worldwide importance of honey bees as pollinators in natural habitats . Proceedings of the Royal Society B: Biological Sciences 2018 ; 285 , DOI: 10.1098/RSPB.2017.2140 . OpenUrl CrossRef 2. ↵ Romero S , Nastasa A , Chapman A et al. The honey bee gut microbiota: strategies for study and characterization . Insect Mol Biol 2019 ; 28 : 455 – 72 . OpenUrl CrossRef PubMed 3. ↵ Klein AM , Vaissière BE , Cane JH et al. Importance of pollinators in changing landscapes for world crops . Proceedings of the Royal Society B: Biological Sciences 2007 ; 274 : 303 – 13 . OpenUrl CrossRef PubMed Web of Science 4. ↵ Castelli L , Balbuena S , Branchiccela B et al. Impact of chronic exposure to sublethal doses of glyphosate on honey bee immunity, gut microbiota and infection by pathogens . Microorganisms 2021 ; 9 : 845 . OpenUrl CrossRef PubMed 5. Goulson D , Nicholls E , Botías C et al. Bee declines driven by combined Stress from parasites, pesticides, and lack of flowers . Science (1979) 2015 ; 347 , DOI: 10.1126/SCIENCE.1255957/ASSET/25ACDBB7-03E0-4D87-BEC1-0935B850C29F/ASSETS/GRAPHIC/347_1255957_FA.JPEG . OpenUrl CrossRef 6. ↵ Langowska A , Zawilak M , Sparks TH et al. Long-term effect of temperature on honey yield and honeybee phenology . , DOI: 10.1007/s00484-016-1293-x . OpenUrl CrossRef 7. ↵ Colin T , Meikle WG , Paten AM et al. Long-term dynamics of honey bee colonies following exposure to chemical stress . Science of The Total Environment 2019 ; 677 : 660 – 70 . OpenUrl CrossRef PubMed 8. ↵ Owen R . Role of Human Action in the Spread of Honey Bee (Hymenoptera: Apidae) Pathogens . J Econ Entomol 2017 ; 110 : 797 – 801 . OpenUrl CrossRef PubMed 9. ↵ Cameron SA , Lozier JD , Strange JP et al. Patterns of widespread decline in North American bumble bees . Proc Natl Acad Sci U S A 2011 ; 108 : 662 – 7 . OpenUrl Abstract / FREE Full Text 10. ↵ Motta EVS , de Jong TK , Gage A et al. Glyphosate effects on growth and biofilm formation in bee gut symbionts and diverse associated bacteria . Appl Environ Microbiol 2024 ; 90 , DOI: 10.1128/AEM.00515-24/SUPPL_FILE/AEM.00515-24-S0004.XLSX . OpenUrl CrossRef 11. ↵ Kwong WK , Moran NA . Gut microbial communities of social bees . Nature Reviews Microbiology 2016 14:6 2016 ; 14 : 374 – 84 . OpenUrl CrossRef PubMed 12. Amiri N , Keady MM , Lim HC . Honey bees and bumble bees occupying the same landscape have distinct gut microbiomes and amplicon sequence variant-level responses to infections . PeerJ 2023 ; 11 : e15501 . OpenUrl CrossRef PubMed 13. Motta EVS , Moran NA . The honeybee microbiota and its impact on health and disease . Nature Reviews Microbiology 2023 22:3 2023 ; 22 : 122 – 37 . OpenUrl PubMed 14. ↵ Martinson VG , Danforth BN , Minckley RL et al. A simple and distinctive microbiota associated with honey bees and bumble bees . Mol Ecol 2011 ; 20 : 619 – 28 . OpenUrl CrossRef PubMed Web of Science 15. ↵ Raymann K , Moran NA . The role of the gut microbiome in health and disease of adult honey bee workers . Curr Opin Insect Sci 2018 ; 26 : 97 – 104 . OpenUrl CrossRef PubMed 16. ↵ Zheng H , Steele MI , Leonard SP et al. Honey bees as models for gut microbiota research . Lab Animal 2018 47:11 2018 ; 47 : 317 – 25 . OpenUrl CrossRef PubMed 17. ↵ Corby-Harris V , Maes P , Anderson KE . The Bacterial Communities Associated with Honey Bee (Apis mellifera) Foragers . PLoS One 2014 ; 9 : e95056 . OpenUrl CrossRef PubMed 18. ↵ Taylor MA , Robertson AW , Biggs PJ et al. The effect of carbohydrate sources: Sucrose, invert sugar and components of mānuka honey, on core bacteria in the digestive tract of adult honey bees (Apis mellifera) . PLoS One 2019 ; 14 : e0225845 . OpenUrl CrossRef PubMed 19. ↵ Cuesta-Maté A , Renelies-Hamilton J , Kryger P et al. Resistance and Vulnerability of Honeybee (Apis mellifera) Gut Bacteria to Commonly Used Pesticides . Front Microbiol 2021 ; 12 : 717990 . OpenUrl CrossRef PubMed 20. Mullin CA , Frazier M , Frazier JL et al. High Levels of Miticides and Agrochemicals in North American Apiaries: Implications for Honey Bee Health . PLoS One 2010 ; 5 : e9754 . OpenUrl CrossRef PubMed 21. ↵ Traynor KS , Pettis JS , Tarpy DR et al. In-hive Pesticide Exposome: Assessing risks to migratory honey bees from in-hive pesticide contamination in the Eastern United States . Scientific Reports 2016 6:1 2016 ; 6 : 1 – 16 . OpenUrl CrossRef PubMed 22. ↵ Casida JE . Neonicotinoids and Other Insect Nicotinic Receptor Competitive Modulators: Progress and Prospects . Annu Rev Entomol 2018 ; 63 : 125 – 44 . OpenUrl CrossRef PubMed 23. ↵ Liu YJ , Qiao NH , Diao QY et al. Thiacloprid exposure perturbs the gut microbiota and reduces the survival status in honeybees . J Hazard Mater 2020 ; 389 , DOI: 10.1016/J.JHAZMAT.2019.121818 . OpenUrl CrossRef 24. ↵ Khoury S El , Gauthier J , Bouslama S et al. Dietary contamination with a neonicotinoid (Clothianidin) gradient triggers specific dysbiosis signatures of microbiota activity along the honeybee (apis mellifera) digestive tract . Microorganisms 2021 ; 9 : 2283 . OpenUrl CrossRef PubMed 25. ↵ Raymann K , Motta EVS , Girard C et al. Imidacloprid decreases honey bee survival rates but does not affect the gut microbiome . Appl Environ Microbiol 2018 ; 84 , DOI: 10.1128/AEM.00545-18/SUPPL_FILE/ZAM013188579S1.PDF . OpenUrl CrossRef 26. ↵ Rademacher E , Harz M , Schneider S . Effects of Oxalic Acid on Apis mellifera (Hymenoptera: Apidae) . Insects 2017 , Vol 8, Page 84 2017 ; 8 : 84 . OpenUrl CrossRef PubMed 27. ↵ Farina WM , Balbuena MS , Herbert LT et al. Effects of the Herbicide Glyphosate on Honey Bee Sensory and Cognitive Abilities: Individual Impairments with Implications for the Hive . Insects 2019 , Vol 10, Page 354 2019 ; 10 : 354 . OpenUrl CrossRef PubMed 28. ↵ Rainio MJ , Ruuskanen S , Helander M et al. Adaptation of bacteria to glyphosate: a microevolutionary perspective of the enzyme 5-enolpyruvylshikimate-3-phosphate synthase . Environ Microbiol Rep 2021 ; 13 : 309 – 16 . OpenUrl CrossRef PubMed 29. ↵ Motta EVS , Raymann K , Moran NA . Glyphosate perturbs the gut microbiota of honey bees . Proc Natl Acad Sci U S A 2018 ; 115 : 10305 – 10 . OpenUrl Abstract / FREE Full Text 30. ↵ Motta EVS , Moran NA . Impact of Glyphosate on the Honey Bee Gut Microbiota: Effects of Intensity, Duration, and Timing of Exposure . mSystems 2020 ; 5 , DOI: 10.1128/MSYSTEMS.00268-20/SUPPL_FILE/REVIEWER-COMMENTS.PDF . OpenUrl CrossRef 31. ↵ Alberoni D , Favaro R , Baffoni L et al. Neonicotinoids in the agroecosystem: In-field long-term assessment on honeybee colony strength and microbiome . Science of The Total Environment 2021 ; 762 : 144116 . OpenUrl CrossRef PubMed 32. ↵ Almasri H , Liberti J , Brunet JL et al. Mild chronic exposure to pesticides alters physiological markers of honey bee health without perturbing the core gut microbiota . Scientific Reports 2022 12:1 2022 ; 12 : 1 – 15 . OpenUrl CrossRef PubMed 33. ↵ Blot N , Veillat L , Rouzé R et al. Glyphosate, but not its metabolite AMPA, alters the honeybee gut microbiota . PLoS One 2019 ; 14 : e0215466 . OpenUrl CrossRef PubMed 34. ↵ Kim D , Hofstaedter CE , Zhao C et al. Optimizing methods and dodging pitfalls in microbiome research . Microbiome 2017 5:1 2017 ; 5 : 1 – 14 . OpenUrl CrossRef PubMed 35. ↵ Hrovat K , Dutilh BE , Medema MH et al. Taxonomic resolution of different 16S rRNA variable regions varies strongly across plant-associated bacteria . ISME Communications 2024 ; 4 : 34 . OpenUrl 36. ↵ Powell JE , Martinson VG , Urban-Mead K et al. Routes of acquisition of the gut microbiota of the honey bee Apis mellifera . Appl Environ Microbiol 2014 ; 80 : 7378 – 87 . OpenUrl Abstract / FREE Full Text 37. ↵ Callahan BJ , McMurdie PJ , Rosen MJ et al. DADA2: High-resolution sample inference from Illumina amplicon data . Nat Methods 2016 ; 13 : 581 – 3 . OpenUrl CrossRef PubMed 38. ↵ McMurdie PJ , Holmes S . Phyloseq: A bioconductor package for handling and analysis of high-throughput phylogenetic sequence data . Pacific Symposium on Biocomputing 2012 : 235 – 46 . 39. ↵ Bodenhofer U , Bonatesta E , Horejš-Kainrath C et al. msa: an R package for multiple sequence alignment . Bioinformatics 2015 ; 31 : 3997 – 9 . OpenUrl CrossRef PubMed 40. ↵ Glöckner FO . The SILVA Database Project: An ELIXIR core data resource for high-quality ribosomal RNA sequences . Biodiversity Information Science and Standards 3 : e36125 2019 ; 3 : e36125 -. OpenUrl CrossRef 41. ↵ Wood DE , Lu J , Langmead B . Improved metagenomic analysis with Kraken 2 . Genome Biol 2019 ; 20 : 1 – 13 . OpenUrl CrossRef PubMed 42. ↵ McMurdie PJ , Holmes S. phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data . PLoS One 2013 ; 8 : e61217 . OpenUrl CrossRef PubMed 43. ↵ Barnett DJ m., Arts IC w., Penders J . microViz: an R package for microbiome data visualization and statistics . J Open Source Softw 2021 ; 6 : 3201 . OpenUrl CrossRef 44. ↵ Peddada S , Lin H. Multi-group Analysis of Compositions of Microbiomes with Covariate Adjustments and Repeated Measures . 2023 , DOI: 10.21203/RS.3.RS-2778207/V1 . OpenUrl CrossRef 45. ↵ Shannon P , Markiel A , Ozier O et al. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks . Genome Res 2003 ; 13 : 2498 – 504 . OpenUrl Abstract / FREE Full Text 46. ↵ Uritskiy G V. , Diruggiero J , Taylor J . MetaWRAP - A flexible pipeline for genome-resolved metagenomic data analysis 08 Information and Computing Sciences 0803 Computer Software 08 Information and Computing Sciences 0806 Information Systems . Microbiome 2018 ; 6 : 1 – 13 . OpenUrl CrossRef PubMed 47. ↵ Chaumeil PA , Mussig AJ , Hugenholtz P et al. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database . Bioinformatics 2020 ; 36 : 1925 – 7 . OpenUrl CrossRef 48. ↵ Douglas GM , Maffei VJ , Zaneveld JR et al. PICRUSt2 for prediction of metagenome functions . Nature Biotechnology 2020 38:6 2020 ; 38 : 685 – 8 . OpenUrl CrossRef PubMed 49. ↵ Huerta-Cepas J , Szklarczyk D , Forslund K et al. EGGNOG 4.5: A hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences . Nucleic Acids Res 2016 ; 44 : D286 – 93 . OpenUrl CrossRef PubMed 50. ↵ Yu G , Wang LG , Han Y et al. ClusterProfiler: An R package for comparing biological themes among gene clusters . OMICS 2012 ; 16 : 284 – 7 . OpenUrl CrossRef PubMed Web of Science 51. ↵ Bonilla-Rosso G , Engel P . Functional roles and metabolic niches in the honey bee gut microbiota . Curr Opin Microbiol 2018 ; 43 : 69 – 76 . OpenUrl CrossRef PubMed 52. ↵ Ellegaard KM , Engel P . Genomic diversity landscape of the honey bee gut microbiota . Nature Communications 2019 10:1 2019 ; 10 : 1 – 13 . OpenUrl CrossRef PubMed 53. ↵ Motta EVS , Moran NA . The honeybee microbiota and its impact on health and disease . Nat Rev Microbiol 2024 ; 22 : 122 – 37 . OpenUrl CrossRef PubMed 54. ↵ Moran NA , Hansen AK , Powell JE et al. Distinctive Gut Microbiota of Honey Bees Assessed Using Deep Sampling from Individual Worker Bees . PLoS One 2012 ; 7 : e36393 . OpenUrl CrossRef PubMed 55. ↵ Rouzé R , Moné A , Delbac F et al. The Honeybee Gut Microbiota Is Altered after Chronic Exposure to Different Families of Insecticides and Infection by Nosema ceranae . Microbes Environ 2019 ; 34 : 226 – 33 . OpenUrl CrossRef PubMed 56. ↵ Wang J , Su Q , Zhang X et al. Entomomonas moraniae gen. Nov., sp. nov., a member of the family Pseudomonadaceae isolated from asian honey bee gut, possesses a highly reduced genome . Int J Syst Evol Microbiol 2020 ; 70 : 165 – 71 . OpenUrl CrossRef PubMed 57. ↵ Kačániová M , Terentjeva M , Žiarovská J et al. In Vitro Antagonistic Effect of Gut Bacteriota Isolated from Indigenous Honey Bees and Essential Oils against Paenibacillus Larvae . International Journal of Molecular Sciences 2020 , Vol 21, Page 6736 2020 ; 21 : 6736 . OpenUrl CrossRef PubMed 58. ↵ Laconi A , Tolosi R , Mughini-Gras L et al. Beehive products as bioindicators of antimicrobial resistance contamination in the environment . Science of the Total Environment 2022 ; 823 , DOI: 10.1016/j.scitotenv.2021.151131 . OpenUrl CrossRef 59. ↵ Smutin D , Lebedev E , Selitskiy M et al. Micro”bee”ota: Honey Bee Normal Microbiota as a Part of Superorganism . Microorganisms 2022 , Vol 10, Page 2359 2022 ; 10 : 2359 . OpenUrl CrossRef PubMed 60. Scheiner R , Strauß S , Thamm M et al. The Bacterium Pantoea ananatis Modifies Behavioral Responses to Sugar Solutions in Honeybees . Insects 2020 , Vol 11, Page 692 2020 ; 11 : 692 . OpenUrl CrossRef PubMed 61. ↵ Nguyen PN , Rehan SM . Wild bee and pollen microbiomes across an urban-rural divide . FEMS Microbiol Ecol 2023 ; 99 , DOI: 10.1093/FEMSEC/FIAD158 . OpenUrl CrossRef 62. Shell WA , Rehan SM . Comparative metagenomics reveals expanded insights into intra- and interspecific variation among wild bee microbiomes . Communications Biology 2022 5:1 2022 ; 5 : 1 – 12 . OpenUrl CrossRef PubMed 63. ↵ El Khoury S , Gauthier J , Mercier PL et al. Honeybee gut bacterial strain improved survival and gut microbiota homeostasis in Apis mellifera exposed in vivo to clothianidin . Microbiol Spectr 2024 , DOI: 10.1128/SPECTRUM.00578-24/SUPPL_FILE/REVIEWER-COMMENTS.PDF . OpenUrl CrossRef 64. ↵ Kešnerová L , Emery O , Troilo M et al. Gut microbiota structure differs between honeybees in winter and summer . ISME J 2020 ; 14 : 801 – 14 . OpenUrl CrossRef PubMed 65. ↵ Jabal-Uriel C , Alba C , Higes M et al. Effect of Nosema ceranae infection and season on the gut bacteriome composition of the European honeybee (Apis mellifera) . Scientific Reports 2022 12:1 2022 ; 12 : 1 – 13 . OpenUrl CrossRef PubMed 66. ↵ Li C , Tang M , Li X et al. Community Dynamics in Structure and Function of Honey Bee Gut Bacteria in Response to Winter Dietary Shift . mBio 2022 ; 13 , DOI: 10.1128/MBIO.01131-22/SUPPL_FILE/MBIO.01131-22-S0005.PDF . OpenUrl CrossRef 67. ↵ Luo S , Zhang X , Zhou X . Temporospatial dynamics and host specificity of honeybee gut bacteria . Cell Rep 2024 ; 43 : 114408 . OpenUrl CrossRef PubMed 68. ↵ Muñoz-Colmenero M , Baroja-Careaga I , Kovačić M et al. Differences in honey bee bacterial diversity and composition in agricultural and pristine environments – a field study . Apidologie 2020 ; 51 : 1018 – 37 . OpenUrl CrossRef 69. ↵ Hénaff E , Najjar D , Perez M , et al. Holobiont Urbanism: sampling urban beehives reveals cities’ metagenomes . Environ Microbiome 2023 ; 18 : 1 – 12 . OpenUrl CrossRef PubMed 70. ↵ Nguyen PN , Rehan SM . The effects of urban land use gradients on wild bee microbiomes . Front Microbiol 2022 ; 13 : 992660 . OpenUrl CrossRef PubMed 71. ↵ Fernandes AD , Macklaim JM , Linn TG et al. ANOVA-Like Differential Expression (ALDEx) Analysis for Mixed Population RNA-Seq . PLoS One 2013 ; 8 : e67019 . OpenUrl CrossRef PubMed 72. ↵ Simmons WR , Angelini DR . Chronic exposure to a neonicotinoid increases expression of antimicrobial peptide genes in the bumblebee Bombus impatiens . Scientific Reports 2017 7:1 2017 ; 7 : 1 – 9 . OpenUrl CrossRef PubMed 73. ↵ Zhao H , Tao K , Zhu J et al. Bioremediation potential of glyphosate-degrading pseudomonas spp. Strains isolated from contaminated soil . Journal of General and Applied Microbiology 2015 ; 61 : 165 – 70 . OpenUrl 74. Zhan H , Feng Y , Fan X et al. Recent advances in glyphosate biodegradation . Appl Microbiol Biotechnol 2018 ; 102 : 5033 – 43 . OpenUrl CrossRef 75. Feng D , Soric A , Boutin O . Treatment technologies and degradation pathways of glyphosate: A critical review . Science of The Total Environment 2020 ; 742 : 140559 . OpenUrl CrossRef PubMed 76. ↵ Singh S , Kumar V , Gill JPK et al. Herbicide Glyphosate: Toxicity and Microbial Degradation . International Journal of Environmental Research and Public Health 2020, Vol 17, Page 7519 2020 ; 17 : 7519 . OpenUrl CrossRef 77. ↵ Hove-Jensen B , Zechel DL , Jochimsen B . Utilization of Glyphosate as Phosphate Source: Biochemistry and Genetics of Bacterial Carbon-Phosphorus Lyase . Microbiology and Molecular Biology Reviews 2014 ; 78 : 176 – 97 . OpenUrl Abstract / FREE Full Text 78. ↵ Motta EVS , Powell JE , Moran NA . Glyphosate induces immune dysregulation in honey bees . Anim Microbiome 2022 ; 4 : 1 – 14 . OpenUrl CrossRef PubMed 79. ↵ Wilks JC , Kitko RD , Cleeton SH et al. Acid and base stress and transcriptomic responses in Bacillus subtilis . Appl Environ Microbiol 2009 ; 75 : 981 – 90 . OpenUrl Abstract / FREE Full Text 80. ↵ Grąz M . Role of oxalic acid in fungal and bacterial metabolism and its biotechnological potential . World J Microbiol Biotechnol 2024 ; 40 : 1 – 11 . OpenUrl CrossRef 81. ↵ Lv Y , Li J , Yan K et al. Functional characterization of ABC transporters mediates multiple neonicotinoid resistance in a field population of Aphis gossypii Glover . Pestic Biochem Physiol 2022 ; 188 : 105264 . OpenUrl CrossRef PubMed 82. ↵ Wu C , Chakrabarty S , Jin M et al. Insect ATP-Binding Cassette (ABC) Transporters: Roles in Xenobiotic Detoxification and Bt Insecticidal Activity . International Journal of Molecular Sciences 2019 , Vol 20, Page 2829 2019 ; 20 : 2829 . OpenUrl CrossRef PubMed 83. ↵ Almeida EL , Ribiere C , Frei W et al. Geographical and Seasonal Analysis of the Honeybee Microbiome . Microb Ecol 2023 ; 85 : 765 – 78 . OpenUrl CrossRef PubMed 84. Su Q , Tang M , Hu J et al. Significant compositional and functional variation reveals the patterns of gut microbiota evolution among the widespread Asian honeybee populations . Front Microbiol 2022 ; 13 : 934459 . OpenUrl CrossRef PubMed 85. ↵ Oliveira Soares K , Ferreira T , Rocha D et al. Comparing the impact of landscape on the gut microbiome of Apis mellifera in Atlantic Forest and Caatinga Biomes . Scientific Reports 2025 15:1 2025 ; 15 : 1 – 9 . OpenUrl CrossRef PubMed 86. ↵ Liu X , Zhang K , Liu Y et al. Oxalic Acid From Sesbania rostrata Seed Exudates Mediates the Chemotactic Response of Azorhizobium caulinodans ORS571 Using Multiple Strategies . Front Microbiol 2019 ; 10 : 490471 . OpenUrl 87. Raymann K , Coon KL , Shaffer Z et al. Pathogenicity of serratia marcescens strains in honey bees . mBio 2018 ; 9 , DOI: 10.1128/MBIO.01649-18/SUPPL_FILE/MBO005184101SF7.PDF . OpenUrl CrossRef 88. Sun H , Mu X , Zhang K et al. Geographical resistome profiling in the honeybee microbiome reveals resistance gene transfer conferred by mobilizable plasmids . Microbiome 2022 ; 10 : 1 – 14 . OpenUrl CrossRef PubMed 89. Regan T , Barnett MW , Laetsch DR et al. Characterisation of the British honey bee metagenome . Nature Communications 2018 9:1 2018 ; 9 : 1 – 13 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted February 25, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Global distribution of honeybee gut microbiome and pesticide-driven adaptations in opportunistic microbial species Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Global distribution of honeybee gut microbiome and pesticide-driven adaptations in opportunistic microbial species Nelli Vardazaryan , Lusine Adunts , Inga Bazukyan , Magdalina Zakharyan , Honglian Liu , Chrats Melkonian , Lilit Nersisyan bioRxiv 2025.02.25.640077; doi: https://doi.org/10.1101/2025.02.25.640077 Share This Article: Copy Citation Tools Global distribution of honeybee gut microbiome and pesticide-driven adaptations in opportunistic microbial species Nelli Vardazaryan , Lusine Adunts , Inga Bazukyan , Magdalina Zakharyan , Honglian Liu , Chrats Melkonian , Lilit Nersisyan bioRxiv 2025.02.25.640077; doi: https://doi.org/10.1101/2025.02.25.640077 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Microbiology Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17691) Bioengineering (13892) Bioinformatics (41937) Biophysics (21452) Cancer Biology (18589) Cell Biology (25504) Clinical Trials (138) Developmental Biology (13378) Ecology (19899) Epidemiology (2067) Evolutionary Biology (24320) Genetics (15609) Genomics (22506) Immunology (17736) Microbiology (40394) Molecular Biology (17181) Neuroscience (88605) Paleontology (666) Pathology (2832) Pharmacology and Toxicology (4824) Physiology (7641) Plant Biology (15156) Scientific Communication and Education (2045) Synthetic Biology (4294) Systems Biology (9825) Zoology (2271)