Distinct Roles of Honeybee Gut Bacteria on Host Metabolism and Neurological Processes

To reveal microbiota-induced circulating metabolome changes, we colonized newly emerged bees with six honeybee core gut bacteria, Gilliamella apicola (Gi), Bifidobacterium asteroides (Bi), Snodgrassella alvi (Sn), Lactobacillus Firm4 (F4) and Firm5 (F5), and Bartonella apis (Ba) (Fig. 1A and Fig. S1 in the supplemental material). Then, we performed quasi-targeted metabolomics analysis of hemolymph samples from gnotobiotic bees. In total, we identified 326 metabolites among different bee groups (Table S1 in the supplemental material). Orthogonal partial least squares discriminant analysis (OPLS-DA) showed that the metabolic signatures of hemolymph samples were significantly different between groups (Fig. 1B). To associate clusters of metabolites highly correlated with particular gut members, we performed weighted correlation network analysis (WGCNA). Bees inoculated with different gut microbes were used as the sample traits. Based on interaction patterns among the 326 metabolites, WGCNA clustered them into eight modules (M) (Table S2). Interestingly, only three modules (red, pink, and black) were significantly correlated with bees mono-colonized by Gilliamella, Lactobacillus Firm4, and Lactobacillus Firm5 (P < 0.01; Fig. 1C).

The relatedness of bacteria-colonized groups with significantly altered metabolites often comprises an intramodular connectivity within the network module. The red and pink M were significantly associated with the Gilliamella group. The metabolites that comprised the red and pink M were mainly involved in carbohydrate, glycerophospholipid, and amino acid metabolism (Fig. 2A). Specifically, the primary driving metabolites from the red M were glycerophosphocholine, choline, and glycerol-3-phosphocholine, which are involved in the glycerophospholipid metabolism pathways (Fig. 2B). In the pink M, the significantly related metabolites were sucrose, trehalose, galactinol, glucose, and isomaltose, which are enriched in carbohydrate metabolism (Fig. 2C). This is consistent with the potential of G. apicola for carbohydrate metabolism in the gut (17, 18). When we reanalyzed published metabolomic data obtained from the guts of mono-colonized bees (16), we found that certain pathways, including purine metabolism, glycolysis/gluconeogenesis, and alanine, aspartate, glutamate, glycine, serine, and threonine metabolism, were stimulated in both hemolymph and gut samples by the colonization of Gilliamella (Fig. 2D and Fig. S2A).

Both Lactobacillus Firm4 and Firm5 were significantly correlated with the black M of metabolites enriched in amino acid metabolism (Fig. 3A). Interestingly, the most-related circulating metabolite, homovanillic acid (Fig. 3B and C), is an index of brain dopamine metabolism for assessing neurologic and psychiatric illnesses, such as Parkinson’s and Huntington’s diseases (21–23). Pathways including lysine degradation and pyrimidine, tryptophan, and purine metabolism were consistently upregulated in the hemolymph and gut of F4- and F5-colonized bees (Fig. 3D and Fig. S2B–C). Altogether, these data indicate that distinctive gut members regulate specific metabolic pathways.

We have shown that metabolites in the circulatory system are regulated by honeybee gut bacteria. These small molecules can influence gene expression and neuronal function in the brain (24). Since the hemolymph metabolome was only modulated in bees colonized with Gi, F4, and F5, these three groups were used for brain transcriptome analysis. Our data revealed that 713 genes in total were differentially expressed in bees colonized with gut members compared to MF bees (|log₂-fold change| >1, false discovery rate [FDR] < 0.05, Table S3), and different bee groups exhibited distinctive brain gene expression profiles (Fig. S3A). Specifically, the heat shock protein gene (loc410620) was upregulated in the Gi group (Fig. 4A). The odorant-binding protein gene Obp14 was upregulated both in the F4 and F5 group (Fig. 4B and C). Insect odorant-binding proteins are essential for detection and distinguishing of specific odors (25). In addition, several mrjp family genes of the major royal jelly protein (MRJP) encoded in the A. mellifera genome were significantly upregulated in the F4 and F5 groups, while bees colonized with Gilliamella exhibited decreased expression of mrjp genes (Fig. 4A–C). MRJPs have polyfunctional properties and participate in all major aspects of eusocial behavior in honeybees, such as caste determination and age polyethism (26).

Enrichment analysis of differentially expressed genes determined that KEGG pathways, including those for alpha-linolenic, arachidonic acid, linoleic, vitamin, and glycerophospholipid metabolism, were upregulated in the brains of F4 and F5 groups (Fig. 4D). The glycolysis/gluconeogenesis and glucagon signaling pathway, which is critical for brain physiology as it provides the fuel for brain functions (27), was only upregulated in bees colonized with Gilliamella. A KEGG pathway involving RNA processing by the spliceosome was upregulated in Gi and F4 bees. These results showed that the transcriptomic programs were differentially altered in the bacteria-colonized groups.

We then performed Pearson correlation analysis to analyze the relationship between gene expression and the metabolites which are significantly regulated by Gi, F4, and F5 separately (Fig. S3B–D). In the Gi group, carbohydrate metabolites showed significantly positive correlations (P < 0.05) with upregulated genes in the brain, while metabolites involved in glycerophospholipid metabolic pathways showed significantly negative correlations with those genes (Fig. S3B). In the F4 and F5 groups, fexofenadine, cholesteryl sulfate, and 2′-hydroxy-5′-methylacetophenone showed significantly positive correlations with upregulated genes in the brain, while homovanillic acid and gamma-caprolactone showed significantly negative correlations with those genes (Fig. S3C–D). These results suggest that upregulated genes in the brain were correlated with specific metabolites altered by different gut bacteria, indicating the potential roles of circulating metabolisms in host-bacteria interactions.

Shpigler et al. (28) have reported that the gene expression signatures of honeybees with disordered social behaviors are significantly enriched for human autism spectrum disorder (ASD)-related genes. Likewise, the differentially expressed genes in bacteria-colonized bees overlapped with those from human ASD patients (Fig. 4E), implying the involvement of gut microbiota in host behaviors. However, the differentially expressed genes present in the SFARI data set did not overlap the gene set for bees with disordered behaviors, as defined by Shpigler et al. (28).

Voineagu et al. (29) have shown that dysregulation of alternative splicing (AS) in ASD-related genes is also associated with the psychiatric disorder. Given evidence for spliceosome alterations, we investigated whether gut bacteria-colonized bee brains showed different AS events compared with those of MF bees. rMATS analysis of alternative splicing events of brain genes detected a total of 18,985 events in 4,929 genes, and skipped exon (SE) was the most abundant of the different types of AS. About 10 to 25% of events for each type of AS showed significantly different inclusion rates in bacteria-colonized bees (Fig. S4A). The relative abundance of different AS event types was similar across bee groups (Fig. S4A). However, the UpSet plot shows that the vast majority of events do not intersect between sets, indicating that multiple AS events can occur in a single gene and that the gut members cause different events (Fig. S4B).

Next, we examined the overlap of genes showing significantly differential AS events between MF and bacteria-colonized bees with the ASD risk genes from the SPARK and SFARI gene data sets (30). Fifty-nine of the 2,128 differentially spliced genes in MF bees are associated with human autism (Fig. 5A). Interestingly, almost all identified homologs belong to the high-confidence SFARI gene list (category 1) implicated in ASD (Table S4). Specifically, genes related to the pathophysiology of ASD are differentially spliced in MF bees compared with bacteria-colonized bees. For example, genes belonging to the functional terms of postsynaptic membrane (Ank2, Fmr1) and voltage-gated ion channel (Scn1a) showed different inclusion rates, and these are mostly regulated in the Gi groups compared to in the MF bees (Fig. 5B). Taking these findings together, we identified that gut microbes induce differential gene expression profiles and mediate AS, resulting in specific gene isoforms in the honeybee brain. The genes essential for bee social behaviors and related to human ASD disease are affected especially by Gilliamella, confirming the similarities of genes associated with social responsiveness in humans and honeybees (28).

The enrichment analysis of differentially spliced genes identified that KEGG pathways, including the neurotrophin, MAPK, and calcium signaling pathways, as well as ion channels and the GABAergic, dopaminergic, and serotonergic synapses, were regulated in the brains of the Gi, F4, and F5 groups (Fig. 5C). Therefore, we investigated changes in the brain neurochemistry of MF and mono-colonized bees. The concentrations of three major neurotransmitters, GABA, 5-HT, and dopamine, which are important modulators of honeybee feeding behavior, were determined (31–33). The concentration of the inhibitory transmitter GABA was significantly higher in the brains of F4 bees (Fig. 6A). In contrast, the concentrations of dopamine and 5-HT were significantly lower in bees mono-colonized with Gi, F4, and F5 than in the MF bees (Fig. 6B and C).

So far, our results have illustrated the impact of gut bacteria on the brain transcriptomic profile and neurochemistry of the honeybee. Finally, we tested whether the colonization of specific gut species alters host behavior. Proboscis extension response (PER) is a taste-related behavior that is fundamental for olfactory discrimination and colony performance in honeybees (34). We then measured the PER of MF bees, conventional (CV) bees, and bees mono-colonized with Gi, F4, and F5 to determine whether specific gut members influence host behavior. Compared with that of the MF group, sucrose sensitivity was significantly elevated in bacteria-colonized bees (Fig. 6D). Interestingly, CV bees were more sensitive to low sucrose concentrations compared to mono-colonized bees, implying an integrative effect of gut bacteria on honeybee behavior. Our findings indicate that colonization of either the normal gut microbiota or single core gut member can affect neurotransmitter levels in the brain, which might be associated with altered olfactory sensitivity.

Microbiota-free bees were obtained as described by Zheng et al. (18) with modifications (Fig. S1). Late-stage pupae were removed manually from brood frames and placed in sterile plastic bins. The pupae emerged in an incubator at 35°C, with humidity of 50%. Newly emerged MF bees (day 0) were kept in axenic cup cages with sterilized sucrose syrup (50%, wt/vol) for 24 h and divided into three groups: (i) MF, (ii) mono-colonized (MC), and (iii) CV bees. For each setup, 20 to 25 MF bees (day 1) were placed into one cup cage and fed on the corresponding solutions or suspensions for 24 h. For the MF group, 1 mL of 1× phosphate-buffered saline (PBS) was mixed with 1 mL of sterilized sucrose solution (50%, wt/vol) and 0.3 g sterilized pollen. For the MC group, stocks of Gilliamella apicola (W8127), Snodgrassella alvi (W6238G3), Bifidobacterium asteroides (W8113), Bartonella apis (B10834G6), Lactobacillus Firm4 (W8089), and Lactobacillus Firm5 (W8172) in 25% glycerol stock at −80°C were resuspended in 1 mL 1×PBS (Solarbio, Beijing, China) at a final optical density at 600 nm (OD₆₀₀) of 1, and then mixed with 1 mL sterilized sucrose solution (50%, wt/vol) and 0.3 g sterilized pollen. For the CV group, 5-μL homogenates of freshly dissected hindguts of nurse bees from their hives of origin were mixed with 1 mL 1×PBS, 1 mL sterilized sucrose solution (50%, wt/vol), and 0.3 g sterilized pollen. Then, MF, MC, and CV bees were provided sterilized sucrose (0.5 M) with sterile pollen and kept in an incubator (35°C, RH 50%) until day 7. Brains, guts, and hemolymph of bees were collected on day 7 for further analysis. All bees used for behavior experiments and tissue collection came from the same colony. Bees coming from the same cup cage for each group were considered a single replicate. Here, we used n = 3 replicates per group, 20 to 30 bees per replicate. Colonization levels of MF and MC bees were determined by CFU from dissected guts, as described by Kwong et al. (46).

Whole guts were dissected by tweezers sterilized with 75% alcohol. Dissected guts were directly crushed in 25% (vol/vol) glycerol on ice for bacterial load quantification or collected into an empty 1.5-mL centrifuge tube for metagenomic sequencing and metabolomics analysis. All gut samples were frozen at −80°C until analysis. Honeybee brains were collected using a dissecting microscope (Canon). Individual bees were fixed on beeswax using two insect needles through the thorax. After removal of the head cuticle, the whole brain was placed on a glass slide and soaked in RNAlater (Thermo; Waltham, MA, USA) or proteinase inhibitor (Roche; Mannheim, Germany) for gene expression profiling, proteome analysis, and neurotransmitter concentration quantification. Then, the hypopharyngeal glands, salivary glands, three simple eyes, and two compound eyes were carefully removed. Dissected brains were kept frozen at −80°C. Hemolymph was collected using a 10-μL pipettor (Eppendorf; Hamburg, Germany) from an incision above the median ocellus. A minimum of 50 μL hemolymph was collected from 10 bees into a 1.5-mL centrifuge tube. During the collection process, tubes were temporarily preserved on dry ice and subsequently stored at −80°C until analysis. Brain and hemolymph samples were collected from different individual bees following the same colonization step.

Hemolymph metabolites were determined by quasi-targeted metabolomics using high-pressure liquid chromatography (HPLC)-tandem mass spectrometry (MS-MS). First, 50 μL of hemolymph samples were mixed with 400 μL prechilled methanol by vortexing. All samples were incubated on ice for 5 min and then centrifuged at 15,000 × g at 4°C for 10 min. The supernatant was diluted to a final concentration containing 53% methanol by liquid chromatography-mass spectrometry (LC-MS) grade water. The samples were then transferred to a fresh vial and centrifuged at 15,000 × g at 4°C for 20 min. Finally, the supernatant was injected into the LC-MS-MS system, and the analyses were performed using an ExionLC AD system (SCIEX) coupled with a QTRAP 6500+ mass spectrometer (SCIEX). Samples were injected onto a BEH C8 column (100 mm × 2.1 mm × 1.9 μm) using a 30-min linear gradient at a flow rate of 0.35 mL/min for the positive-polarity mode. Eluent A was 0.1% formic acid-water, and eluent B was 0.1% formic acid-acetonitrile. The solvent gradient was set as follows: 5% B, 1 min; 5 to 100% B, 24.0 min; 100% B, 28.0 min; 100 to 5% B, 28.1 min; 5% B, 30 min. The QTRAP 6500+ mass spectrometer was operated in positive-polarity mode with curtain gas at 35 lb/in², collision gas at medium, ion-spray voltage at 5,500 V, temperature at 500°C, ion source gas at 1:55, and ion source gas at 2:55. For negative-ion mode, samples were injected onto aHSS T3 Column (100 mm × 2.1 mm) using a 25-min linear gradient at a flow rate of 0.35 mL/min. The solvent gradient was set as follows: 2% B, 1 min; 2 to 100% B, 18.0 min; 100% B, 22.0 min; 100 to 5% B, 22.1 min; 5% B, 25 min. The QTRAP 6500+ mass spectrometer was operated in negative polarity mode with curtain gas at 35 lb/in², collision gas at medium, ion-spray voltage at −4,500 V, temperature at 500°C, ion source gas at 1:55, and ion source gas at 2:55.

Detection of the experimental samples using MRM was based on the Novogene in-house database. Q3 (daughter) was used for the quantification. Q1 (parent ion), Q3, retention time, declustering potential, and collision energy were used for metabolite identification. Data files generated by HPLC-MS/MS were processed with SCIEX OS (version 1.4) to integrate and correct the peaks. A total of 326 compounds were identified in the hemolymph samples. Metabolomics data analysis was then performed using MetaboAnalyst 4.0 (47).

R software package WGCNA 1.69 (48) was used to identify key phenotype-related metabolic modules based on correlation patterns. The Pearson correlation matrix was calculated for all possible metabolite pairs and then transformed into an adjacency matrix with soft thresholding power set to 5 for the best topological overlap matrix. A dynamic tree cut algorithm was used to detect groups of highly correlated metabolites. The minimum module size was set to 14, and the threshold for merging module was set to 0.25 as default. Each module was assigned a unique color and contained a unique set of metabolites. After modules had been obtained from each group, module eigenmetabolite was calculated with the “ModuleEigengenes” function. Association analysis between a module and the trait of each group was performed using the “corPvalueStudent” function based on the module eigenmetabolite. Statistical significance was set at P < 0.01. Metabolites in each module were annotated on the KEGG database and classified into major categories using MetaboAnalyst 4.0 (47) for enrichment analysis. Finally, the network connections among metabolites in modules were visualized using Cytoscape 3.7.0 (49).

Total RNA was extracted from individual brains using the Quick-RNA MiniPrep kit (Zymo; Irvine, CA, USA). RNA degradation and contamination were monitored on 1% agarose gels, and the purity was checked with the NanoPhotometer spectrophotometer (IMPLEN; CA, USA). RNA integrity was assessed using the RNA Nano 6000 assay kit of the Bioanalyzer 2100 System (Agilent Technologies; Santa Clara, CA, USA). RNA sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England BioLabs; Ipswich, MA, USA), and index codes were added to attribute sequences to each sample. Clustering of the index-coded samples was performed on the cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina; San Diego, CA, USA), and the library preparations were then sequenced on an Illumina NovaSeq 6000 platform (Illumina; San Diego, CA, USA) and 150-bp paired-end reads were generated.

The sequencing quality of individual samples was assessed using FastQC version 0.11.5 with default parameters. An index of the bee reference genome (Amel_HA version 3.1) was built using HISAT2 version 2.0.5 (50), and the FastQC trimmed reads were then aligned to the built index using HISAT2 version 2.1.0 with default parameters. Gene expression was quantified using HTSeq version 0.7.2 (51) with ‘union’ mode, and only reads mapping unambiguously to a single gene were counted. In contrast, reads aligned to multiple positions or overlapping more than one gene were discarded. If it were counted for both genes, the extra read from the differentially expressed gene might have caused the other gene to be wrongly identified as differentially expressed, so we chose ‘union’ mode.

Differential gene expression analysis was performed using the DESeq2 package (52) in R. We modeled read counts following a negative binomial distribution with normalized counts and dispersion. The proportion of the gene counts in the sample to the cDNA concentration was scaled by a normalization factor using the median-of-ratios method. The variability between replicates was modeled by the dispersion parameter using empirical Bayes shrinkage estimation. For each gene, we fit a generalized linear model to obtain the overall expression strength of the gene and the log₂-fold change between the MC and MF groups. For significance testing, differential gene expression was determined by the Wald test. The resulting P values were corrected for multiple comparisons using the Benjamini-Hochberg FDR method (53). Genes with an adjusted P value of <0.05 and a |log₂-fold change| of >1 were assigned as differentially expressed.

To get a better annotation of the honeybee reference genome, we reannotated it using eggNOG-mapper version 5.0 (54). A total of 6,269 out of 12,375 honeybee genes were successfully assigned to a KO entry with ‘diamond’ mode, and hierarchy information for the KEGG metabolic pathway was extracted. Functional analysis of differentially expressed genes was performed based on KEGG Orthologue (KO) markers. The percentages of KO markers belonging to each category (KEGG class at level 3) out of the total MC- and MF-enriched KO markers were designated as comparison parameters. The significance level was calculated by Fisher’s exact test using clusterProfiler version 3.10.1 (55). Correlation analysis between regulated hemolymph metabolites and upregulated genes was performed using Pearson’s correlation analysis. The correlation coefficients and the corresponding P values were calculated in R.

Analysis of event-level differential splicing was performed using rMATS version 4.0.2 (56) based on the bee reference genome. An exon-based ratio metric, commonly defined as percent-spliced-in (PSI) value, was employed to measure the alternative splicing events. An effective length of l was used for normalization. The PSI value was calculated for several classes of alternative splicing events, including skipped exon (SE), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), mutually exclusive exons (MXE), and retained introns (RI). Events with a P of <0.05 were considered to be differentially spliced across gnotobiotic bees and microbiota-free bees.

To find overlaps between the differentially expressed or spliced genes of bee brain and those from human autism spectrum disorders, a total of 3,531 high-quality reference protein sequences corresponding to 948 known autism risk genes (SFARI: https://gene.sfari.org/, SPARK for autism: http://spark-sf.s3.amazonaws.com/SPARK_gene_list.pdf) were aligned against protein sequences of the honeybee genome using BLASTP (57) with a two-way best matching strategy. In total, 649 autism protein sequences obtained a match (similarity > 30% and E value < 0.000394). Then we calculated the intersection of the autism risk genes and the differentially expressed or spliced genes between bacteria-colonized and MF bees (P < 0.05).

Brain tissues dissected from MF and MC bees were sent to Biotree Biotech Co. Ltd. (Shanghai, China) for targeted metabolomics analysis of GABA, dopamine, and serotonin. Six brain tissues from one treatment group were put into one tube and centrifuged (2,400 g × 1 min at 4°C). A 100-μL volume of acetonitrile containing 0.1% formic acid and 20 μL ultrapure water was added and the tubes were vortexed thoroughly. Metabolites were sonicated in an ice-water bath for 30 min, followed by subsiding at −20°C for 2 h. Supernatants were collected after centrifugation (14,000 g × 10 min at 4°C). Next, 20 μL of supernatant was transferred to a new vial followed by incubation for 30 min after the addition of 10 μL sodium carbonate solution (100 mM) and 10 μL 2% benzoyl chloride acetonitrile. Then, 1.6 μL internal standard and 20 μL 0.1% formic acid was added and the samples were centrifuged (14,000 g × 5 min at 4°C). A total of 40 μL of the supernatants was transferred to an auto-sampler vial for downstream UHPLC-MS/MS analysis. Serotonin hydrochloride, γ-aminobutyric acid, and dopamine hydrochloride (Aladdin Bio-Chem Technology; Shanghai, China) derivatized with benzoyl chloride (Sigma-Aldrich; Saint Louis, MO, USA) were used for the construction of the calibration standard curve. The internal standards mixtures (γ-aminobutyric acid, dopamine hydrochloride, and serotonin hydrochloride) derivatized with benzoyl chloride-d5 (Sigma-Aldrich; Saint Louis, MO, USA) (58) of the corresponding concentrations were prepared.

The UHPLC separation was carried out using an ExionLC System (AB SCIEX; Framingham, MA, USA), and samples were analyzed on the QTRAP 6500 LC-MS/MS system (AB SCIEX). Two μL of samples was directly injected onto an ACQUITY UPLC HSS T3 column (100 mm × 2.1 mm × 1.8 μm; Waters Corporation; Milford, Ma, USA). The column temperature was set at 40°C, and the auto-sampler temperature was set at 4°C. Chromatographic separation was achieved using a 0.30 mL/min flow rate and a linear gradient of 0% to 2% B in 2 min; 2 to 98% B in 9 min, 98% B for 2 min, and equilibration for 2 min. Solvent A was 0.1% formic acid and solvent B was acetonitrile. For all multiple-reaction monitoring (MRM) experiments, 6500 QTrap acquisition parameters were as follows: ion-spray voltage at 5,000 V, curtain gas setting of 35, nebulizer gas setting of 60, and temperature at 400°C. Raw data were analyzed using Skyline (59).

Seven-day-old MF, MC, and CV bees were used to measure responses to different concentrations of sucrose solution as previously described, with some modifications (9). Before the test, bees were starved for 2 h in the incubator by removing sugar syrup and bee bread from the cup cage. Bees were then mounted to modified 2.0-mL centrifuge tubes using Parafilm M (Bemis; Sheboygan Falls, WI, USA), and they could only move their heads and propodeum for antennae sanitation. Individual responsiveness was measured by presenting a series of sucrose solution concentrations (0.1, 0.3, 1, 3, 10, 30, and 50%; wt/vol) to the antennae of bees (60). Before each sucrose solution presentation, all bees were tested for their response to pure water to control for the potential effects of repeated sucrose stimulations, which could have led to either sensitization or habituation (61). The interstimulus interval between water and sucrose solution was 4 min. When a bee’s antenna is stimulated with a sucrose solution of sufficient concentration, the bee reflexively extends its proboscis. The lowest sucrose concentration at which an individual responded by extending its proboscis was recorded and interpreted as its sugar response threshold. At the end of the experiment, a gustatory response score was obtained for each bee based on the sucrose concentrations to which the bees responded. The response was arbitrarily quantified with scores of 1 to 7, where 1 represented a bee that only responded to the highest sucrose concentration and 7 represented an individual that responded to all concentrations tested. Bees that did not respond to any of the sucrose concentrations were excluded from further analyses. Additionally, bees that responded to all concentrations of sucrose solutions and all water presentations were also excluded, as they appeared to be unable to discriminate between sucrose and water (61).

Raw sequence reads have been deposited in the NCBI SRA database under the BioProject accession no. PRJNA783076.