Gut microbiota contribute to variations in honey bee foraging intensity

Abstract Gut microbiomes are increasingly recognized for mediating diverse biological aspects of their hosts, including complex behavioral phenotypes. Although many studies have reported that experimental disruptions to the gut microbial community result in atypical host behavior, studies that address how gut microbes contribute to adaptive behavioral trait variation are rare. Eusocial insects represent a powerful model to test this, because of their simple gut microbiota and complex division of labor characterized by colony-level variation in behavioral phenotypes. Although previous studies report correlational differences in gut microbial community associated with division of labor, here, we provide evidence that gut microbes play a causal role in defining differences in foraging behavior between European honey bees (Apis mellifera). We found that gut microbial community structure differed between hive-based nurse bees and bees that leave the hive to forage for floral resources. These differences were associated with variation in the abundance of individual microbes, including Bifidobacterium asteroides, Bombilactobacillus mellis, and Lactobacillus melliventris. Manipulations of colony demography and individual foraging experience suggested that differences in gut microbial community composition were associated with task experience. Moreover, single-microbe inoculations with B. asteroides, B. mellis, and L. melliventris caused effects on foraging intensity. These results demonstrate that gut microbes contribute to division of labor in a social insect, and support a role of gut microbes in modulating host behavioral trait variation.


*
replicates indicated that single microbe inoculated bee guts were mostly composed of the intended honey bee-associated treatment bacteria, while microbiota-depleted bee guts were mostly composed of chloroplast and mitochondrial reads, as well as non-honey bee associated and unassigned bacteria.(B-C) Quantitative PCR analysis indicated that single microbe inoculated bee guts had a higher bacterial load than microbiota-depleted bees, measured as absolute abundance of 16S rRNA gene copies normalized to Actin gene copies for each sample.used to construct the colony; half with a paint mark on the thorax, and half with a plastic tag attached to the thorax (~3 mm diameter, ~1 mm thick; "big-back" bees).Four days later, ~400 newly eclosed bees from ~10 mixed source colonies (headed by naturally mated queens) were introduced to increase the proportion of precocious foragers in the first cohort [30,43,45].The entrance to the hive was blocked by a piece of Plexiglas with holes in it to prevent the big-back bees from leaving the hive, but to allow paint-marked bees to come and go freely.We observed that guard bees patrolled the hive entrance on the inside of the Plexiglas piece.
Bees were collected at 10 days of age: returning active foragers and nurses collected as described above, and big-back/inactive forager bees collected as they were attempting to leave the hive via the holes in the plastic, indicating they too were in the behavioral state of foraging but had never left the hive.10 bees per behavioral task group (nurse, active forager, inactive forager) were collected from all three colonies, with the exception that only seven nurses from the focal group in colony 2 were found and collected.ASVs that were taxonomically identified as a bee specific genus by the BEExact database, but were unclassified at the species level, were subsequently classified to species level if possible, using NCBI megaBLAST.To do this, representative 16S rRNA gene V4 sequences for each ASV obtained from the QIIME2 rep-seqs.qzaobject were BLASTed against the entire NCBI nucleotide database, and unclassified ASVs were called as the first identified full species with the lowest Escore, if query cover > 80% and percent identity > 92%.Prior to statistical analyses, ASVs that were identified as mitochondrial or chloroplast were removed from the data.Mitochondrial and chloroplast reads were retained for visualization of single-microbe inoculated and microbiotadepleted bee gut microbial communities.
To estimate the abundance of individual honey bee-associated microbial species in each sample, the read counts for all ASVs that matched the same species were combined.We retained the following honey bee-associated species:  [20,98,99].We also retained other unclassified Lactobacillus species, Limosilactobacillus spp., Arsenophonus spp., Klebsiella spp., and Mixta spp.because they made up a significant portion of reads (>500 reads in at least 1 sample).In the case of single-microbe inoculated and microbiota-depleted bees (see below), Enterococcus spp., Staphylococcus spp., and Unassigned (not able to be classified below Kingdom of Bacteria) were additionally retained for visualization.All other ASVs, each comprising 0-1% of total reads per sample, and representing microbes that were not classified in the aforementioned groups and/or are not typically associated with honey bees, were combined and labeled as "other."The variation between samples in the abundance of "other" microbes is similar to previously published datasets [22,24].To calculate the proportion of each species in each sample, that species' read count was divided by the total read count for each sample.Raw read counts and/or proportion data were used in analyses as described below.For the depiction of single-microbe inoculated bee gut microbial communities, "other" ASVs were combined based on genus (or the lowest taxonomic level identified), and those representing >500 reads in at least 1 sample are depicted in a stacked barplot, whereas all genera are included in supplementary data on Mendeley Data (DOI: 10.17632/f2s47y3nhn).
Due to the compositional nature of 16S rRNA gene amplicon sequencing data, typical measures of beta diversity and differential abundance analyses have limitations [48,50,100].To estimate absolute bacterial species abundances in individual samples ("absolute abundance"), we quantified the bacterial load in each sample using quantitative PCR (qPCR), as in [22,23].We performed standard curves using serial dilutions of plasmids (TOPO pCR2.1 (Invitrogen)) containing the target sequence (10 8 -10 3 copies per 3 µl), which were calculated from the molecular weight of the plasmid and the DNA concentration of the plasmid.Primer efficiencies were measured using: E = 10 (-1/slope) [101], and the copy number of each target in 3 µl of DNA sample was calculated from the sample's Ct score, primer efficiency, and standard curve using: n = E (intercept-Ct) x(DNA extraction elution volume/3) [22,23].These values were calculated for both the 16S rRNA gene and the Actin gene.To account for differences in DNA extraction efficiency, we calculated a normalized number of 16S rRNA gene copies by dividing the number of 16S rRNA gene copies by that sample's number of Actin gene copies and multiplying this by the median number of Actin copies for that experiment [22,23].To calculate the absolute abundance of each microbial species in each sample, the relative abundance (proportion) of each species in each sample was multiplied by the normalized number of 16S rRNA gene copies in that sample [22,23].We then took the log10 value of each of these numbers and used these in further analyses, replacing 0s with 1s before we did so.
Standard curves using serial dilutions of plasmids (TOPO pCR2.1 (Invitrogen)) containing the target sequence (10 8 -10 3 copies per 3 µl) were calculated from the molecular weight of the plasmid and the DNA concentration of the plasmid.
We found an effect of colony replicate on gut microbial community structure across all experiments, indicating that bees from different colonies have different gut microbial communities and that early environment and/or genetics influences the composition of the gut microbial community.This was the case when comparing bees from different typical colonies (Supplementary Fig. 3A-B), among SCC bees originating from different colony sources (Supplementary Fig. 3C-D), and across bees from different big-back colonies (Supplementary Fig. 3E).Colony differences in honey bee gut microbial community composition have been previously reported [44].

Single-microbe inoculations
To produce single-microbe inoculated and microbiota-depleted bees, modified methods from [38,52] were used.Specifically, tan-colored pupae with dark eyes were gently removed from brood frames from three or four source colonies per trial using ethanol sterilized forceps, placed dorsal side down in sterilized 3D-printed dental grade resin modular pupation plates (Supplementary asteroides) or MRS + 20g/L fructose (B.mellis and L. melliventris) media plates and were grown under anoxic conditions at 35°C for 3 nights (Oxoid Anaerojar with Anaerogen bags, Thermo Scientific, Waltham, MA, USA).A single colony of bacteria from each plate was cultured under anoxic conditions (Anaerobic Hungate Culture Tubes with air displaced by CO2, VWR, Radnor, PA, USA), shaking at 35°C for 3 nights and was then retained at 4°C.Every day, starting 2 days prior to the onset of inoculations, new cultures from this same original stock were prepared by culturing 50 µl of the stock in 5 mL of fresh MRS/MRS + 20g/L fructose broth under anoxic conditions, shaking at 35°C for 2 nights.These cultures were spun down (3000 rpm for 5 min) and resuspended in 1x PBS to an OD of 1 and 50 µl of a new, fresh culture solution was added to 1.7 mL sterile 25% sucrose water in a new inverted microtube for each treatment box every morning over the course of a 5-day inoculation period.For microbiota-depleted bees, 50 µl of sterile 1x PBS was added to 1.7 mL 25% sucrose water in a new inverted microtube for each treatment box every morning over the course of a 5-day inoculation period.Bees were then fed sterile 50% sucrose water for 2 days following the 5-day inoculation period.Very few bees (1-2 per box) died during this inoculation period.Bees were kept in treatment boxes until 7 days old.
Then they were either used in behavioral assays (described below) or flash frozen for 16S rRNA gene amplicon sequencing to confirm single-microbe inoculation and microbiota-depleted status (Supplementary Figure 2A-C).Unfortunately, we did not keep track of the treatment box each bee came from, rather bees for each treatment were pooled and randomly chosen for behavioral assays or 16S rRNA gene amplicon sequencing.Microbe survival in 25% sucrose water was ensured by placing 50 µl of each microbial culture (OD ~1) in 1.5 mL 25% sucrose water overnight and then plating this solution using species specific culturing conditions.

Foraging assays and analysis using barcoded bees
After single-microbe inoculations, 7-day old inoculated and microbiota-depleted bees were coldanaesthetized, a unique barcode was affixed to their thorax using super glue, and a spot of colored paint associated with their treatment group was placed on their abdomen.An equal number (~100) of bees from one single-microbe inoculated group and a corresponding microbiota-depleted group (same mix of colony backgrounds, same age) were given barcodes from different sets, and were placed together in an experimental double-cohort colony (DCC) in order to assess the relative effects of the two different inoculation treatments on behavioral maturation rate and foraging intensity in a common colony environment.DCCs were established similarly to SCCs with a new, mated queen, 1 plastic honeycomb frame provisioned with honey and pollen paste, 1 empty plastic honeycomb frame, ~200 7-day old barcoded treatment bees (100 of bees from one single-microbe inoculated treatment group and 100 from a corresponding microbiota-depleted group) and ~800 newly eclosed unbarcoded "background" bees (mixed colony backgrounds), whose addition contributed to accelerated behavioral maturation of the older cohort [43,45].In the presence of the older cohort, younger "background" bees are not expected to contribute to the foraging force [104] and did not throughout the course of the experiment.Four colony replicates per single-microbe treatment were performed (B.asteroides 1-4, B. mellis 1-4, L. melliventris 1-4), for a total of 12 experimental colonies.Two experimental days, corresponding with B. mellis colony replicate 3 days 1-2 and L. melliventris colony replicate 2 days 1-2, experienced bad weather, and thus no foraging occurred on those days.
Experimental DCCs were kept in a dark, temperature (32°C) and humidity (~50%) controlled building, with access to the outside environment through a plastic tube.To track flight activity, an entrance monitor with a video camera (described below) was attached to the hive entrance, on the outside-facing end of this tube (Supplementary Fig. 5).The camera recorded videos of bees entering and leaving the hive from 05:00 to 21:00 daily for a total of six days.
Barcodes in these videos were detected as in [54], and detections were filtered to remove tracking errors and misidentifications.A mean of 97.82%, with a range between 96.24%-99.02%, of detections were retained after these filtering steps.We then applied a flight activity detector [56,57] (described below) to the remaining barcode detections to identify passes through the entrance monitor, which were used to identify foraging trips.Computer code for this detector can be found at: https://github.com/gernat/btools.
Passes through the entrance monitor determined by the flight activity detector were filtered to remove unused barcodes, whose detections either occurred through detection error or due to bees accidentally entering the wrong experimental colony.A mean of 99.19%, with a range of 95.4%-100%, of detections were retained after this filtering.Passes were then classified as "incoming," "outgoing," or "other."Because incoming passes had a lower error rate than outgoing passes (Supplementary Table 8), incoming passes alone were used to denote a foraging trip.
Individual foraging trips that occurred within 5 min of each other were condensed into a single trip.A bee's first day of foraging was defined as the first day on which it performed at least 4 foraging trips, with at least 50% of these trips occurring during peak foraging hours (11:00-15:00 CST), as per previous studies [35,58].These criteria yield similar automated thresholds corresponding to human observations of foraging behaviors [35].All incoming passes from this first day of foraging were counted as foraging trips, as were all subsequent incoming passes [35].
Likewise, an individual bee was considered a "forager" from this day forward.
To compare behavioral maturation rate between treatment groups (single-microbe

Entrance monitor
Flight activity was recorded with an improved version of the entrance monitor described in [64].
Improvements aimed to make the entrance monitor easier to traverse by the bees, increase bCode detection rate and recording frame rate, and make the entrance monitor enclosure cheaper and easier to manufacture.The changes we made to achieve these goals are described below.

Enclosure
When passing through the entrance monitor enclosure, bees need to traverse a maze.This maze slows them down so they can be recorded multiple times while exiting or returning to the hive.
The longer dimension of this maze was shortened to approximately half its original size.In addition, we removed two of the three inner walls.These changes helped the bees navigate the maze and made it less likely that they considered it part of their hive and congregated in it.The roof of the maze was changed to glass with an antireflective coating.This coating eliminated reflections of the camera and other enclosure components on the glass and thus helped to increase the barcode detection rate.
To manufacture the enclosure, we created a three-dimensional model of it (Supplementary Fig. 5) in Onshape (PTC Inc., Boston, Massachusetts, USA), an online computer-aided design software system.This model was printed on a Form 2 printer (Formlabs, Somerville, Massachusetts, USA), using Clear Resin (Formlabs, Somerville, Massachusetts, USA).Clear Resin cures to optical translucency, which reduced motion blur by permitting more light to pass through the enclosure walls than the opaque material we used before.After printing, the enclosure was cleaned in a Form Wash (Formlabs, Somerville, Massachusetts, USA) and cured in a Form Cure (Formlabs, Somerville, Massachusetts, USA), using manufacturer-recommended settings.

Camera
The entrance monitor camera was upgraded to a Raspberry Pi camera module v2.1 (Raspberry Pi Ltd, Cambridge, UK).The higher resolution of this camera made it possible to use pixel binning to record the maze at 1640 x 1232 pixels.This improved the camera's performance under lowlight conditions and led to bigger barcodes in the recorded footage, which contributed to the higher barcode detection rate.For data storage, we switched from capturing images in burst mode to recording video.This enabled the camera to automatically adjust to changes in illumination.It also reduced the amount of data and made data storage more time-efficient.These changes allowed us to step up the recording frame rate to 10 Hz, a five-fold increase that enabled us to capture each bee multiple times despite them passing faster through the shorter, simplified maze.

Performance
To obtain an estimate of the entrance monitor identification rate, we manually annotated all bees that were not automatically identified in 17,125 images that were randomly sampled from entrance monitor videos recorded during an unrelated experiment.Of the bees with a barcode that was at consisting of 70% and 30% of the trajectories, respectively.The training set was subsampled to ensure that all three classes were represented equally, which reduced its size to 62% of the ground truth trajectories.

Pass detection and classification
Our flight activity detector first groups a bee's barcode detections that are at most c=30 s apart into a pass.The cutoff c corresponds to the 99.9th percentile of the time between barcode detections in the training set, but was otherwise chosen arbitrarily.For each pass, the detector then computes the features listed in Supplementary Table 6 on the subset of successive barcode detections.Feature calculations were limited to this subset to ensure that pass features are not calculated across barcode detection gaps.Next, the detector predicts the pass class (i.e., whether it is an incoming, outgoing, or "other" pass), using a random forest consisting of 50 decision trees.
This random forest was trained with default parameters on the training set trajectories, using the R package randomForest [66].

Performance
If the flight activity detector split an annotated pass into multiple detected passes, we counted only one of the detections as a positive.The remaining detections were considered to be false positives and thus decreased the detector's positive predictive value.Similarly, if a detected pass spanned multiple annotated passes, only one of these passes was considered to be detected.
The other passes were treated as false negatives and thus decreased the detector's sensitivity.
Finally, since the purpose of the flight activity detector is to identify hive exits and returns, its performance (Supplementary Table 7) on the "other" class was not evaluated and detections classified as "other" were ignored when calculating the performance on the incoming and outgoing class.
Therefore we followed analyses outlined in[48].In short, for beta-diversity analyses, we used clrtransformations of raw read counts for each sample in statistical tests and Aitchison distance for visualization (see "Statistical analysis" section below).To estimate the differential relative abundance of individual microbial species between samples through tools that have a compositional foundation (referred to as "relative abundance" throughout the text), we used Analysis of Composition of Microbiomes with Bias Correction (ANCOM-BC)[49][50][51] on raw read counts.This test relies on the total sample read count and relative abundance of individual microbes in each sample to estimate the true abundance of individual microbes in each gut ecosystem while accounting for the compositional nature of sequencing data[49].

Fig. 4 ,
Fig.4, design file on Mendeley Data), which were fitted into sterilized Plexiglas cages, and kept

Supplementary Figure 1. Bees from single cohort colonies differ in gut microbial community. (A-B)
Age-matched typical-age nurses and precocious foragers from a second single cohort colony did not differ in gut microbial community structure at one week of age (A), but age-matched over-age nurses and typical-age foragers significantly differed in gut microbial community structure at three weeks of age (B). 1 week: Permutation MANOVA using Aitchison Distance, F1,19 = 1.4,R 2 = 0.07, p = 0.158; n = 10 bees. 3 weeks: Permutation MANOVA using Aitchison Distance, F1,19 = 2.6, R 2 = 0.128, p = 0.001; n = 10 bees.Depicted as PCA plots.Lowercase letters in legends denote statistically significant groups.(C-D) Age-matched typicalage nurses and precocious foragers from a second SCC differed in relative abundance of one individual microbial species (C), and age-matched over-age nurses and typical-age foragers differed in relative abundance of two individual species (D).Depicted as stacked bar plots, with each bar representing a single bee's microbiome.Asterisks in legend: *, p£0.05, **, p£0.01,ANCOM-BC between nurses and foragers.See Supplementary Table 3 for all p values.(E-F) Age-matched typical-age nurses and precocious foragers from a second SCC did not differ in absolute abundance of individual microbial species or in the total normalized number of 16S rRNA gene copies (E), while age-matched over-age nurses and typical-age foragers differed in absolute abundance of four individual species and the total normalized number of 16S rRNA gene copies (F). 10 x number of 16S rRNA gene copies, calculated by multiplying the relative abundance each microbe in each sample (determined through 16S rRNA gene sequencing) by the normalized number of 16S rRNA gene copies in the sample (determined through qPCR).Depicted as dot plots with all data points plotted, line represents median, n = 10 bees.*, p£0.05, **, p£0.01,Permutation ANOVA Test between nurses and foragers.See Supplementary Table 3 for all p values.

Figure 2. Single microbe inoculated bees differ in gut microbial community
. (A) 16S rRNA gene sequencing of a subset of single microbe inoculated bee

Table 7 :
3 colonies.Depicted as PCA plots.Lowercase letters in legends denote statistically significant groups as determined by Pairwise Permutation MANOVA.16S rRNA gene amplicon sequencing reads per sample.Absolute abundance, 10 x median number of 16S rRNA gene copies per sample.This was calculated by multiplying the relative abundance of each microbe in each sample (determined through 16S rRNA gene amplicon sequencing) by the normalized number of 16S rRNA gene copies in the sample (determined through qPCR).Relative abundances were analyzed via ANCOM-BC, while absolute abundances were analyzed via Permutation ANOVA.p values were adjusted to account for multiple comparisons through FDR adjustment.NAs represent taxa that were not present in enough samples to reliably analyze.n = 10 bees/task at each age, 1 source colony.single microbe inoculation experiment, with inoculation treatment and day as main factors, and replicate and individual as random factors, followed by pairwise comparison results between inoculation treatments on each day.n = 4 colonies.Flight activity detector performance.
we also quantified the degree of skew in foraging intensity among all workers and calculated the Gini coefficient[35, 59]for each experimental colony across all days and each experimental colony on each day.The Gini coefficient values we obtained (Table2) are similar to those in previous studies[35, 59], indicating differences in foraging intensity between individuals in each colony.Finally, to determine which specific bees performed the majority of the foraging for each experimental colony on each day, we ranked individuals based on the proportion of total foraging trips they performed for the colony on each day, and defined a group of "elite foragers" as the subset of bees (variation, 5-40% of the foragers) performing >50% of the foraging trips for each colony on each day[35, 59].Due to the large size of these data sets (owing to high dimensional automated behavioral tracking) number of foragers per treatment group on each experimental day, the number of foraging trips per bee and treatment group, and the number of elite foragers per treatment group on each day are available on Mendeley Data (DOI: