We first investigated the impact of the banana and bee bread on beetle survival in the lab condition. After 3 weeks of rearing, 30 beetles survived in the Banana group, and 27 survived in the Bee_Bread group (Supplementary Tables S1), hence the variance of the banana and bread on SHBs survival was minor (Pearson’s Chi-squared test, df = 2, P = 0.516).

As the microbes associated with SHB are largely unknown, we then assembled the microbiome genomes. We randomly selected 10 SHBs in the Banana and Bee_Bread groups and eight SHBs collected in bee hive for metagenomic analyses to quantify the associated gene content and microbes. On average, 5.5 ± 1.3 million reads (150 bp per read) were assigned to microbes in each library, and the ratios of microbe to host reads did not differ significantly between libraries (t-test, P > 0.05, Supplementary Table S2). We assembled 196,214 metagenome contigs (Supplementary File S2). In those contigs, 436,918 genes were predicted, and 200,587 protein-coding genes retrieved functional annotation in KEGG (Supplementary File S3-S4; Figure 1). The symbionts showed a substantial number of genes involved in nucleotide, energy, carbohydrate, and amino acid metabolism. We also noticed that the symbionts harbor a few xenobiotics biodegradation genes, presumably to degrade toxic compounds (Table 1; Supplementary Table S3).

To quantify the gene copy number, we aligned the reads back to the metagenome. In a pairwise analysis, 15,748 genes were significantly differentially represented between the Banana and Wild groups, 10,678 genes between the Banana and Bee_Bread groups, and 7,994 genes between the Wild and Bee_Bread groups (Table 2; Supplementary Table S4). The Banana group showed a substantially higher number of over-represented genes than the Bee_bread and Wild groups in carbohydrate, lipid, and amino acid metabolism (Chi-squared test, FDR < 0.001). Among the differentially represented genes, 143 carbohydrate-active enzymes (CAZy) were over-represented in the Banana group, compared with 130 in the Wild group and 35 in the Bee_Bread group (Supplementary Table S5). The Banana group showed more CAZy than random among the three treatment groups (Pearson’s Chi-squared test, df = 2, p < 0.001).

To determine microbial diversity, we aligned the microbial reads to the Kraken2 standard database (Supplementary Figure S2-S4). Despite high symbiont diversity overall, 13 bacterial genera dominated the microbial community (> 99% of relative abundance, Figure 2A). The diet separated the microbes and accounted for 87% of the variance (Figure 2B). The Banana group maintained 95 microbial species/strains, followed by 67 in the Bee_Bread group and 56 in the Wild group (Figure 3A). The Banana group showed the highest alpha diversity (2.1 ± 0.31), compared with the Bee_Bread group (1.83 ± 0.40) and the Hive group (1.70 ± 0.30; Kruskal-Wallis test, df = 2, P < 0.05, Figure 3B). We also compared the beta diversity between the paired groups. The Wild and Banana group showed the highest beta diversity (0.425), followed by 0.301 between the Wild and Bee_Bread groups. The Banana and Bee_Bread showed the lowest beta diversity of 0.196 (Table 3).

We found plant-associated bacteria in the Banana group, such as Corynebacterium glyciniphilum, which was the same species initially isolated from the banana for fermentation. The Banana group was significantly enriched in the acetic acid bacteria. Microbes associated with fermentation (Gluconobacter albidus, Mammaliicoccus sciuri, Corynebacterium nuruki) were also highly enriched in the Bee_bread group. We found a few honey bee symbionts in the Wild group, including Lactococcus lactis and the core symbiont Snodgrassella alvi (Supplementary File S5-S7).

We further tested whether SHB can support the colonization of honey bee gut symbionts. We inoculated newly emerged beetles with a genetically tagged clone of the symbiont S. alvi (wkB2:pBTK570). All the beetles survived by the end of the experiment. Tagged symbionts were found in all inoculated beetles, with the average CFU (colony forming unit) of 5,380 at 1 dpi (day post inoculation) on selective plates. Even though the CFU increased to 18,300, the colonization rate dropped to 20% at 3 dpi. When counting the CFU, we observed some microbial colonies morphologically different from S. alvi wkB2:pBTK570 in all beetles. To identify the microbial species, we sequenced one such colony, and the ITS2 region aligned with K. ohmeri (97% identity, P = 7e-97), forming a cluster in the phylogenetic tree (Supplementary Figure S5).

We directly collected adult SHBs from collapsed beehives (Apis cerana) in a commercial apiary. These beetles were maintained in an incubator (28 ± 1°C temperature and 65% ± 10% humidity), and then we transferred their pupal offspring to a new container until hatching (Neumann et al., 2013). We collected pollen from the honey bee Apis mellifera hive entrance. We randomly assigned newly hatched SHBs into two diet groups. SHBs fed on bananas were defined as the Banana group (N = 45), and those fed with simulated bee bread (an equal mix of pollen and 50% w/v sugar water to avoid bee hive microbes) comprised the Bee_Bread group (N = 45). We then assigned each of the 15 SHBs to a rearing tube. We additionally collected 13 SHBs from a beehive as the wild group because these SHBs were directly from the bee hive without lab feeding (N = 13). We refreshed the banana and bee bread daily and collected SHBs after 3 weeks of feeding. We rinsed SHB surfaces with distilled water, then extracted genomic DNA from individual SHBs using MagPure Soil DNA KF Kit (MP Biomedicals, USA). We prepared DNA sequencing libraries using the TruSeq Nano DNA LT Sample Preparation Kit (Illumina, USA). Ten adult SHBs in the Banana group, 10 in the Bee_Bread group, and eight in the Wild group were randomly selected for metagenomic sequencing using the Illumina NovaSeq 6,000 Platform, generating 150 bps paired-end reads.

First, we filtered the sequencing reads using Fastp (Version 0.23.2) with default parameters (Chen et al., 2018). Then we aligned the reads to the small hive beetle genome assembly (GCA_024364675.1) using BWA (Version 0.7.17-r1188) with default parameters and retrieved the unmapped reads using samtools (Version 1.7) and converted the bam to fastq file using bedtools (Version 2.26.0) (Li et al., 2009; Li and Durbin, 2009; Quinlan and Hall, 2010). After that, the unmapped reads were concatenated for all samples to assemble contigs using Megahit (Version 1.2.9) with default parameters (Li et al., 2015), and the contigs were further collapsed using redundans (Version 2020.01.28) (Pryszcz and Gabaldón, 2016). The genes were predicted using MetaGeneMark2 with default parameters.¹ The protein sequences were queried with the eggNOG-mapper and KEGG database to retrieve the putative function (Cantalapiedra et al., 2021). The code is provided in the Supplementary File S1.

The unmapped reads in each beetle were re-aligned to the assembled meta-contigs using BWA (Version 0.7.17-r1188) with default parameters. The number of reads aligned to each gene was quantified using bedtools (Version 2.26.0) (Quinlan and Hall, 2010). The number of aligned reads was normalized to the library size, and the over-represented genes were calculated with edgeR (Quinlan and Hall, 2010). The code is provided in the Supplementary File S1.

We performed two steps to bin the sequencing reads to the microbial species. We first aligned the assembled contigs to the most closely related microbes using the BusyBee tool (Laczny et al., 2017). Additionally, we aligned the reads to the Kraken2 (Version 2.1.2) standard database (built on 12/9/2022) (Wood et al., 2019). The number of reads assigned to each microbe was normalized using bracken (Version 2.8) (Lu et al., 2017). The relative abundance of the microbial species was used to calculate alpha (Shannon’s alpha diversity) and beta (Bray–Curtis dissimilarity) diversity using KrakenTools (Lu et al., 2022). The code is provided in the Supplementary File S1.

As the bee symbionts were constantly reported from SHB metagenomes, we inoculated SHBs with a genetically tagged bee symbiont to validate its colonization. The honey bee symbiont S. alvi wkB2:pBTK570 (Addgene accession ID#110615) was previously engineered to be spectinomycin resistant and stored at −80°C (Leonard et al., 2018). This symbiont isolate was activated on Columbia Blood Agar Base (Difco™, 279,220) with 5% sheep blood for 72 h, after which bacterial cells were diluted in 1000 μL PBS and adjusted to OD₆₀₀ = 1. Then, we mixed the homogenized bacterial cells with filter-sterilized 50% sucrose (50%) at a 1:1 ratio. SHBs were fed sucrose with tagged S. alvi wkB2:pBTK570 (N = 50), then SHBs were rinsed in ethanol and homogenized individually. Homogenates were plated on Columbia Blood Agar Base with 5% sheep blood and spectinomycin (60 μg/mL) to count the Colony Forming Unite (CFU) at 1, 3, 5, and 7 days post-inoculation. Detailed procedures are described in the Supplementary File S1.

We observed microbes morphologically distinct from S. alvi in all CBA plates. We collected a colony using an inoculating loop to extract DNA using the Qiagen DNeasy Plant Mini kit (cat#69104). The DNA was amplified using fungal ITS primers (FungITS.F 5’GTTAAAAAGCTCGTAGTTG3’; FungITS.R5’CTCTCAATCTGTCAATCCTTATT 3′) in a 30ul reaction consisting of 0.2ul Taq DNA polymerase (Invitrogen, 18038–240), 0.4ul primers, 0.2ul 10uM dNTP mix. The reactions were run with the cycling parameters: 94°C for 3 min., followed by 35 cycles of (94°C 15 s., 54°C 30s., 72°C 1 min.), 72°C for 5 min, and maintained in 4°C. The products were visualized on a 1.75% agarose gel, producing a single band. PCR products were sent for PCR clean-up and Sanger sequencing at Azenta Life Sciences, Rockville, Maryland. Sequences were searched in the NCBI Blastn MegaBlast. An additional nine sequences of fungal species were downloaded from NCBI and aligned with MUSCLE (Version 5.1) with default parameters. The tree was built with MrBayes (Version 3.2.7) and viewed with FigTree (Version 1.4.4). Detailed procedures are described in the Supplementary File S1.

We performed all statistics with R (Version 4.2.2) in RStudio (Version 2022.12.0) (R Core Team, 2013; RStudio Team, 2020). We compared surviving SHBs using Pearson’s Chi-squared test. The number of differentially enriched genes among the paired comparisons was analyzed with Pearson’s Chi-squared test and viewed with the VennDiagram package (Gao et al., 2021). The number of genes in metabolic categories was compared with random using a Chi-squared test, adjusted with the false discovery rate (FDR). The alpha diversity was first analyzed using the Kruskal-Wallis rank sum test, followed by pairwise comparisons using the Wilcoxon rank sum exact test, adjusted with FDR.