Social status shapes the bacterial and fungal gut communities of the honey bee

Abstract

Despite the fungal abundance in honey and bee bread, little is known about the fungal gut community of the honey bee and its effect on host fitness. Using pyrosequencing of the 16S rRNA gene and ITS2 region amplicons, we analysed the bacterial and fungal gut communities of the honey bee as affected by the host social status. Both communities were significantly affected by the host social status. The bacterial gut community was similar to those characterised in previous studies. The fungal gut communities of most worker bees were highly dominated by Saccharomyces but foraging bees and queens were colonised by diverse fungal species and Zygosaccharomyces, respectively. The high fungal density and positive correlation between Saccharomyces species and Lactobacillus species, known yeast antagonists, were only observed in the nurse bee; this suggested that the conflict between Saccharomyces and Lactobacillus was compromised by the metabolism of the host and/or other gut microbes. PICRUSt analysis revealed significant differences in enriched gene clusters of the bacterial gut communities of the nurse and foraging bees, suggesting that different host social status might induce changes in the gut microbiota, and, that consequently, gut microbial community shifts to adapt to the gut environment.

Figure 1

Comparison of alpha-diversity indices of the gut microbial community in honey bees performing four social roles. Box plots depict the medians (central horizontal lines), inter-quartile ranges (boxes), and 95% confidence intervals (whiskers). NEB, newly-emerged bee; 12-h, 12-h-old bee; NB, nurse bee; FB, foraging bee. An FDR adjusted p-value are from Kruskal-Wallis test. Asterisks indicate statistically significant differences between pairs of values (*P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001).

Figure 2

Composition of the gut microbial community. A shift in the gut microbial community, as affected by the honey bee social status, is represented (a) at the resident bacterial taxa and (b) at the major fungal genus level (>0.5% of all sequences). NEB, newly-emerged bee; 12-h, 12-h-old bee; NB, nurse bee; FB, foraging bee; Q, queen; RJ, royal jelly; Bo, Bombus.

Figure 3

Changes in the gut microbial community depending on the honey bee social status. (a) The honey bee gut bacterial and fungal communities clustered using PCoA of the unweighted Jaccard-based and weighted thetaYC-based matrices. Group names are designated by initials, with different colours representing categories, described in Table S1. Nonparametric ANOVA tests were used to test (b) intra-group dissimilarity and (c) inter-group dissimilarity based on the unweighted Jaccard- and weighted thetaYC-based distances. Box plots depict medians (central horizontal lines), the inter-quartile ranges (boxes), 95% confidence intervals (whiskers), and outliers (black dots). Upper panel, bacteria; lower panel, fungi. NEB, newly-emerged bee; 12-h, 12-h-old bee; NB, nurse bee; FB, foraging bee. An FDR adjusted p-value are from Kruskal-Wallis test. Asterisks indicate statistically significant differences between pairs of values (*P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001).

Figure 4

Comparison of the predictive functions of bacterial communities using 16S rRNA marker gene sequences. (a) PCA representing changes in the total gene family at the KEGG level 3 of bacterial communities, according to the honey bee social roles. (b) Comparison of the profiles of differentially enriched genes between the nurse bee and foraging bee. The threshold of the LDA score was 3.0. NEB, newly-emerged bee; 12-h, 12-h-old bee; NB, nurse bee; FB, foraging bee. P-values are carried out with unpaired t-test with Welch’s correction. Asterisks indicate statistically significant differences between the pairs of values (*P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001).

Figure 5

Microbial density associated with the social role of the honey bee. Comparison of copy number of 16S rRNA and ITS2 region (a) of individuals and (b) of groups. The copy numbers were estimated by qPCR. (c) Normalization of the 16S rRNA gene copies to ITS2 region copies (ratios) in individual samples. (d) The ratio of 16S rRNA gene copies to ITS2 region copies analysed according to the host social role. NEB, newly-emerged bee; 12-h, 12-h-old bee; NB, nurse bee; FB, foraging bee; Q, queen; RJ, royal jelly; Bo,Bombus; NTC, no template control. Data are expressed as the mean ± s.e.m. An FDR adjusted p-value are from ANOVA after Kruskal-Wallis post-hoc test. Blue and red asterisks refer to statistically significant differences in the bacterial and fungal density, respectively. An FDR adjusted p-value; *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.