The gut microbiome defines social group membership in honey bee colonies

See allHide authors and affiliations

Science Advances  14 Oct 2020:
Vol. 6, no. 42, eabd3431
DOI: 10.1126/sciadv.abd3431


In the honey bee, genetically related colony members innately develop colony-specific cuticular hydrocarbon profiles, which serve as pheromonal nestmate recognition cues. Yet, despite high intracolony relatedness, the innate development of colony-specific chemical signatures by individual colony members is largely determined by the colony environment, rather than solely relying on genetic variants shared by nestmates. Therefore, it is puzzling how a nongenic factor could drive the innate development of a quantitative trait that is shared by members of the same colony. Here, we provide one solution to this conundrum by showing that nestmate recognition cues in honey bees are defined, at least in part, by shared characteristics of the gut microbiome across individual colony members. These results illustrate the importance of host-microbiome interactions as a source of variation in animal behavioral traits.


The integrity and robustness of social groups require the ability of individuals to reliably recognize and exhibit cues for group membership. Via processes analogous to how multicellular organisms perform self-recognition at the cellular and molecular levels to identify and prevent invasion by pathogens, animal social groups have evolved mechanisms of group member recognition to prevent loss of resources to invading individuals such as unrelated conspecifics, predators, and parasites (1, 2). However, the mechanisms that allow for member recognition at the social group level are not well understood (2).

Social insect colonies represent an extreme form of social lifestyle, which makes them an ideal system for studying mechanisms that regulate and maintain group membership and integrity. Social insects often rely on colony-specific blends of cuticular hydrocarbons (CHCs) as pheromonal cues that signify colony membership (3). Yet, despite the high genetic relatedness between colony members, nestmate recognition cues in many social insect species seem to be defined independently of colony-specific genetic variants. Instead, they are largely defined by factors derived from the colony and/or shared social environment (35). Current models suggest that nestmate recognition cues arise via “gestalt” mechanisms, whereby individuals produce CHCs that are then transferred and homogenized between individuals, such that every colony member carries a mean CHC profile (57). However, recent evidence in some species suggests that colony-specific CHC profiles are innately produced by the insect, rather than a result of a colony-level gestalt (4). Thus, it is puzzling how environmentally derived factors could drive an innate biosynthetic pathway to produce a robust colony-specific CHC signature that is shared by colony members, independent of their genetic relatedness. Here, we test the hypothesis that colony-specific nestmate recognition cues in honey bee colonies are defined by colony-specific gut microbial communities.

The concept of a “holobiont,” a distinct biological entity that is composed of the host and its microbiome, is emerging as an important and ubiquitous phenomenon in animal biology (8, 9). Under this concept, population-level genetic variation is a product of variations in the genomes of both animal hosts and their associated microbes (8). Thus, observed variations in host phenotypes, such as recognition cues, may actually result from variations in microbiomes. Studies in diverse animal species support the role of symbiotic microorganisms in driving chemical recognition cues of their hosts (1013). Therefore, it is possible that differences in the chemical profiles of individuals from different social insect colonies are primarily driven by community-level variations in microbiota across colonies. This provides a possible mechanistic explanation for why the development of nestmate recognition cues may seem to follow the rules of a genetically determined trait, yet does not rely on the actual relatedness among honey bee colony members.


To begin to test whether colony-specific nestmate recognition cues in honey bee colonies are defined by colony-specific gut microbial communities, we assessed whether forager bees from different honey bee colonies have different gut microbial communities and CHC profiles by performing 16S ribosomal RNA (rRNA) sequencing on gut samples and gas chromatography on cuticular extracts, respectively. We found that foragers from different colonies differ in overall gut microbial community structure [Fig. 1A, permutation multivariate analysis of variance (MANOVA) using Bray-Curtis (BC) dissimilarity index, F2,29 = 2.00, R2 = 0.13, P = 0.014] and CHC profile (Fig. 1B, permutation MANOVA using BC, F2,23 = 8.54, R2 = 0.45, P < 0.001). In addition, while foragers from different colonies largely overlap in gut microbial taxonomic diversity (Fig. 1C, ANOVA using the Shannon diversity index, F2,27 = 0.19, P = 0.83), they do vary in the abundance of individual gut microbes (Fig. 1D and table S1). Next, we sought to establish a causal relationship between the gut microbiome and CHC profile of individuals, independent of the bee’s genetic background. To do this, we used a cross-hive fostering experiment, where groups of newly eclosed bees sourced from two different colonies were raised in either their natal or an unrelated host colony. These studies revealed that both source- and host colony–related factors contribute to variations in the overall gut microbial community of individual bees (Fig. 1E, two-way permutation MANOVA using BC, host: F1,33 = 3.05, R2 = 0.08, P = 0.002; source: F1,33 = 2.46, R2 = 0.07, P = 0.011; and host*source: F1,33 = 1.80, R2 = 0.05, P = 0.072), as we have previously shown for CHC profiles (4). Nonetheless, when we compared abundances of specific microbial phylotypes across individual bees, we found that of the 14 microbial taxa that significantly differed between the four treatments, 6 were similar between bees that shared the same hive environment during behavioral maturation, regardless of genetic relatedness (table S2). In contrast, we did not identify any taxa that were significantly associated with source hive environment (table S2), suggesting that the posteclosion environment is the most important factor in defining the gut microbial community of adult worker bees. Together, these data indicate that although the honey bee gut microbial community consists of a core set of bacterial phylotypes that are common across colonies (14), individuals from different colonies exhibit varying relative abundances of these specific phylotypes, which may drive intercolony differences in CHC profiles.

Fig. 1 Gut microbial communities differ between foragers from different colonies.

(A to C) Forager honey bees across different colonies differ in overall gut microbial community (A) and CHC profile (B), but not in taxonomic richness (C). (D) Individual gut microbial taxa vary in abundance between forager bees from different colonies. (E) Forager bee gut microbial community is determined by colony environment. (A), (C), and (E) shown as nonmetric multidimensional scaling plots depicting BC dissimilarity between samples. (D) Identified microbial taxa shown as stacked bar plots depicting the average relative abundance of honey bee–associated microbes within forager guts from each colony. Bold letters in legends denote statistical significance, P < 0.05 (permutation MANOVA followed by FDR pairwise contrasts). ns, not significant (ANOVA).

To directly test the role of the gut microbiome in defining the CHC profiles of individual worker bees, we next assessed whether experimentally manipulating the gut microbial community between sister honey bees would be sufficient to drive differences in their CHC profiles. We first examined the overall role of the gut microbiome in regulating the CHC profile via various whole gut microbial community manipulations. As has been previously reported for other insects (10), we observed an effect of antibiotic treatment on honey bee gut microbiota—as determined by culturing gut samples on agar plates (Fig. 2, A and B)—and CHC profile (Fig. 2C, permutation MANOVA, F1,15 = 3.66, R2 = 0.21, P = 0.005). Similarly, when sister bees were inoculated with either live (full microbiome) or heat-killed inoculums, they developed different gut microbial communities (Fig. 2D, permutation MANOVA using BC, F1,39 = 6.93, R2 = 0.15, P = 0.001, and table S3) and CHC profiles (Fig. 2E, permutation MANOVA using BC, F1,31 = 6.17, R2 = 0.17, P = 0.006). Furthermore, because newly eclosed bees acquire their gut microbiota via interactions with older bees (15), we investigated the impact of raising day-old bees from a single source colony with older bees from different colonies (13). This revealed a significant effect on the experimental bees’ overall gut microbial communities (Fig. 2F, permutation MANOVA using BC, F1,19 = 8.60, R2 = 0.32, P = 0.001, table S4) and CHC profiles (Fig. 2G, permutation MANOVA using BC, F1,15 = 5.44, R2 = 0.28, P = 0.01). However, when the older bees were first treated with antibiotics, we did not observe a cohousing effect on the CHC profiles of the experimental bees (Fig. 2H, permutation MANOVA using BC, F1,13 = 1.55, R2 = 0.11, P = 0.194). In addition, we were unable to amplify the 16S rRNA gene in these experimental bees’ guts, indicating a negligible amount of gut bacteria in these samples and demonstrating that experimental bees are unable to acquire microbiomes when older bees are axenic. Together, these data indicate that acquiring different gut microbiomes, via social interactions with older bees, is sufficient to induce differences in CHC profiles between related individuals.

Fig. 2 Sister bees inoculated with different honey bee gut microbial communities develop different CHC profiles.

(A and B) Antibiotic-treated bee guts (A) show less growth when plated than those treated with a control sugar solution (B). (C) Sister bees differ in CHC profile when they are treated with antibiotics versus control. (D and E) Sister bees differ in gut microbial community (D) and CHC profile (E) when inoculated with live, full microbiome inoculum versus heat-killed inoculum. (F to H) Sister bees differ in gut microbial community (F) and CHC profile (G) and when they are inoculated by older bees from two different colonies, but not when the older bees are first treated with antibiotics (H). Data are shown as nonmetric multidimensional scaling plots depicting BC dissimilarity (C, E, G, and H) or weighted UniFrac (D and F) between samples. Bold letters in legends denote statistical significance, P < 0.05 (permutation MANOVA followed by FDR pairwise contrasts).

To further test the role of the microbiome in defining CHC profiles in individual worker bees, we manipulated the microbiomes of sister bees by inoculating them with different gut bacteria. Specifically, newly emerged bees were inoculated for 16 days with either Gilliamella apicola, a honey bee–specific symbiont, or Lonsdalea quercina, an opportunistic, environmentally acquired microbe (14, 16). Both of these gut bacteria were naturally present in the guts of foragers included in our studies and are easily culturable in standard LB (Luria broth) medium (Fig. 3A and table S5). Furthermore, G. apicola plays an important role in sugar metabolism and plant polysaccharide digestion in the bee gut (17, 18)—processes that result in the production of acetyl coenzyme A and short-chain fatty acids, which may serve as precursors to CHCs (10, 19). After inoculation with these two bacteria, we found that, in addition to developing different gut microbial communities (Fig. 3B, permutation MANOVA, F1,18 = 3.09, R2 = 0.15, P = 0.004) with different abundances of these specific bacteria (table S6), bees inoculated with these different bacteria developed different CHC profiles (Fig. 3C, permutation MANOVA, F1,15 = 7.33, R2 = 0.34, P = 0.004). Furthermore, in-lab behavioral acceptance assays showed that bees inoculated with the symbiotic G. apicola were able to distinguish sister bees inoculated with the same phylotype versus those inoculated with L. quercina (Fig. 3D, Pearson’s chi-square, χ2 = 7.59, P = 0.01). In contrast, bees inoculated with the opportunistic L. quercina were not able to discriminate between sister bees inoculated with the same phylotype versus those inoculated with G. apicola (Fig. 3E, Pearson’s chi-square, χ2 = 0.02, P = 1). This suggests that symbiotic, but not opportunistic, microbes impart the ability of individuals to generate and perceive nestmate recognition cues. Consequently, these findings imply a possible solution to the long-standing neuroethological puzzle of how a guard bee’s chemosensory system is specifically tuned to environmentally derived colony-specific chemical cues (3, 2025): Some symbiotic microbes may generate pleiotropic factors that physiologically couple the production and perception of nestmate recognition cues. While previous studies support the functional coupling of pheromone production and perception by genetic and hereditary factors at the evolutionary time scale (2631), our data suggest that coupling could also occur at the developmental and/or physiological time scales via, at least in part, the action of symbiotic gut microbes. Nevertheless, in combination with the experiments above, these findings suggest that acquiring different microbiomes is sufficient to drive the development of distinct CHC profiles between genetically related individuals.

Fig. 3 Gut microbiome plays a role in recognition in honey bees.

(A) G. apicola (black arrow) and L. quercina (green arrow) are easily culturable on LB plates in the lab. (B and C) Bees inoculated with G. apicola differ in gut microbial community (B) and CHC profile (C) from their sister bees inoculated with L. quercina. (D) Bees inoculated with G. apicola are able to distinguish sister bees inoculated with the same phylotype versus those inoculated with L. quercina. (E) Bees inoculated with L. quercina cannot distinguish between those inoculated with the same phylotype versus those inoculated with G. apicola. (F) Colony 1 and 2 bees inoculated with G. apicola or L. quercina differ in CHC profile based on inoculation treatment and colony of origin. (G) Colony 1 bees inoculated with G. apicola accept colony 1 and 2 bees inoculated with G. apicola and reject colony 1 and 2 bees inoculated with L. quercina. (A), (C), and (F) are depicted as nonmetric multidimensional scaling plots depicting BC dissimilarity between samples. Bold letters in legends denote statistical significance, P < 0.05 [permutation MANOVA (A, C, and F) or Pearson’s chi-square (B, D, E, and G) followed by FDR pairwise contrasts]. *P < 0.05 (Pearson’s chi-square).

Next, we investigated whether shared microbiomes could induce similar CHC profiles in unrelated individuals by inoculating newly eclosed bees from two different colonies with G. apicola or L. quercina and by analyzing their CHC profiles and behavior using in-lab behavioral acceptance assays. We found that both treatment and source colony had a significant effect on the CHC profiles of these bees (Fig. 3F, two-way permutation MANOVA using BC, treatment: F1,31 = 5.18, R2 = 0.13, P = 0.010; colony: F1,31 = 4.77, R2 = 0.12, P = 0.015; and treatment*colony: F1,31 = 0.71, R2 = 0.02, P = 0.485). However, behavioral studies revealed that bees inoculated with G. apicola accepted other bees inoculated with G. apicola and rejected bees inoculated with L. quercina, independent of genetic relatedness (Fig. 3G, Pearson’s chi-square, χ2 = 13.17, P = 0.003). Together, these data indicate that some members of the gut microbial community are sufficient to drive similar nestmate recognition cues between unrelated bees, further supporting the role of the gut microbiome in the development of chemical signatures that signify group membership in honey bee colonies.

Previous studies have indicated that the taxonomic diversity of the honey bee gut microbiota is low (32, 33), and therefore, community-level variations may not suffice for the development of unique nestmate recognition cues across independent colonies. Therefore, we hypothesized that functional genetic variations between populations of individual microbial taxa could provide additional opportunities for microbial-dependent cue divergence across colonies. Previous studies suggest that G. apicola strains vary in their ability to metabolize various carbohydrates (17, 18). Because microbial carbohydrate metabolism typically yields short-chain fatty acids and other molecules, which can be used by the host to synthesize CHCs, it is possible that genetic variations across strains of this common honey bee gut microbe may lead to the development of different host CHC profiles (10). To test this hypothesis, we first assessed the population-level genetic diversity of 64 bacterial strains of G. apicola in foragers from four different honey bee colonies by sequencing the encoding DNA sequence for elongation factor Tu (tuf), which has previously been used to assess strain diversity in other honey bee–associated taxa (32). Through pairwise sequence alignments between all strains, we found that specific strains of G. apicola are more genetically similar when present in the same bee than across different bees from the same colony (Fig. 4A, Kruskall-Wallis, H = 72.802, P < 0.001). Furthermore, bacterial strains were significantly more distantly related to each other across forager bees from different colonies than across sister foragers from the same colony (Fig. 4A, Kruskall-Wallis Dunn’s test pairwise contrasts, P = 0.003). Together, these data indicate that genetic relatedness of different strains of a resident bacterial phylotype is associated with the social relatedness of their host bees. In addition, when groups of newly eclosed bees from a single source colony were inoculated with one of four of the most genetically distant strains of G. apicola, they developed significantly different CHC profiles (Fig. 4B, permutation MANOVA using BC, F3,28 = 7.74, R2 = 0.48, P < 0.001). Likewise, in behavioral assays, inoculated bees were able to distinguish sister bees inoculated with a different strain from those inoculated with the same strain (Fig. 4C, Pearson’s chi-square, χ2 = 7.39, P = 0.007). Together, these data indicate that strain-level genetic diversity of at least one resident bacterial phylotype in the honey bee gut is sufficient to induce differences in the recognition cues of their hosts, suggesting that symbiont genetics can drive variations in complex social phenotypes of their animal hosts (8).


Fig. 4 Strain level diversity is associated with differences in CHC profile and recognition in honey bees.

(A) G. apicola coding gene tuf has more single-nucleotide polymorphisms (SNPs) between forager bees from different colonies than between those from the same colony. (B) Sister bees inoculated with different strains of G. apicola develop different CHC profiles. (C) Bees inoculated with G. apicola strain 2 or 3 can discriminate sister bees inoculated with the same strain from those inoculated with a different strain. Statistics in (A) using Kruskall-Wallis followed by Dunn’s test pairwise contrasts are shown as violin plots with points representing median and lines represent interquartile region. Statistics in (B) using permutation MANOVA followed by FDR pairwise contrasts are shown as nonmetric multidimensional scaling plots depicting BC dissimilarity between samples. Statistics in (C) using Pearson’s chi-square. Bold letters in the graph (A) and legends (B), and asterisks (C) denote post hoc statistical significance (P < 0.05).

Overall, the data we present here indicate that in the honey bee, nestmate recognition cues are defined by colony-specific gut microbial communities. Several mechanisms have been proposed for how this might occur (10). One option is that the microbiome directly produces the cue used by their insect host. While we cannot exclude this possibility, previous studies have suggested that honey bee nestmate recognition relies on colony-specific CHC profiles (3, 4), which are synthesized by the host’s subcuticular pheromone-producing oenocytes (4, 19). Therefore, because gut microbes likely do not have direct access to their host’s oenocytes, it is more likely that in the honey bee, the microbiome influences quantitative and qualitative aspects of individual CHC profiles by supplying or depleting precursor metabolites and/or by modulating the expression or activity of host genes involved in CHC synthesis, half-life, and transport (10). Future functional studies may seek to test these possible mechanisms.

Our finding that social insects can evolve to use the genetics of symbiotic bacteria as a robust and reliable proxy for colony membership supports the concept of a holobiont, whereby the genomes of the host and its microbiota interact (8). In combination with similar findings in other social hymenopteran species (34), these data suggest that the effect of the microbiome on recognition cues may have a common ancestral origin. While the selective pressures that led to the symbiotic relationship between social Hymenoptera and their microbiomes are unknown, one possibility is that this mechanism arose via a mutualism. For example, while nestmate recognition is clearly important for the social insect host, it is also likely to be beneficial for gut microbes, since host nestmate recognition prevents the invasion of unrelated animals that may host and transmit genetically distant strains of gut bacteria (3, 15). Thus, by driving host nestmate recognition, some gut microbes may reduce possible competition with unrelated strains.

While the majority of previous studies have focused on the role of the microbiome in honey bee health and immunity (14, 3537), our studies suggest that the gut microbiome also plays a fundamental role in regulating social behavioral traits of this economically important insect species. Consequently, as more data from diverse animal groups are becoming available, the effect of microbiomes on host behavior is emerging to be a fundamental rule of life (9, 12, 13), indicating that some aspects of animal behavior in general, and sociality in particular, may have evolved via an obligatory codependency between animal hosts and their symbiotic microbes (12, 13, 38).


Animal husbandry, treatments, and collections

Honey bee (Apis mellifera) colonies were reared and managed using standard beekeeping techniques across two locations near St. Louis, MO: Tyson Research Center and a residential home. To decrease relatedness between colonies used within an individual experiment, honey bee queens were sourced from different locations in the United States, including Georgia, California (Isabees), and New York (Betterbee). For collection of natural foragers from typical honey bee colonies, foragers were identified by pollen loads on their hind legs or having a distended abdomen due to nectar loads and were collected. For all experiments that included in-lab treatments or cross-fostering, capped brood frames were taken from a colony and placed in a 32°C incubator at 75% relative humidity. For cross-fostering experiments (as in Fig. 1E), ~1000 newly eclosed bees from two independent source colonies were marked with a spot of paint (Testors, Vernon Hills, IL, USA) on their thorax, and then half of each group were randomly reintroduced into either their source or an unrelated foster colony. Marked bees were then collected as returning foragers at 18 days after reintroduction. All bees used for chemical analyses were placed in individual 1.7-ml microtubes and flash frozen. All bees used in 16S sequencing analyses were washed once with 12.5% bleach in water and twice with double deionized water (39) and flash frozen. All samples were stored at −80°C until further analysis.

In-lab–treated bees were kept in a 32°C incubator at 75% relative humidity in groups of 50 in Plexiglas boxes (10 × 10 × 7 cm) (40) with a sterilized pollen patty and a hanging, inverted, sterile 1.7-ml microtube containing sterilized 25% (w/v) sugar (sucrose) water with the specific treatment and was replaced daily. To test microbial survival in 25% sugar water, we placed 50 μl of G. apicola culture [OD600 (optical density, 600 nm), ~1] in 1.5 ml of 25% sugar water overnight and then plated it on standard LB plates to ensure growth. In addition, we incubated several aliquots of the full microbial community in 25% sugar water solution for 24 hours and then performed 16S rRNA sequencing (table S7). For antibiotic treatments, a mixture of three antibiotics known to perturb insect microbiota was used (41). Antibiotic stock solutions (1000×) were composed of tetracycline (50 mg/ml) in water, rifampicin (200 mg/ml) in dimethyl sulfoxide, and streptomycin (100 mg/ml) in water. We then added 50 μl of each stock solution to 50 ml of 25% sugar water, and 1.5 ml of this working solution was added to a new inverted microtube for each treatment box every treatment day. For Fig. 2 (A to C), newly emerged bees were placed in their source colony for 3 days to establish a microbiome and were then recollected and placed in treatment boxes where they received either 25% sugar water or antibiotic treatment for 15 days. For Fig. 2 (G and H), forager bees were collected from their source colony and were placed in treatment boxes, where they received either 25% sugar water or antibiotic treatment for 3 days. Groups of 10 of these bees were then transferred to new treatment boxes with 40 newly eclosed bees and were fed 25% sugar water for 16 days. For treatments with live inoculum (Fig. 2, D and E), six forager honey bee mid- and hindguts from a single colony were dissected under sterile conditions (39), homogenized in 1 ml of sterile 25% sugar water, and centrifuged at 2700 revolutions per minute for 1 min, and 250 μl of the supernatant was added to 1.3 ml of 25% sugar water in a new inverted microtube for each treatment box. For heat-killed treatments, the remainder of the supernatant was heated at 95°C for 20 min and chilled on ice, and 250 μl was added to 1.3 ml of 25% sugar water in a new inverted microtube for each treatment box. This was repeated every treatment day (17 days). For treatments with specific cultures (Figs. 3 and 4), single colonies of bacteria were cultured in standard LB overnight to an OD of ~1 and placed in the refrigerator. Every day for the treatment period (~15 to 20 days), 50 μl of this culture was added to 1.5 ml of 25% sugar water in a new inverted microtube for each treatment box. In all cases, dead bees were removed daily from the treatment boxes. Bees were kept in treatment boxes until 14 to 21 days old, depending on the experiment and survival rate, when they were used for behavioral assays or were washed once with 12.5% bleach in water and twice with double deionized water (39), flash frozen, and kept at −80°C until further analysis.

In-lab behavioral acceptance assay

To assess recognition behaviors, we performed a modified version of the intruder assay (42, 43). In short, groups of bees, from a single colony (Fig. 3, D and E) or from two colonies (Fig. 3G), were inoculated with either G. apicola or L. quercina for 16 days following the protocol described above. After the inoculation period, bees from one of each type of treatment group were placed in groups of three in a petri dish in a humidified 32°C incubator overnight. On the day of the trial, these petri dishes were moved to a humidified 25°C room 1 hour before behavioral assays. During this time, bees from the other treatment groups (those not used to populate petri dishes) were marked with a spot of paint and placed in individual 15-ml tubes until the behavioral assay. These bees served as “intruders,” and each one was introduced into one petri dish with bees inoculated with either the same phylotype or a different phylotype. Behavioral interactions were videotaped for 10 min, and videos were subsequently scored by a blinded researcher. Intruder bees were scored as “rejected” or “accepted,” where they were considered rejected if they were bit stung and/or dragged by at least one bee, and accepted if they never received aggression from the other bees.

CHC extractions and GC analysis

CHCs were extracted from whole bees by placing individual bees into 6-ml glass vials fitted with 16-mm polytetrafluoroethylene (PTFE)/silica septa screw caps (Agilent CrossLab, Santa Clara, CA, USA). Bee CHCs were extracted in 500 μl of hexane containing octadecane (C18; 10 ng/μl) and hexacosane (C26; 10 ng/μl) (MilliporeSigma, St. Louis, MO), which served as injection standards. For extraction, each vial was gently vortexed (Fisher Scientific, Waltham, MA, USA) for 2 min at minimum speed. Extracts were immediately transferred to new 2-ml glass vials fitted with 9-mm PTFE lined caps (Agilent CrossLab, Santa Clara, CA, USA). In cases where experiments involved forager honey bees, all bees (including nonforagers) had their hind legs removed before extraction to ensure removal of pollen. One hundred microliters of each extract was transferred to a new 2-ml glass vial and stored at −20°C for further analysis; the remaining 400 μl was stored at −80°C as backup.

Representative pooled samples of foragers and nurses of known age were first analyzed by combined gas chromatography/mass spectrometry (GC/MS) for compound identification. Samples were run from 150° (3-min hold) to 300° at 5°/min. Compounds were identified by their fragmentation pattern as compared with synthetic compounds. For profile characterizations of individual bees, samples were analyzed using an Agilent 7890A gas chromatograph system with a flame ionization detector (GC/FID) and programmable temperature vaporizer injector (cool-on-column mode) and outfitted with a DB-1 20-m × 0.18-mm Agilent 121-1022 fused silica capillary column (Agilent Technologies Inc., Santa Clara, CA, USA). Sample volumes of 1.0 μl were injected onto the column. Helium was the carrier gas and applied at a constant flow rate of 1 ml/min. Analysis of the extract was carried out with a column temperature profile that began at 50°C (held for 1 min) and was ramped at 36.6°C/min to 150°C and then at 5°C/min to 280°C, where it was held for 10 min. The injector and FID temperatures were programmed to 280° and 300°C, respectively. Agilent OpenLab CDS (EZChrom Edition) software was used to calculate the retention time and total area of each peak. Data were normalized to known quantity (ng) of internal standard hexacosane.

Gut microbiome DNA extraction, 16S rRNA sequencing, and analysis

Frozen honey bee guts were dissected under sterile conditions on dry ice. Individual honey bee mid- and hindguts were homogenized by maceration with a disposable mixer sterile pestle (VWR Products). The homogenate was added to a PowerSoil Bead Solution tube (MO BIO), and DNA was extracted using a DNeasy PowerSoil DNA isolation kit (MO BIO, Carlsbad, CA) following the manufacturer’s instructions. The hypervariable V4 region of the 16S rRNA gene was amplified by polymerase chain reaction (PCR) in triplicates using primers and barcodes designed in (44). Before sequencing, the PCR products were visualized on 1.5% agarose gels, and samples that did not have a band were removed from further analysis. Note that these samples were primarily those that received antibiotic treatment or were raised with older bees that had received antibiotic treatments. All samples were pooled on the basis of concentrations, and each pool was sequenced on a single lane on an Illumina MiSeq with 2 × 250–base pair paired-end reads. The samples for this paper were split among three sequencing pools: #1, Fig. 1 (A, C, and D); #2, Figs. 1E and 2D; #3, Figs. 2F and 3B and fig. S1.

After obtaining sample sequences, sequences were demultiplexed using QIIME 2, and paired-end reads were truncated at the first base with a quality score of <Q3 using DADA2. Paired-end reads were then merged, and “amplicon sequence variants” (ASV) were identified using DADA2. Chimeric ASVs were removed, and the remaining ASVs were taxonomically classified using the Greengenes Database (45). Sample data were subsequently rarified to the lowest reasonable sample read count. Samples with a read count lower than the rarified amount were removed from the dataset. ASVs that were found in at least five samples per dataset were further taxonomically classified by aligning a representative ASV sequence to the genomes of honey bee–associated microbes via the National Center for Biotechnology Information Basic Local Alignment Search Tool (NCBI BLAST).

For sequencing pool #1 (relevant samples found at NCBI accession number PRJNA630281), 1,723,236 (861,618 pairs) total sequence reads were obtained; 798,864 pairs were filtered and denoised; 710,716 pairs were merged and identified as nonchimeric (82.4%); 233 ASVs were identified; and data were rarified to 11,675 reads per sample. For sequencing pool #2 (relevant samples found at NCBI accession number PRJNA630292), 1,429,358 (714,679 pairs) total sequence reads were obtained; 692,814 pairs were filtered and denoised; 672,361 pairs were merged and identified as nonchimeric (94.1%); 633 ASVs were identified; and data were rarified to 1794 reads per sample. For sequencing pool #3 (relevant samples found at NCBI accession number PRJNA630294), 1,280,076 (640,038 pairs) total sequence reads were obtained; 606,799 pairs were filtered and denoised; 582,712 pairs were merged and identified as nonchimeric (91%); 183 ASVs were identified; and data were rarified to 1126 reads per sample.

G. apicola strain diversity measurement

G. apicola colonies were cultured as described above and included four colonies per individual honey bee, and four bees per colony across four colonies (a total of 64 cultures). DNA was extracted from these cultures using a Qiagen DNeasy Blood and Tissue kit (MO BIO, Carlsbad, CA), and a PCR was performed targeting the elongation factor Tu (tuf) gene using primers CGATACACCAACTCGTCACT and AACAACACCAGCACCAACAG. PCR products were subsequently run on a 1.5% agarose gel, the amplicon band was excised, and DNA was extracted from the gel piece using a Wizard SV Gel and PCR Clean-Up System (Promega, Madison, WI) following the manufacturer’s directions with an extended (3 min) elution step using 60°C heated nuclease-free water, which was repeated by reusing the first elution flowthrough. The extracted DNA was then sent for Sanger sequencing by GENEWIZ (, using the forward and reverse tuf primers separately to obtain 2× coverage. Sequencing was repeated for samples that had low-quality scores (<30), and if the repeated sequence still had low-quality scores, then the sample was removed from the dataset. The resulting sequences were merged using USEARCH v10.0. Individual sample sequences were globally aligned in a pairwise manner across all samples using USEARCH v10.0, and the proportion of sequence overlap for each comparison was calculated and used in downstream analysis.

Statistical analysis

All CHC analyses included a set of 19 peaks that represent well-established honey bee CHCs, identified by comparing GC traces to published data (46). For the comparisons of CHC profiles, the relative proportion of each compound in each sample was calculated. These data were subsequently rescaled within each CHC compound to a value between 0 and 1 using the calculation 𝑧𝑖 = (𝑥𝑖 − min(𝑥))/(max(𝑥) − min(𝑥)), where 𝑥 = (𝑥1,...,𝑥𝑛) and 𝑧𝑖 is the rescaled value to limit the influence of highly expressed CHCs to the overall CHC profile. These rescaled values were used in further statistical analysis. For comparisons of gut microbiome data, ASV counts were used. For each dataset, a permutation MANOVA was run using the adonis function in the vegan package of R with BC dissimilarity measures (47). Pairwise comparisons with false discovery rate (FDR) P value correction were subsequently run on experiments where more than two groups were compared. Data were visualized using nonmetric multidimensional scaling (metaMDS function in the vegan package of R) (47) using BC dissimilarity for CHC data or BC dissimilarity or weighted UniFrac (48, 49), depending on which visualization method best depicts the statistical comparisons, for microbiome data. To identify individual ASVs that differed in abundance between groups, we performed a nonparametric Kruskall-Wallis ANOVA test followed by a pairwise Dunn’s test with FDR adjustment (when comparing three or more groups), or Mann-Whitney U test (when comparing two groups), using read counts for each ASV. For ASVs that differed between groups, a representative sequence was aligned to the genomes of honey bee–associated microbes via NCBI BLAST, and the “likely microbe” represented by each ASV was noted. For behavioral data, the proportion of intruder bees accepted by the group of bees was calculated, and a Pearson’s chi-square was run with subsequent pairwise comparisons. For G. apicola alignment data, comparisons using an ANOVA were made using the proportion of overlap between samples from the same bee, between samples from the same honey bee colony, and between samples from different honey bee colonies. These data were visualized using violin plots (50).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank Tyson Research Center and E. Hartz for beekeeping assistance; A. Cantu and B. Henson for bee collection assistance; B. Wang for sequencing library preparation support; J. Hoisington-López, M. Crosby, E. Martin, and B. Koebbe for next-generation sequencing and high-throughput computing support at the Edison Family Center for Genome Sciences and Systems Biology at Washington University in St Louis School of Medicine; and G. E. Robinson and the Robinson Lab members for comments on the manuscript. Funding: NSF 1545778, 1707221, and 1754264 to Y.B.-S.; NIH R01AT009741 to G.D.; NSF GRFP DGE-1745038 to C.L.V.; and CIHR, CRC, NSERC, and CIFAR to J.L. Author contributions: Conceptualization: C.L.V., G.D., and Y.B.-S. Methodology: C.L.V., B.A.-O., J.J.K., J.L., G.D., and Y.B.-S. Investigation: C.L.V., I.M.C., and J.J.K. Resources: J.L., G.D., and Y.B.-S. Supervision: G.D. and Y.B.-S. Writing, original draft: C.L.V. and Y.B-S. Writing, review and editing: C.L.V., I.M.C., B.A.-O., J.J.K., J.L., G.D., and Y.B-S. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. All sequencing reads are available on the NCBI sequence read archive, accession numbers PRJNA630281, PRJNA630292, and PRJNA630294.

Stay Connected to Science Advances

Navigate This Article