Nested Connections: Local Phage and Broad Plasmid Sharing in the Honey Bee Mobilome

Nested Connections: Local Phage and Broad Plasmid Sharing in the Honey Bee Mobilome | 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 Nested Connections: Local Phage and Broad Plasmid Sharing in the Honey Bee Mobilome View ORCID Profile Delaney L. Miller , Lílian Caesar , Carrie Ganote , Sergio López-Madrigal , Danny W. Rice , Chris R. P. Robinson , Alyssa M. Seawright , Jillian P. Lewis , Irene L. G. Newton doi: https://doi.org/10.1101/2025.10.31.685865 Delaney L. Miller 1 University of Wisconsin Madison Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Delaney L. Miller Lílian Caesar 2 Indiana University Bloomington Find this author on Google Scholar Find this author on PubMed Search for this author on this site Carrie Ganote 2 Indiana University Bloomington Find this author on Google Scholar Find this author on PubMed Search for this author on this site Sergio López-Madrigal 2 Indiana University Bloomington Find this author on Google Scholar Find this author on PubMed Search for this author on this site Danny W. Rice 2 Indiana University Bloomington Find this author on Google Scholar Find this author on PubMed Search for this author on this site Chris R. P. Robinson 2 Indiana University Bloomington Find this author on Google Scholar Find this author on PubMed Search for this author on this site Alyssa M. Seawright 2 Indiana University Bloomington Find this author on Google Scholar Find this author on PubMed Search for this author on this site Jillian P. Lewis 3 University of Connecticut Find this author on Google Scholar Find this author on PubMed Search for this author on this site Irene L. G. Newton 2 Indiana University Bloomington Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: irnewton{at}iu.edu Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Honey bees rely on bacterial symbionts for their nutritional needs and for protection against invading pathogens. Genetic diversity among strains within the colony has the potential to impact symbiont function and subsequently the benefits that honey bees receive from them. Mobilized genes vectored by mobile genetic elements (MGEs) like phages, transposons, conjugative elements and plasmids, are known to rapidly alter bacterial phenotypes. We identified phages and plasmids in genomes of the symbiont Bombella apis between two colonies, with the goal of understanding which MGEs contribute to strain diversification as well as MGE distribution across colonies and between microbial species. Interestingly, we found some B. apis strains carry plasmids while all harbor a diversity of integrated phages, with only one phage clusters conserved across all. Identified B. apis phages are not found outside of the Bombella and Saccaribacter species, suggesting some host specificity for these MGEs. Of the five plasmids discovered, two appear to be phage-plasmids with high similarity to phages found in previously sequenced B. apis genomes. Interestingly, three plasmids in B. apis shared significant average nucleotide identity with known plasmids from acetic acid bacteria isolated from flowers, plants, and fermented foods. This result suggests that B. apis has acquired MGEs, either vertically or horizontally, from plant- and fermented-food associated AABs. Overall, our findings suggest that MGE content varies between colonies and has the potential to shape genetic and phenotypic variation between strains. Introduction Bacteria are ubiquitous among eukaryotic hosts and their role in shaping host health is a tenet of modern biology, and yet we still understand only a fraction of their genetic and functional diversity [ 1 , 2 ]. Pangenome analyses of bacteria, which aim to capture the breadth of genes that vary between strains within a species, show that symbionts with free-living life stages have highly variable or “open” pangenomes, meaning they are more variable due to gene gains and losses [ 3 ]. This fluidity can have consequences for host health. Strain-level variation among symbionts is known to affect disease outcomes, nutrition availability, and host development [ 4 – 6 ]. It is therefore important to use population-level approaches to identify gene gains and losses that may impact symbiont function. Among bacteria, mobile genetic elements (MGEs) such as plasmids, insertion sequences, integrative conjugative elements (ICEs), and phages facilitate rapid gene gain and loss [ 7 ]. Among symbionts, the most well-studied examples of mobilized traits can be found in plant and insect symbionts. A classic example is rhizobial symbionts, which cannot fix nitrogen for their legume hosts without the nitrogenases encoded on MGEs (plasmids or ICEs depending on the species) [ 8 – 10 ]. In insect symbionts, the gene clusters that encode for antimicrobials, which defend their host from disease, are often mobile, horizontally transmitting between species [ 4 , 11 ]. Another striking example of host traits mobilized by bacterial MGEs is from the symbiont Wolbachia pipientis and its prophage WO which carries a two-gene module responsible for cytoplasmic incompatibility, a mechanism for reproductive manipulation [ 12 ]. Clearly, symbionts and their myriad MGEs have the potential to radically change their host’s health and ecology, but MGEs remain understudied, both due to the challenges recovering them from genome assemblies and the difficulty classifying them, especially in non-model bacterial hosts [ 13 – 15 ]. Honey bees are agriculturally important pollinators and their microbiomes impact bee health and productivity [ 16 – 18 ]. Within the adult honey bee microbiome, strain variation impacts important traits, from pathways that detoxify sugars in the bee’s diet to the T6SS, which underpins interbacterial competition [ 19 , 20 ]. These strains often carry many different mobile genetic elements (MGEs) [ 21 – 24 ] such as phages and plasmids, which are often shared between bacterial species [ 24 ] and carry cargo genes with metabolic functions predicted to benefit honey bees [ 25 ]. In other colony environments, such as the shared food reserves or the guts of larvae or queens[ 26 ], it is less understood what MGEs are present and how they impact symbiont function. Here, we survey MGEs found in B. apis , a symbiont that colonizes multiple colony environments, from larvae and queens to the colony’s food reserves. Since B. apis is exposed to different microbial communities in these environments, we hypothesized that strains from different environments would also carry different MGEs. In particular, we expect that genetic exchange between honey bee symbionts and environmental microbes is most likely to occur at the interface between host digestive tracts and foraged resources: the food reserves ( Figure 1B ). Nectar and pollen are collected from nearby flowers by foragers, stored in cells at the periphery of the hive ( Figure 1A ), and fed to the rest of the colony. At the center of the colony, where the queen lays her eggs, nurse bees feed ( Figure 1C ) the nectar to the queen and larvae by regurgitating from their foregut, also called the crop [ 27 ]. In addition to nectar and pollen, queens and young larvae are fed royal jelly, a secretion of transformed pollen from the nurses’ hypopharyngeal gland ( Figure 1D & E ) [ 28 ]. This flow of resources from the food reserves at the periphery of the hive, to nurses and subsequently the larvae and queen is thought to reduce exposure of larvae and queens to environmental microbes. As nectar and pollen flow to the center of the colony, they are transformed in ways that likely exclude environmental microbes. In the food reserves, nectar stored long-term is dehydrated into honey, and pollen is fermented into bee bread [ 27 , 29 ]. These changes in osmotic pressure and acidity create environments inhospitable to many microbes. Any microbe capable of surviving in the food reserves may be picked up by a nurse bee and fed to the queen or larvae, along with a serving of royal jelly, which is highly toxic to most bacteria due to antimicrobial peptides, low pH, and viscosity [ 30 , 31 ]. This gauntlet of transmission to increasingly inhospitable colony environments means that most transient, environmental microbes do not make it past the food reserves [ 32 – 34 ]. Therefore, it is likely that honey bee symbionts will have the best chance to acquire MGEs when in the food stores, where foragers deposit floral resources with foreign microbes. Among the microbes that are capable of surviving in nectar, and the fermented food that honey bees and many other pollinators eat, are acetic acid bacteria (AAB) and lactic acid bacteria (LAB) [ 35 , 36 ]. Download figure Open in new tab Figure 1: Distribution of Bombella apis across colony environments: (A) Foragers visit flowers and transport nectar (via crop) and pollen (via pollen baskets) back to the colony. The influx of nectar and pollen (B) brings in exogenous microbes: nectar and pollen associates, as well as microbes deposited by other visiting pollinators. These microbes transiently persist in the nectar and pollen, which are incorporated into the diet fed to larvae (D) by nurse bees (C) and to (E) queens by her attendants. Microbial diversity drops in the diet, likely due to the addition of antimicrobial royal jelly (synthesized by nurses and attendants). This combination of factors, described above, that eliminate transient microbes, may explain the low species diversity of the nectar, larval and queen microbiomes. All three microbial communities are dominated by AAB and LAB, including Bombella apis and Apilactobacillus kunkeeii [ 33 , 34 , 37 ] in the nectar and larvae, and B. apis and a consortium of LAB in queens [ 26 , 38 ]. Both B. apis and A. kunkeeii have been shown to have defensive properties against pathogenic fungi and bacteria respectively [ 17 , 39 ]. B. apis is also a nutritional symbiont associated with both queens and larvae, and survives in the presence of high concentrations of royal jelly, which are toxic to A. kunkeeii and other microbiome members [ 18 ]. B. apis is unique in its distribution within the colony, as it is predominantly found in the queen digestive tract, larval gut, nurse hypopharyngeal glands, nurse crops, and foraged resources. B. apis therefore exists in niches that have either high or low exposure to environmental microbes. In this study we isolated new Bombella strains across colony environments – from nectar to crops to queens and larvae – with the aim to understand how MGEs have shaped their genetic diversity. We sequenced the isolates’ genomes, generated complete and closed assemblies, and used those genomes to ask the following questions about plasmids and phages we identified: (1) What MGEs are shared between B. apis strains? (2) What cargo genes do they encode that might shape B. apis ’ phenotypic diversity? (3) Lastly, which bacterial species, both inside and outside the colony, share similar MGEs? Ultimately, we identified many phage clusters unique to our apiary, and one phage cluster conserved across most B. apis genomes. Plasmids were rare, found in only five genomes and not shared between B. apis strains. Despite sharing low genetic similarity, the plasmids all carried toxin-antitoxin systems, and two plasmids (pJPL21 and pAML2) encode phages with high similarity to those found in other B. apis genomes. Phage-plasmids, pJPL21 and pAML2, have CDSs similar to other Bombella and Saccharibacter genomes, whereas plasmids (pJPL24_01, pJPL24_02, and pSME1) had backbones and cargo genes more similar to plasmids from more distantly-related AAB. In short, phages vary colony to colony and have the potential to shape genetic and phenotypic variation between Bombella strains. Plasmids, on the other hand, are rare among B. apis strains and their similarity to other AAB isolated from plants, other insects, and fermented foods, suggests they may have been acquired from AAB from outside the colony. Materials and methods Strain collection and identification Late instar larvae, nurse crops, and queen guts from two bee colonies in the Indiana University research apiary were suspended in 1X PBS, plated on MRS media and incubated for 48 hours at 34° C. Isolated bacteria were cryoarchived, extracted with Qiagen’s MagAttract HMW DNA kit, PCR amplified (using 27F, (AGAGTTTGATCMTGGCTCAG) and 1492R (TACGGYTACCTTGTTACGACTT) primers) and Sanger sequenced at QuintaraBio. Top NCBI BLASTn [ 40 ] hits (>99% sequence identity) was used to identify species. gDNA extraction and long-read sequencing Nineteen new isolates and one previously identified isolate were chosen for long-read, whole genome sequencing. gDNA extractions followed Quiagen’s MagAttract HMW DNA’s gram-negative protocol. Quality control of extractions was determined with a Qubit fluorometer and Nanodrop spectrophotometer. gDNA was sheared to approximately 13kb fragments with a Megaruptor 3 sonicator and libraries were prepared using the SMRTBell Express Template Prep Kit 2.0. One SMRTcell 8M was used to sequence libraries on the PacBio Sequel lle on CCS sequencing mode and a 30hs movie time. Circular Consensus Sequence (CCS) analysis was performed on the instrument using SMRTLink V11 (parameters: ccs –min passes 3 –min rq 0.99) to obtain HiFi reads. Sequencing QC, assembly, and annotation HiFi reads were checked for quality using longQC (parameter: -x pb-sequence) [ 41 ]. Any samples with more than one GC peak were removed. Reads were assembled with Canu (parameters: genomeSize=2M -pacbio-hifi) [ 42 ] three times to produce three distinct assemblies (since Canu is non-deterministic) which were used to generate a consensus assembly with Trycycler [ 43 ]. Default parameters were used unless otherwise stated. Contig clusters were removed if found in only one assembly. During reconciling, clusters with mismatched lengths were excluded. Consensus assemblies were used for all downstream analysis. Annotation was done with NCBI’s PGAP software version 2022-08-11.build 627; this annotation was used in later analysis unless otherwise stated [ 44 ]. For already published genomes, the NCBI entries annotated with PGAP were also used so that annotations would be consistent across genomes for ortholog analysis. This reduces false positive gains or losses due to different gene calling software. Core ortholog phylogeny Proteins were clustered into ortholog groups by requiring all proteins in a cluster to be reciprocal best BLASTP [ 40 ] hits with all other proteins in the cluster (i.e. complete linkage clustering of reciprocal best hits). In addition, if two or more proteins within a genome were more similar to each other than they were to any protein in any other genome (likely duplications within the lineage leading to that genome), all such proteins were placed into the same ortholog group. Protein alignments of ortholog group members were calculated with LINSI [ 45 ]. The corresponding codon alignments were derived from these protein alignments and concatenated into a 1.8 MB alignment. Before concatenation, a filter was used to exclude some ortholog groups based on the following criteria: 1) The group had to include sequences from 20 or more of the genomes, and 2) The median all vs all pairwise amino acid identity had to be 60% or higher. A pairwise amino acid identity was defined as (number of identical amino acids) / max(seq1_len, seq2_len) * 100. This last criterion was to reduce misalignments of divergent orthologs, groups with large disparities in group member length, and ortholog groups containing non orthologous genes. All available nucleotide evolution models were tested using IQTREE [ 46 ]. The best model found based on the Bayesian information criterion was the GTR model with invariant sites and a gamma rate distribution using 30 categories. The final tree was calculated with this model along with 1000 bootstrap replicates in IQTREE. Gene gain/loss analysis Gene gains and losses were inferred for each ortholog group relative to the species tree based on maximum parsimony with the cost of a gain set to 2.001 times the cost of a loss. The gainLoss program used by the GLOOME web server [ 47 ] was run locally for this analysis. COG category assignment A representative sequence (based on median length) from each ortholog group was BLAST-ed [ 40 ] against COG member proteins (3,213,025 sequences) from the 2020 COG database [ 48 ]. Only COG member proteins in “COG membership classes” 0 and 1 were considered, as these member proteins cover most of the COG. For a Bombella ortholog to be assigned to a COG category it had to overlap ≥ 60% of the COG members range that aligned to its COG, meaning that in most cases the Bombella protein would align with ≥ 60% of the COG itself, give the “COG memberships class” constraint above. The highest BLASTP bit scores were considered first. If lower scoring matches to different COGs existed, the Bombella protein could be assigned to multiple COG categories if the above constraints were met and the Bombella region matching such a COG did not overlap with regions already assigned to a COG by more than 60%. Of the 5924 Bombella orthologs that were assigned to some COG only 17 were assigned to 2 COGs and none to more than 2. Curation of plasmid annotation Plasmid annotation was manually refined as follows: the sequence similarity of all loci and intergenic sequences (IGSs) were searched in NCBI’s non-redundant protein sequences (nr) database through Blastx using default parameters [ 40 ]. Automatic annotation artifacts (e.g. split loci, unannotated loci) were corrected when detected. Endogenous replication initiator proteins (i.e. either repA, repB, or repC, depending on the plasmid) were set as the first nucleotide in the assembly. Final GenBank files were fed to the Proksee web server [ 49 ] in order to generate plasmids’ graphical maps. Categorization of coding sequences by COG functional category was completed with eggNOG [ 50 ]. Plasmids visualization in electrophoresis All Bombella strains were grown overnight in 4 mL MRS (34° C, 250 rpm). Liquid cultures were split downstream for either total DNA extraction or plasmid DNA extraction with Qiagen kits (DNeasy or MiniPrep). Approximately 150 ng DNA resulting from each extraction protocol were analyzed in a 0.5% agarose gel electrophoresis (80 Volts, 100 min). Plasmid and phage coding sequence similarity analysis Plasmid and phage coding sequences were compared to all publicly available genomes on NCBI using cblaster [ 51 ], setting a minimum query coverage and percent identity to 50%. To increase the stringency of our search we further modified the parameters. Namely we required that at least 25% of plasmid coding sequences were present in a neighborhood less than or equal to 1.5 * length of the query plasmid or phage. This threshold was determined for the following reasons: (1) plasmid sequences misassembled into bacterial genomes will be linear sequences, and therefore genes on opposite ends of the contig may actually be neighbors; (2) plasmids and phages expand and contract in length frequently as they acquire or lose new DNA. By allowing for a 50% increase in size, we can capture similar plasmids or phages that have expanded in length. The amino acid similarity and synteny of plasmids with similarity to B. apis plasmids were visualized with gene arrow maps using clinker [ 52 ]. Plasmid nucleotide sequence similarity analysis Sourmash [ 53 ] was used to calculate pairwise Jaccard similarity between query plasmids and curated databases of plasmid sequences. Plasmid databases were compiled from one of three sources: (1) PLSDB sequences isolated from honey bees, (2) the top 100 tBLASTn [ 40 ] hits to each plasmid’s replication initiation protein, and (3) cblaster hits from verified plasmid sequences. In sourmash, k-mers of 31 bp were used, as recommended for the genetic relatedness between our input sequences. Identification, annotation, and validation of prophages All B. apis genomes, including genomes recovered from NCBI, were screened for the presence of prophages with VIBRANT [ 54 ]. Recovered prophages were annotated using DRAMv [ 54 ] Annotations were then queried for the presence of canonical genes for prophage excision [ 55 ] Prophages without any canonical gene for prophage excision were classified as degenerate while prophages with canonical excision genes and at least one structural gene were classified as complete. Only prophages with sufficient genome quality (Medium or higher) were considered.[ 56 ] Recovered prophages were then clustered into genus and species-level OTUs with dRep [ 57 ] dRep dereplicate -l 2000 --ignoreGenomeQuality -pa 0.8 -sa 0.95 -nc 0.85 -comW 0 -conW 0 -strW 0 -N50W 0 -sizeW 1 -centW 0). These parameters are informed by observations of known phage biological and ecological diversity [ 58 ]. Prophage-encoded auxiliary metabolic genes (AMGs) were inspected by comparing the output of VIBRANT and DRAM-v. Results Closely related strains co-occur within the same colony To compare the genotypic diversity of B. apis strains between bee colonies and within colony environments, bacteria were isolated from either nectar ( Fig. 1B ), nurse crops ( Fig. 1C ), larvae ( Fig. 1D ), or queen guts ( Fig. 1E ), between two bee colonies and identified as B. apis by 16S sequencing. The resulting nineteen new isolates ( Table 1 ) - as well as SME1, a nectar isolate collected four years prior for which we only had short-reads [ 59 ] - were long-read sequenced on PacBio’s Sequel II platform with an average of ∼300x coverage and assembled with a combination of Canu and Trycyler [ 42 , 43 ] ( Table 1 ). To determine the phylogenetic placement of our new genomes, an ortholog tree was inferred ( Fig. 2A-B ), using all publicly available Bombella genomes and other closely related insect and flower-associated AAB as outgroups. Of the nineteen new Bombella strains, two larval strains were placed within the Bombella favorum clade, a species previously identified in honey bee nectar provisions [ 60 ] ( Fig. 2A ). B. apis genomes were GC rich (∼59%), as expected, with lengths around 2 Mb ( Table 1 ). Genome length across the dataset varied by up to 100,000 bp, ranging from 1.98 to 2.2 Mb in total, hinting at the gain or loss of large regions. Download figure Open in new tab Figure 2: Bombella strains isolated from the same colony and/or environment are phylogenetically distant: A. Core ortholog phylogeny of all sequenced Bombella and (B) Bombella apis genomes and closely-related acetic acid bacteria from insects, plants, and fermented foods. Phylogeny was inferred from 2,711 orthologs (see Methods) using IQ-tree. New genomes generated from this study are highlighted in black. Isolation source is denoted by colored boxes, and, when applicable, the host species and colony environment are labeled with acronyms. View this table: View inline View popup Download powerpoint Table 1: Genome ID, metrics, and isolation environments for newly-sequenced genomes Some strains isolated from the same colony (colony 1: JPL; colony 2: DLM) were closely related and formed monophyletic clades (i.e. DLM19-18) ( Fig. 2B ). Other strains were more closely related to isolates from Europe (i.e. ESL0387 from Switzerland and JPL5; TMW2.1891 from Germany and DLM17). In short, strains sampled from our apiary were diverse, and spanned the B. apis phylogeny. It is worth noting that B. apis strains are all highly related - regardless of isolation source. The average nucleotide identity (ANI) across all strains is ∼99%. Recent gene gains and losses are driven by mobile genetic elements and defense islands To better understand the diversification of strains within the same colony and/or apiary, we identified genes gained or lost within the B. apis clade using GLOOME [ 47 ]. Specifically, we were interested in the most recent gene gains - genes that are unique to each new strain we sequenced. To better visualize the broad pattern of what functions are newly gained and lost, we classified gained/lost genes by their COG (Cluster of Orthologous Groups of proteins) category (24) ( Fig. 3 ; Fig. S1). Across the B. apis tree, 24.3% of gene gains were classified by COG as mobile elements. Gene gains unique to our strains were largely involved in replication and repair (L), defense mechanisms (V), and mobile elements (X) ( Fig.3A ). Gene gains due to mobile elements were largely driven by phages ( Fig. 3C , suggesting that phages are contributing to genetic diversity among strains of the same colony. Download figure Open in new tab Figure 3: Strain-specific gains are biased towards mobile elements and defense islands: (A) A. OG categories of genes gained within the B. apis clade. Gene gains and losses were defined with LOOME and categorized by their COG family. Pie charts denote the COG categories of genes ained at that node or leaf. Pie chart colors represent COG categories (U: intracellular trafficking, ecretion and vesicular transport; O: posttranslational modification, protein turnover and chaperones; Q: secondary metabolite biosynthesis, transport, and catabolism; G: carbohydrate transport and etabolism; E: amino acid transport and metabolism; F: nucleotide transport and metabolism; H: oenzyme transport and metabolism; I: lipid transport and metabolism; P: inorganic ion transport and etabolism; J: translation, ribosomal structure and biogenesis; D: cell cycle control, cell division, and hromosome partitioning; K: transcription; L: replication, recombination, and repair; V: defense echanisms; T: signal transduction mechanisms; M: cell wall/membrane/envelope biogenesis; N: cell otility; W: extracellular structures; X: mobilome, prophage, and transposons; C: energy production nd conversion) and the diameter of each pie represents the number of genes gained. (B) Plasmid resence. Presence of at least one plasmid identified during genome assembly is denoted in filled quares. (C) Presence of phage within each genome. Phage presence and absence was annotated ith Vibrant. Phages with incomplete integrase and structural capsid genes are considered egenerate. VIRDIC was used to bin phage sequences at the genus level, identifying 14 distinct enera present in this clade. Identification of phages within B. apis genomes sampled from different colonies To identify groups of related phages – termed clusters – shared between Bombella genomes from both colonies, we used VIBRANT[ 56 ] for prophage prediction and dRep [ 57 ] to group phage sequences into related clusters. With this method, we identified 91 prophages of sufficient quality, which were then clustered into seven genus-level clusters (80% ANI), five of which are found in B. apis genomes isolated from our colonies ( Fig. 4C ). Sixty percent of phage clusters are found in both colonies (Fig. S4A, C), with two clusters (2 and 4) unique to colonies D and J respectively. To examine phage distribution at a finer genetic resolution, we grouped phages with at least 95% ANI into subclusters (34 distinct subclusters). Very few subclusters (23.5%) were identified in both colonies (Fig. S4B). Rather, most phage subclusters were unique to a single colony (76.5%), and in many cases, unique to one B. apis genome (Fig. S4D). Together, this suggests that although phage clusters are often shared between colonies, rapid diversification of phages may lead to distinct subclusters, unique to colonies. Download figure Open in new tab Figure 4: B. apis phages have high genic similarity with CDSs from other Bombella species and Saccharibacter genomes. (A) The proportion of phage-encoded genes present in bacterial genomes. One representative phage from each cluster was used as a tBLASTn query against all publicly available genomes on NCBI, requiring 50% query coverage and percent amino acid identity. Only three representative phages returned hits. To determine if CDSs were syntenic -as would be expected for phage-encoded genes -cblaster was used, requiring at least 25% of phage encoded genes to be present in a neighborhood less than or equal to the phage length* 1.5 to allow for phage size expansion. If all the above criteria were met, the bacterial genome was considered a hit. The proportion of query phage proteins present in each hit genome is shown in the first panel (A). For each phage cluster that returned hits, the hits were classified by the (B) genus and (C) species of its host genome assembly. Assemblies with no designation at the genus level are labeled ‘na’ and those without designated species are labeled ‘sp’. Similarity of B. apis phages to CDSs from other Bombella and Saccharibacter strains To determine if B. apis may be acquiring phages from other bacterial species, either other honey bee symbionts or environmental bacteria, we searched for any bacterial genomes with high amino acid similarity to phages identified in our colonies. Querying all publicly available genomes on NCBI using cblaster [ 51 ], we identified gene clusters with at least 50% amino acid identity and 50% query coverage to representative phages from each phage cluster. Three of seven phage clusters returned hits ( Fig. 4A ), suggesting that either phages or phage-encoded cargo genes are shared with these genomes. All hits were to AAB, specifically Bombella and Saccharibacter , a clade of flower- and solitary bee-associated AAB sister to Bombella [ 61 – 63 ] ( Fig 4B ), and most hits were from B. apis or other Bombella species isolated from honey or bumble bees ( Fig. 4C ). What functions the shared CDSs encode is unclear, given that most sequences shared between phages and cblaster hits are of unknown function (Fig. S3). Identification of plasmids from nectar, nurse crop, and queen gut isolates During genome assembly, we identified five closed extrachromosomal elements in four genome assemblies (AML2, JPL21, JPL24, and SME1) ( Fig. 5B ). To determine if these closed contigs were plasmids we manually annotated with BlastX to identify replication machinery. In all five extrachromosomal elements, replication proteins were present ( Fig. 5A ), belonging to either RepA, RepB, or RepC families (Table S1). Additionally, we visualized DNA fragments of approximately the correct length in plasmid preps from each strain (Figure S2). Interestingly, all plasmids were found in isolates from nectar, nurse crops, and queens ( Fig. 5A ), but never in larval isolates ( Figure 5A ). These are the first sequenced plasmids from Bombella , although it is likely that short-read sequencing may be obscuring plasmids in older assemblies. To quantify overall nucleotide similarity between plasmids, we measured Jaccard similarity with sourmash [ 53 ] . Two pairs of plasmids (pJPL24_01 and pSME1; pAML2 and pJPL21) are similar, with a Jaccard Index of 0.1 ( Fig. 5B ). Further annotation of these plasmids revealed that three of five plasmids have conjugation machinery, suggesting a possible HGT mechanism. Download figure Open in new tab Figure 5: Five unique plasmids were identified in isolates from nectar, nurse crops, and queen guts. (A). Four strains harbored plasmids (JPL21, JPL24, AML1, SME1). Plasmid replication machinery and cargo genes are colored by annotation. (B). Overall nucleotide similarity of B. apis plasmids to one another was measured with Jaccard similarity using sourmash. Scores above 0.1 indicate that plasmids are significantly similar. (C) COG categories of plasmid-encoded genes. Plasmid coding sequences were annotated by eggNOG. Colors denote COG functional categories. (D) Alignment of putative phage-plasmids pJPL21 and pAML2 showing percent amino acid similarity of phage and plasmid genes as annotated by Pfam. Plasmids encode toxin-antitoxin systems, antimicrobial resistance cassettes, and prophage regions To gain insight into the function of the plasmid cargo genes, we annotated them with eggNOG (Table S2) in addition to PGAP. Most proteins were not classified by COG categories ( Fig.4C ). Of these unclassified proteins in pAML2 and pJPL21, many had high homology to phage proteins in Pfam and annotation with VIBRANT confirmed the presence of similar phages genes on both pAML2 and pJPL21. Since both encode plasmid replication proteins in addition to this phage, pAML2 and pJPL21 would more appropriately be classified as phage-plasmids (P-Ps). When compared to phages identified in B. apis genomes, both P-Ps were placed in cluster 1, a cluster not otherwise represented in our colonies ( Fig. 3C ). In addition to sharing similar phages, both P-Ps have high genetic similarity (Jaccard index: 0.09) (Fig. S6B); the larger of the two, pJPL21, encodes 86% of genes found on the smaller P-P, pAML2, but has a stretch of approximately 17,000 bps unique to it, including a tetracycline resistance cassette (Fig.S6A). In short, both P-Ps are similar to each other, but otherwise bear no similarity to other phages or plasmids found within our colonies, instead clustering with phages found in Bombella isolates from Europe. One unifying functional feature across all plasmids was the presence of toxin-antitoxin (T-A) systems ( Figure 5A ; Table S1). All T-A systems identified are Type II, producing a stable toxin protein and an unstable antitoxin. Likely these systems are involved in plasmid maintenance but could also play a role in phage defense [ 64 ]. Plasmid coding sequences share similarity with acetic acid bacteria (AAB), including other Bombella apis strains Since plasmids are mosaic, and frequently take up and lose their cargo genes, taxonomically classifying plasmids based on sequence similarity is challenging. To approach this classification issue, we took two different approaches: (1) one to find syntenic regions with similar coding sequences to our plasmids and (2) to identify plasmids with similar backbones. The first approach repurposed software developed for biosynthetic gene clusters to identify contiguous clusters of coding sequences (CDS) with high similarity to CDSs from our plasmids. This approach has the advantage of identifying clusters of similar CDSs regardless of whether they are located on chromosomes, plasmids, etc. The second approach used nucleotide similarity across short regions (k=31 bp) to identify closely related sequences. This approach is more stringent and has a higher chance of recovering plasmid sequences with the same backbone, or non-coding sequences. Using the first approach, we identified similar gene clusters in AAB and in one LAB from honey bees ( Figure 6A ). Acetobacter , Gluconobacter , and Komagaetibacter CDSs cluster with three of five plasmids. On average, around half the CDSs encoded on each plasmid were found in other genomes. The B. apis plasmids pJPL24_01, pJPL24_02, and pSME1 all had similar CDSs from known plasmids ( Figure 6A and B ). On the other hand the P-Ps, pAML2 and pJPL21, only shared similar CDSs with other Bombella (formerly Parasaccharibacter ) and Saccharibacter genomes ( Figure 6A and B ). Two Bombella genome assemblies (LMG1352 and TMW 2.1884) shared 91% and 88% of pAML2’s CDSs respectively. Saccharibacter assemblies (EH611, EH60 and EH70) by comparison shared only 55-58% of pAML2’s CDSs ( Figure 6A ). All hits to our plasmids were either AAB or LAB and isolated from acidic environments such as insect guts, nectar, or fermented foods ( Figure 6C ). P-Ps, pAML2 and pJPL2, were most similar to other bacterial genomes isolated from honey bees, whereas all other plasmids were most similar to AAB isolated from fermented food, fruits, flowers, and plants ( Figure 6C ). Download figure Open in new tab Figure 6: Plasmid-encoded genes have high similarity to CDS found among other acetic acid bacteria isolated from insects, flowers, and fermented foods. (A) The proportion of plasmid-encoded genes present in bacterial genome assemblies. The presence of CDSs with high similarity to B. apis plasmid CDSs in all NCBI genomes was determined using tBLASTn with 50% query coverage and 50% amino acid identity. To determine if CDSs were present in clusters instead of dispersed across the genome, cblaster was used, requiring that at least 25% of plasmid CDSs were present within a range equal to the (plasmid length) * 1.5 to allow for plasmid expansion. If all the above criteria were met, the assembly was considered a hit. (B) Plasmid CDS hits were characterized by genera or (C) isolation source as reported on NCBI. Plasmid replication and conjugation machinery is highly conserved among a small number of AAB plasmids Our second approach to identify similar plasmids quantified overall nucleotide sequence similarity, which encompasses both the plasmid backbone and coding regions. Using sourmash, plasmids were subset into k-mers of 31 bps, and pairwise nucleotide similarity, as defined by Jaccard Similarity, was calculated. We defined similar plasmid backbones as having a Jaccard Index ≥ 0.1. This threshold has previously been used to classify similar plasmids into “clusters” or “cliques” [ 65 , 66 ]. Using these criteria we queried three different databases of plasmid sequences. First, we queried all the known plasmids that have been isolated from honey bee symbionts and pathogens and are available in the Plasmid Database (PLSDB) but found no hits (Figure S5). Secondly, we created a database of plasmid sequences identified as hits during our CDS cluster similarity analysis ( Figure 6A ). Five plasmids from this database had a Jaccard index value greater or equal to 0.1, indicating similarity ( Figure 7A ). Thirdly, we built a database of plasmid sequences with similarity to our plasmid’s replication proteins (top 100 hits from a tBlastN query). Sequence similarity analysis of this database yielded three new hits with significant Jaccard index values and one hit already identified by CDS cluster similarity ( Komagataeibacter nataicola strain DS12 plasmid pKNA10) ( Figure 7A ). Replication machinery was highly conserved between pJPL24_02 and its top hit, pKHC110_1 ( Figure 7B ). pSME1 followed a similar trend, sharing strong homology with the RepA protein and recombinase from its top hit, 2P ( Figure 7C ). In the case of pJPL24_01, however, a large region (∼6 kb) is highly conserved with pKNA10, a plasmid from Komagataeibacter nataicola . This region encodes RepA, conjugal transfer proteins, TraD and TraA, as well as a VapB antitoxin ( Figure 7D ). Together this suggests that three B. apis plasmids isolated from our apiary possess backbones, replication machinery, and some cargo genes homologous to known plasmids in AAB from fermented foods, but are not similar to any known plasmids from honey bee-associated bacteria. Download figure Open in new tab Figure 7: B. apis plasmids have low overall similarity to plasmids from other AAB but have highly conserved replication and conjugation regions compared to AAB plasmids (A) Similarity of plasmids to other AAB plasmids as determined by Jaccard-Index. Alignment of conserved plasmid regions from (B) pJPL24_02, (C) pSME1, and (D) pJPL24_01. Colors of genes represent PGAP functional annotation. Unaligned genes are shown in gray. Pairwise amino acid identity between CDS is shown in greyscale. Discussion The Acetobacteraceae B. apis is an important symbiont of honey bees – both buffering against poor nutrition and reducing the likelihood of fungal infections during larval development [ 17 , 18 ]. The diversity of B. apis strains within and between colonies and their distribution across colony environments has not previously been examined, since most of our understanding o f B. apis ’ ecology is based on 16S rRNA amplicon studies that lack the resolution to study strain variation [ 33 , 34 , 38 ]. Here we isolated and sequenced nineteen new B. apis strains and identified MGEs, plasmids, phages, and phage-plasmids (P-Ps) that contribute to their genetic diversification. Additionally, we used two complementary approaches to identify MGE relatives in other bacteria, with the goal of identifying possible sources from which B. apis acquired MGEs. In this process we identified and isolated two strains of Bombella favorum , a close-relative of the bumblebee symbiont, Bombella intestini . This finding, in corroboration with prior studies, suggests that other Bombella species are present in the colony [ 60 ]. Whether this association is transient or perpetual is unclear. Given the high similarity of the 16S rRNA operon between Bombella species, it is possible that OTUs ascribed to B. apis in past studies actually represent a consortium of Bombella species, but further studies are needed to determine their association with honey bees. Within the B. apis clade there is no striking correlation between geography and genetic relatedness. Genetic variation in other symbiont populations, such as those associated with bobtail squid and deep sea snails, is largely explained by geographic location [ 67 , 68 ]. This striking lack of correlation is perhaps representative of honey bee management practices. The European honey bee is, after all, a domesticated insect, having been managed by humans for 4,000 years and transported across continents for its pollination services [ 69 , 70 ]. New colonies, which assemble naturally by swarming, are now generally founded by beekeepers with queens bred in queen-breeding facilities [ 71 ]. In the US, only a handful of queen-breeding companies supply most commercial and private beekeepers. Given B. apis’ association with the queen gut, it seems likely that queen-breeding facilities serve as reservoirs, and that colonies across the US are seeded with strains derived from these queens. Strains derived from our colonies did not cluster by their isolation environment: larvae, queens, nurse crops, or nectar. Instead, larval isolates are clustered with strains from all other hive environments. The bias in our sampling towards larval strains impedes our ability to statistically test for correlation between isolation environment and phylogenetic relatedness and draw a more meaningful conclusion. To test for spatial segregation or niche partitioning of strains (as observed in other symbionts[ 72 – 75 ] within the hive, both deeper sampling and longitudinal data are needed. The clearest distinction between larval and non-larval strains was the distribution of plasmids among them. Despite our larger sampling of larval strains, none carried any detectable plasmids, but all other strains had at least one plasmid. Based on this observation, and the fact that geography can dictate plasmid transmission and evolution [ 76 ], it may be worthwhile to examine whether B apis ’ spatial distribution within the colony correlates with plasmid abundance and distribution in future studies. Plasmid transmission and evolutionary history in natural microbial communities are still poorly understood. Due to the mosaic composition of plasmids, and other mobile genetic elements, we need different evolutionary frameworks to disentangle their patterns of diversification [ 24 , 76 ]. To search publicly available genomes for similar plasmids or plasmid cargo genes, we adopted a tool, cblaster, developed for other horizontally transferred genes: biosynthetic gene clusters [ 51 ]. This method allowed us to survey draft genome assemblies in an unbiased fashion. Some of our hits were also confirmed by our alternative method (Jaccard similarity) when used on a subset of known plasmid sequences. Whether the CDS present in other AAB assemblies are chromosomal or on plasmids is not known but would be an interesting avenue of future study. The coding sequences present in other AAB are largely unannotated ( Figure 4C ), but their prevalence across diverse AAB may be indicative of a function important for this phylogenetic group. Alternatively, if these CDSs are localized to plasmids, their distribution throughout AAB may be the result of the plasmid’s broad host range and horizontal transfer. How did Bombella acquire these plasmids? The strong similarity in both synteny and gene content between the two P-Ps and other B. apis and S. floricola assemblies suggests that highly similar P-Ps may be stably associated with honey bees. The other three plasmids, however, have low similarity to plasmids from AAB isolated from plants and fermented foods. Typically, CDSs shared between B. apis and AAB plasmids were involved in replication, recombination and conjugation and did not include cargo genes carried on the plasmid backbone. Perhaps the similarities between them are indicative of a common plasmid backbone with a broad host range, that has infected AAB sporadically throughout their evolutionary history. Likely, B. apis acquired these three plasmids from a plant-associated AAB, but whether vertically or horizontally is challenging to disambiguate. Both insect- and fermented food associated-AAB diverged from a plant-associated clade [ 36 ]. Therefore the similarity of Bombella ’s plasmids to genera associated with flowers, namely Neokomagataea [ 77 , 78 ], Gluconobacter , and Saccharibacter [ 62 ] could be indicative of an ancient association or a recent acquisition. In the latter scenario, flowers may act as transmission hubs for the exchange of AAB symbionts and their MGEs [ 35 ]. This phenomenon has been documented for pollinator pathogens, which can swap hosts via shared floral resources [ 79 ]. Regardless of how B. apis acquired these plasmids, it is clear that they belong to a group of similar plasmids that can be harbored by AAB across diverse environments. Most studies of MGEs in honey bees focus on the gut microbiome of adult workers -which is compositionally distinct from the communities in which B. apis resides. Connecting our findings to these studies of phages and plasmids in the adult microbiome, however, reveals some commonalities in the distribution of MGEs within and across colonies. In the adult gut, plasmids and phage clusters are shared between different strains and species within the same colony, as well as across colonies in the same apiary [ 24 ]. For phages, the same clusters can be identified in bacterial genomes from colonies worldwide, suggesting a widespread association. It is likely that some of these phages are also P-Ps, given that many also encode RepC proteins. In the adult worker microbiome, phages often encode metabolic genes with obvious links to host nutrition, but the same is not true for B. apis phages [ 25 ]. Future studies interrogating the function of the phage cargo genes could elucidate their role in shaping B. apis ’ phenotypic diversity. Overall, this study has identified novel phages, plasmids and phage-plasmids associated with the important honey bee symbiont B. apis and documented their distribution across colonies and colony environments. Our work highlights the role of phages in driving genetic diversity between B. apis strains within the same colony and apiary - a finding reflected in studies of other honey-bee microbiome members. In contrast, the plasmids and P-Ps we identified were rare within our sampling and only found in a few strains isolated from queens, nurses, and food stores. While both P-Ps clustered with other B. apis phages, the plasmids had low genic similarity to AAB plasmids from plants and fermented foods. To better interpret these findings, B. apis ’ transmission within and between colonies should be explored, as well as the influx of AAB into the colony. The interconnectedness of AAB populations, linked by the trophic interactions of their hosts, represents a fascinating area of study for tracing the transmission of mobilized traits across different symbiont populations. Funder Information Declared NSF , NSF DBI 2022049 USDA , 2022-67011-36577 References 1. ↵ McFall-Ngai M et al. Animals in a bacterial world, a new imperative for the life sciences . Proceedings of the National Academy of Sciences 2013 ; 110 : 3229 – 3236 . doi: 10.1073/pnas.1218525110 OpenUrl Abstract / FREE Full Text 2. ↵ Louca S et al. Function and functional redundancy in microbial systems . Nat Ecol Evol 2018 ; 2 : 936 – 943 . OpenUrl 3. ↵ Dewar AE et al. Bacterial lifestyle shapes pangenomes . Proceedings of the National Academy of Sciences 2024 ; 121 : e2320170121 . doi: 10.1073/pnas.2320170121 OpenUrl CrossRef 4. ↵ Flórez LV et al. An antifungal polyketide associated with horizontally acquired genes supports symbiont-mediated defense in Lagria villosa beetles . Nat Commun 2018 ; 9 : 2478 . OpenUrl CrossRef PubMed 5. Heath KD et al. MGEs as the MVPs of Partner Quality Variation in Legume-Rhizobium Symbiosis . mBio 2022 ; 13 : e00888 – 22 . doi: 10.1128/mbio.00888-22 OpenUrl CrossRef PubMed 6. ↵ Hillman K , Goodrich-Blair H . Are you my symbiont? Microbial polymorphic toxins and antimicrobial compounds as honest signals of beneficial symbiotic defensive traits . Curr Opin Microbiol 2016 ; 31 : 184 – 190 . doi: 10.1016/j.mib.2016.04.010 OpenUrl CrossRef 7. ↵ Carr VR et al. Probing the Mobilome: Discoveries in the Dynamic Microbiome . Trends in Microbiology 2021 ; 29 : 158 – 170 . doi: 10.1016/j.tim.2020.05.003 OpenUrl CrossRef 8. ↵ Hooykaas PJJ , Snijdewint FGM , Schilperoort RA . Identification of the sym plasmid of Rhizobium leguminosarum strain 1001 and its transfer to and expression in other Rhizobia and Agrobacterium tumefaciens . Plasmid 1982 ; 8 : 73 – 82 . doi: 10.1016/0147-619X(82)90042-7 OpenUrl CrossRef PubMed Web of Science 9. Weisberg AJ et al. Pangenome Evolution Reconciles Robustness and Instability of Rhizobial Symbiosis . mBio 2022 ; 13 : e00074 – 22 . doi: 10.1128/mbio.00074-22 OpenUrl CrossRef PubMed 10. ↵ Wardell GE et al. Why are rhizobial symbiosis genes mobile? Philos Trans R Soc Lond B Biol Sci 2022 ; 377 : 20200471 . OpenUrl CrossRef PubMed 11. ↵ Van Arnam EB , Currie CR , Clardy J . Defense contracts: molecular protection in insect-microbe symbioses . Chem Soc Rev 2018 ; 47 : 1638 – 1651 . doi: 10.1039/C7CS00340D OpenUrl CrossRef 12. ↵ LePage DP et al. Prophage WO Genes Recapitulate and Enhance Wolbachia-induced Cytoplasmic Incompatibility . Nature 2017 ; 543 : 243 – 247 . doi: 10.1038/nature21391 OpenUrl CrossRef PubMed 13. ↵ Johnson J , Soehnlen M , Blankenship HM . Long read genome assemblers struggle with small plasmids . Microb Genom 2023 ; 9 : mgen001024 . doi: 10.1099/mgen.0.001024 OpenUrl CrossRef 14. Wick RR , Holt KE . Benchmarking of long-read assemblers for prokaryote whole genome sequencing . F1000Res 2019 ; 8 : 2138 . OpenUrl 15. ↵ Ho SFS et al. Gauge your phage: benchmarking of bacteriophage identification tools in metagenomic sequencing data . Microbiome 2023 ; 11 : 84 . doi: 10.1186/s40168-023-01533-x OpenUrl CrossRef PubMed 16. ↵ Motta EVS , Moran NA . The honeybee microbiota and its impact on health and disease . Nat Rev Microbiol 2024 ; 22 : 122 – 137 . doi: 10.1038/s41579-023-00990-3 OpenUrl CrossRef PubMed 17. ↵ Miller DL , Smith EA , Newton ILG . A Bacterial Symbiont Protects Honey Bees from Fungal Disease . mBio 2021 . doi: 10.1128/mbio.00503-21 OpenUrl CrossRef 18. ↵ Parish AJ et al. Honey bee symbiont buffers larvae against nutritional stress and supplements lysine . ISME J 2022 ; 16 : 2160 – 2168 . doi: 10.1038/s41396-022-01268-x OpenUrl CrossRef 19. ↵ Zheng H et al. Metabolism of toxic sugars by strains of the bee gut symbiont Gilliamella apicola . MBio 2016 ; 7 . 20. ↵ Steele MI et al. Diversification of type VI secretion system toxins reveals ancient antagonism among bee gut microbes . mBio 2017 ; 8 . doi: 10.1128/mBio.01630-17 OpenUrl Abstract / FREE Full Text 21. ↵ Honey bees harbor a diverse gut virome engaging in nested strain-level interactions with the microbiota | PNAS . https://www.pnas.org/doi/10.1073/pnas.2000228117 . . 22. Busby TJ et al. Global Composition of the Bacteriophage Community in Honey Bees . mSystems ; 7 : e01195 – 21 . doi: 10.1128/msystems.01195-21 OpenUrl CrossRef 23. Bueren EK et al. Characterization of prophages in bacterial genomes from the honey bee (Apis mellifera) gut microbiome . PeerJ 2023 ; 11 : e15383 . OpenUrl CrossRef PubMed 24. ↵ Robinson CRP et al. Mobile genetic elements exhibit associated patterns of host range variation and sequence diversity within the gut microbiome of the European Honey bee. 2025 . bioRxiv , 2025 . , 2025.01.31.635958 25. ↵ Bueren EK et al. Characterization of prophages in bacterial genomes from the honey bee (Apis mellifera) gut microbiome . PeerJ 2023 ; 11 : e15383 . doi: 10.7717/peerj.15383 OpenUrl CrossRef 26. ↵ Caesar L et al. Metagenomic analysis of the honey bee queen microbiome reveals low bacterial diversity and Caudoviricetes phages . mSystems 2024 ; 9 : e01182 – 23 . doi: 10.1128/msystems.01182-23 OpenUrl CrossRef PubMed 27. ↵ Wright GA , Nicolson SW , Shafir S . Nutritional Physiology and Ecology of Honey Bees . Annual Review of Entomology 2017 ; 63 : 327 – 344 . doi: 10.1146/annurev-ento-020117-043423 OpenUrl CrossRef PubMed 28. ↵ Crailsheim K . The flow of jelly within a honeybee colony . J Comp Physiol B 1992 ; 162 : 681 – 689 . doi: 10.1007/BF00301617 OpenUrl CrossRef Web of Science 29. ↵ Di Cagno R et al. Novel solid-state fermentation of bee-collected pollen emulating the natural fermentation process of bee bread . Food Microbiology 2019 ; 82 : 218 – 230 . doi: 10.1016/J.FM.2019.02.007 OpenUrl CrossRef 30. ↵ Fontana R et al. Jelleines: a family of antimicrobial peptides from the Royal Jelly of honeybees (Apis mellifera) . Peptides 2004 ; 25 : 919 – 928 . doi: 10.1016/j.peptides.2004.03.016 OpenUrl CrossRef PubMed Web of Science 31. ↵ Xue X , Wu L , Wang K. Chemical Composition of Royal Jelly . Bee Products -Chemical and Biological Properties . Cham : Springer International Publishing , 2017 , 181 – 190 . 32. ↵ Yun J-H et al. Social status shapes the bacterial and fungal gut communities of the honey bee . Scientific Reports 2018 ; 8 : 2019 . doi: 10.1038/s41598-018-19860-7 OpenUrl CrossRef 33. ↵ Vojvodic S , Rehan SM , Anderson KE . Microbial Gut Diversity of Africanized and European Honey Bee Larval Instars . PLoS ONE 2013 ; 8 : 72106 . doi: 10.1371/journal.pone.0072106 OpenUrl CrossRef 34. ↵ Anderson KE et al. Microbial Ecology of the Hive and Pollination Landscape: Bacterial Associates from Floral Nectar, the Alimentary Tract and Stored Food of Honey Bees (Apis mellifera) . 2013 . doi: 10.1371/journal.pone.0083125 OpenUrl CrossRef PubMed 35. ↵ Crotti E et al. Acetic acid bacteria, newly emerging symbionts of insects . Applied and Environmental Microbiology 2010 ; 76 : 6963 – 6970 . doi: 10.1128/AEM.01336-10 OpenUrl Abstract / FREE Full Text 36. ↵ Cuesta-Maté A et al. Evolutionary genomics reveals plant origins of acetic acid bacteria in fermented food. 2025 . bioRxiv , 2025 . , 2025.07.09.663944 37. ↵ Kowallik V , Mikheyev AS . Honey Bee Larval and Adult Microbiome Life Stages Are Effectively Decoupled with Vertical Transmission Overcoming Early Life Perturbations . mBio 2021 ; 12 : e02966 – 21 . doi: 10.1128/mBio.02966-21 OpenUrl CrossRef 38. ↵ Tarpy DR , Mattila HR , Newton ILG . Development of the honey bee gut microbiome throughout the queen-rearing process . Applied and environmental microbiology 2015 ; 81 : 3182 – 91 . doi: 10.1128/AEM.00307-15 OpenUrl Abstract / FREE Full Text 39. ↵ Zendo T et al. Kunkecin A, a New Nisin Variant Bacteriocin Produced by the Fructophilic Lactic Acid Bacterium, Apilactobacillus kunkeei FF30-6 Isolated From Honey Bees . Front Microbiol 2020 ; 11 . doi: 10.3389/fmicb.2020.571903 OpenUrl CrossRef 40. ↵ McGinnis S , Madden TL . BLAST: at the core of a powerful and diverse set of sequence analysis tools . Nucleic Acids Res 2004 ; 32 : W20 – W25 . doi: 10.1093/nar/gkh435 OpenUrl CrossRef PubMed Web of Science 41. ↵ Fukasawa Y et al. LongQC: A Quality Control Tool for Third Generation Sequencing Long Read Data . G3 (Bethesda) 2020 ; 10 : 1193 – 1196 . doi: 10.1534/g3.119.400864 OpenUrl Abstract / FREE Full Text 42. ↵ Koren S et al. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation . Genome Res 2017 ; 27 : 722 – 736 . doi: 10.1101/gr.215087.116 OpenUrl Abstract / FREE Full Text 43. ↵ Wick RR et al. Trycycler: consensus long-read assemblies for bacterial genomes . Genome Biology 2021 ; 22 : 266 . doi: 10.1186/s13059-021-02483-z OpenUrl CrossRef PubMed 44. ↵ Li W et al. RefSeq: expanding the Prokaryotic Genome Annotation Pipeline reach with protein family model curation . Nucleic Acids Res 2020 ; 49 : D1020 – D1028 . doi: 10.1093/nar/gkaa1105 OpenUrl CrossRef 45. ↵ Katoh K . MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform . Nucleic Acids Research 2002 ; 30 : 3059 – 3066 . doi: 10.1093/nar/gkf436 OpenUrl CrossRef PubMed Web of Science 46. ↵ Minh BQ et al. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era . Mol Biol Evol 2020 ; 37 : 1530 – 1534 . doi: 10.1093/molbev/msaa015 OpenUrl CrossRef PubMed 47. ↵ Cohen O et al. GLOOME: gain loss mapping engine . Bioinformatics 2010 ; 26 : 2914 – 2915 . OpenUrl CrossRef PubMed Web of Science 48. ↵ Galperin MY et al. Microbial genome analysis: the COG approach . Brief Bioinform 2019 ; 20 : 1063 – 1070 . doi: 10.1093/bib/bbx117 OpenUrl CrossRef PubMed 49. ↵ Grant JR et al. Proksee: in-depth characterization and visualization of bacterial genomes . Nucleic Acids Res 2023 ; 51 : W484 – W492 . doi: 10.1093/nar/gkad326 OpenUrl CrossRef PubMed 50. ↵ Hernández-Plaza A et al. eggNOG 6.0: enabling comparative genomics across 12 535 organisms . Nucleic Acids Res 2023 ; 51 : D389 – D394 . OpenUrl CrossRef PubMed 51. ↵ Gilchrist CLM et al. cblaster: a remote search tool for rapid identification and visualization of homologous gene clusters . Bioinform Adv 2021 ; 1 : vbab016 . OpenUrl 52. ↵ Gilchrist CLM , Chooi Y-H. clinker & clustermap.js: automatic generation of gene cluster comparison figures . Bioinformatics 2021 ; 37 : 2473 – 2475 . doi: 10.1093/bioinformatics/btab007 OpenUrl CrossRef PubMed 53. ↵ Pierce NT et al. Large-scale sequence comparisons with sourmash . 2019. F1000Research , 2019 . 54. ↵ DRAM for distilling microbial metabolism to automate the curation of microbiome function | Nucleic Acids Research | Oxford Academic . https://academic.oup.com/nar/article/48/16/8883/5884738 . . 55. ↵ Mavrich TN , Hatfull GF . Bacteriophage evolution differs by host, lifestyle and genome . Nat Microbiol 2017 ; 2 : 17112 . doi: 10.1038/nmicrobiol.2017.112 OpenUrl CrossRef PubMed 56. ↵ Kieft K , Zhou Z , Anantharaman K . VIBRANT: automated recovery, annotation and curation of microbial viruses, and evaluation of viral community function from genomic sequences . Microbiome 2020 ; 8 : 90 . OpenUrl CrossRef PubMed 57. ↵ Olm MR et al. dRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication . ISME J 2017 ; 11 : 2864 – 2868 . doi: 10.1038/ismej.2017.126 OpenUrl CrossRef PubMed 58. ↵ Gregory AC et al. Marine DNA Viral Macro-and Microdiversity from Pole to Pole . Cell 2019 ; 177 : 1109 – 1123 .e14. doi: 10.1016/j.cell.2019.03.040 OpenUrl CrossRef 59. ↵ Draft Genome Sequence of a Bombella apis Strain Isolated from Honey Bees | Microbiology Resource Announcements . https://journals.asm.org/doi/full/10.1128/mra.01329-19 . . 60. ↵ Hilgarth M et al. Bombella favorum sp. nov. and Bombella mellum sp. nov., two novel species isolated from the honeycombs of Apis mellifera . Int J Syst Evol Microbiol 2021 ; 71 . doi: 10.1099/ijsem.0.004633 OpenUrl CrossRef 61. ↵ Smith EA , Newton ILG . Genomic signatures of honey bee association in an acetic acid symbiont . bioRxiv 2018 ; 367490 . doi: 10.1101/367490 OpenUrl Abstract / FREE Full Text 62. ↵ Jojima Y et al. Saccharibacter floricola gen. nov., sp. nov., a novel osmophilic acetic acid bacterium isolated from pollen . International Journal of Systematic and Evolutionary Microbiology 2004 ; 54 : 2263 – 2267 . doi: 10.1099/ijs.0.02911-0 OpenUrl CrossRef PubMed Web of Science 63. ↵ Smith EA et al. Draft Genome Sequences of Four Saccharibacter sp. Strains Isolated from Native Bees . Microbiology Resource Announcements 2020 ; 9 . doi: 10.1128/mra.00022-20 OpenUrl CrossRef 64. ↵ LeRoux M , Laub MT . Toxin-Antitoxin Systems as Phage Defense Elements . Annu Rev Microbiol 2022 ; 76 : 21 – 43 . doi: 10.1146/annurev-micro-020722-013730 OpenUrl CrossRef PubMed 65. ↵ Vereau Gorbitz D et al. Plasmid transmission dynamics and evolution of partner quality in a natural population of Rhizobium leguminosarum . bioRxiv 2024 ;2024.10.17.618979. 66. ↵ Acman M et al. Large-scale network analysis captures biological features of bacterial plasmids . Nat Commun 2020 ; 11 : 2452 . doi: 10.1038/s41467-020-16282-w OpenUrl CrossRef PubMed 67. ↵ Rotman ER et al. Natural Strain Variation Reveals Diverse Biofilm Regulation in Squid-Colonizing Vibrio fischeri . J Bacteriol 2019 ; 201 . 68. ↵ Hauer MA et al. Geography, not lifestyle, explains the population structure of free-living and host-associated deep-sea hydrothermal vent snail symbionts . Microbiome 2023 ; 11 : 106 . doi: 10.1186/s40168-023-01493-2 OpenUrl CrossRef 69. ↵ Dams M , Dams L . Spanish rock art depicting honey gathering during the Mesolithic . Nature 1977 . doi: 10.1038/268228a0 OpenUrl CrossRef 70. ↵ vanEngelsdorp D , Meixner MD. A historical review of managed honey bee populations in Europe and the United States and the factors that may affect them . Journal of Invertebrate Pathology 2010 ; 103 . doi: 10.1016/j.jip.2009.06.011 OpenUrl CrossRef PubMed 71. ↵ Büchler R et al. Standard methods for rearing and selection of Apis mellifera queens 2.0 . Journal of Apicultural Research ; 0 : 1 – 57 . doi: 10.1080/00218839.2023.2295180 OpenUrl CrossRef 72. ↵ Brochet S et al. Niche partitioning facilitates coexistence of closely related honey bee gut bacteria . eLife 2021 ; 10 : e68583 . doi: 10.7554/eLife.68583 OpenUrl CrossRef PubMed 73. Dodge R et al. A symbiotic physical niche in Drosophila melanogaster regulates stable association of a multi-species gut microbiota . Nat Commun 2023 ; 14 : 1557 . OpenUrl CrossRef PubMed 74. Meier DV et al. Niche partitioning of diverse sulfur-oxidizing bacteria at hydrothermal vents . ISME J 2017 ; 11 : 1545 – 1558 . doi: 10.1038/ismej.2017.37 OpenUrl CrossRef 75. ↵ Caesar L et al. Spatial segregation and cross-kingdom interactions drive stingless bee hive microbiome assembly. 2025 . bioRxiv , 2025 . , 2025.03.07.642116 76. ↵ Yu MK , Fogarty EC , Eren AM . Diverse plasmid systems and their ecology across human gut metagenomes revealed by PlasX and MobMess . Nat Microbiol 2024 ; 9 : 830 – 847 . OpenUrl 77. ↵ Vannette RL , Gauthier M-PL , Fukami T . Nectar bacteria, but not yeast, weaken a plant– pollinator mutualism . Proceedings of the Royal Society B: Biological Sciences 2013 ; 280 : 20122601 . doi: 10.1098/rspb.2012.2601 OpenUrl CrossRef PubMed 78. ↵ Vannette RL . The Floral Microbiome: Plant, Pollinator, and Microbial Perspectives . Annu Rev Ecol Evol Syst 2020 ; 51 : 363 – 386 . OpenUrl CrossRef 79. ↵ Graystock P et al. The Trojan hives: Pollinator pathogens, imported and distributed in bumblebee colonies . Journal of Applied Ecology 2013 ; 50 : 1207 – 1215 . doi: 10.1111/1365-2664.12134 OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted October 31, 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 Nested Connections: Local Phage and Broad Plasmid Sharing in the Honey Bee Mobilome 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 Nested Connections: Local Phage and Broad Plasmid Sharing in the Honey Bee Mobilome Delaney L. Miller , Lílian Caesar , Carrie Ganote , Sergio López-Madrigal , Danny W. Rice , Chris R. P. Robinson , Alyssa M. Seawright , Jillian P. Lewis , Irene L. G. Newton bioRxiv 2025.10.31.685865; doi: https://doi.org/10.1101/2025.10.31.685865 Share This Article: Copy Citation Tools Nested Connections: Local Phage and Broad Plasmid Sharing in the Honey Bee Mobilome Delaney L. Miller , Lílian Caesar , Carrie Ganote , Sergio López-Madrigal , Danny W. Rice , Chris R. P. Robinson , Alyssa M. Seawright , Jillian P. Lewis , Irene L. G. Newton bioRxiv 2025.10.31.685865; doi: https://doi.org/10.1101/2025.10.31.685865 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 (7633) Biochemistry (17680) Bioengineering (13889) Bioinformatics (41927) Biophysics (21445) Cancer Biology (18585) Cell Biology (25491) Clinical Trials (138) Developmental Biology (13373) Ecology (19897) Epidemiology (2067) Evolutionary Biology (24308) Genetics (15606) Genomics (22494) Immunology (17736) Microbiology (40385) Molecular Biology (17175) Neuroscience (88583) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4822) Physiology (7641) Plant Biology (15149) Scientific Communication and Education (2045) Synthetic Biology (4293) Systems Biology (9822) Zoology (2271)