Gut microbiota influences onset of foraging-related behavior but not physiological hallmarks of division of labor in honeybees

ABSTRACT Gut microbes can impact cognition and behavior, but whether they regulate the division of labor in animal societies is unknown. We addressed this question using honeybees since they exhibit division of labor between nurses and foragers and because their gut microbiota can be manipulated. Using automated behavioral tracking and controlling for co-housing effects, we show that gut microbes influence the age at which bees start expressing foraging-like behaviors in the laboratory but have no effects on the time spent in a foraging arena and number of foraging trips. Moreover, the gut microbiota did not influence hallmarks of behavioral maturation such as body weight, cuticular hydrocarbon profile, hypopharyngeal gland size, gene expression, and the proportion of bees maturing into foragers. Overall, this study shows that the honeybee gut microbiota plays a role in controlling the onset of foraging-related behavior without permanent consequences on colony-level division of labor and several physiological hallmarks of behavioral maturation. IMPORTANCE The honeybee is emerging as a model system for studying gut microbiota-host interactions. Previous studies reported gut microbiota effects on multiple worker bee phenotypes, all of which change during behavioral maturation—the transition from nursing to foraging. We tested whether the documented effects may stem from an effect of the microbiota on behavioral maturation. The gut microbiota only subtly affected maturation: it accelerated the onset of foraging without affecting the overall proportion of foragers or their average output. We also found no effect of the microbiota on host weight, cuticular hydrocarbon (CHC) profile, hypopharyngeal gland size, and the expression of behavioral maturation-related genes. These results are inconsistent with previous studies reporting effects of the gut microbiota on bee weight and CHC profile. Our experiments revealed that co-housed bees tend to converge in behavior and physiology, suggesting that spurious associations may emerge when rearing environments are not replicated sufficiently or accounted for analytically.

at specific stages in an animal's life.However, behavior can change dramatically with (st)age, and some species even exhibit transitions and reversals between distinct behavioral states.In rodent models, it was shown that gut bacteria can influence both the early canalization of behavioral development, with widespread consequences on cognitive ability later in life (10,11), and the physiological mechanisms that determine behavioral variation within social groups, such as dominance hierarchies (12).So far, few studies have attempted to map the influence of symbiotic organisms onto developmen tal axes of behavior.
Eusocial insects (ants, termites, some bees, and wasps) live in complex societies in which individuals specialize on different tasks during adult life.Morphologically distinct queen and worker "castes" are typically determined early during development, and their developmental programs cannot be reversed (13).However, adult workers can some times transition between defined physiological/behavioral states (i.e., polyethism) (14,15).Eusocial insects are therefore studied to understand how morphological, physiolog ical, and behavioral diversity can derive from the same genetic makeup (16).While developmental trajectories are known to be regulated by (epi)genetic mechanisms in response to dietary and environmental cues (17,18), individuals from different (sub-)castes often show differences in gut microbiota composition or structure (19)(20)(21)(22)(23).These differences are generally assumed to be a consequence of different host physiology or dietary preferences.However, whether the gut microbiota could in turn play a regulatory role in division of labor remains unknown (24,25).
Among eusocial insects, the honeybee has emerged as a model to address these questions (24,26,27) because (i) it has a well-characterized, simple and stable "core" gut microbiota (28), (ii) it is highly social, exhibiting behaviors that the gut microbiota may influence, and (iii) its microbiota can be manipulated by generating gut microbiotadepleted (MD) bees and inoculating them with defined communities (7,27,29).The gut microbiota of worker honeybees has been suggested to influence various host phenotypes, including aspects of neurophysiology and consequent cognitive abilities (7,(29)(30)(31), collective behavior (7), weight gain (32), and cuticular hydrocarbon (CHC) profiles (33), which are used in nestmate recognition and to indicate behavioral sub-caste (34).These phenotypes all covary with behavioral state: honeybee workers generally spend their first 2-3 weeks caring for brood inside the hive ("nursing") and performing other in-hive tasks, then undergo a rapid behavioral transition to foraging-regularly leaving the nest in search of food.This transition is regulated by hormones and is associated with profound physiological and behavioral changes, including in CHC profile, weight, gene expression, dietary preference, and gut microbiota composition (21,23,(35)(36)(37)(38).We therefore hypothesized that the detected effects of the gut microbiota on different aspects of honeybee physiology may be indirect and mediated by an effect of the gut microbiota on behavioral maturation.For example, all documented effects would be expected if the gut microbiota accelerated or retarded behavioral maturation.
To test this hypothesis, we conducted a series of experiments to assess the effect of the gut microbiota on behavioral maturation.We addressed this at the behavioral level with an automated tracking system in the laboratory, calculating the age at which bees made the first trip to a foraging arena, the proportion of time they spent in the arena, and the total number of foraging trips performed.We also measured several physiolog ical hallmarks of behavioral maturation, such as CHC profile, weight, hypopharyngeal (HP) gland size (these glands degenerate during maturation [39]), and gene expression.Overall, our results suggest an effect of the gut microbiota on the timing of the first trip to the foraging arena but not on the total foraging rate of a colony nor on any of the physiological hallmarks associated with behavioral maturation.
Therefore, the previously documented effects of the honeybee gut microbiota cannot be explained by an effect of the microbiota on the transition between nurse and forager physiology.Our results are also incongruent with previous studies which suggested that the honeybee gut microbiota modifies the host CHC profile with consequences on nestmate recognition (33) and promotes host weight gain (32).A possible explanation for these discrepancies may be that previous studies used several individuals from the same cage for statistical analyses and did not recognize that bees can precociously transition to a forager state under these experimental conditions.Individuals within a cage engage in social interactions, and hence, they are not independent from each other in aspects of behavior and physiology, and cages show skewed distributions of nurses and foragers.Treating co-housed bees as individual data points in statistical analyses can result in spurious associations between gut microbiota composition and host phenotypes (40,41), particularly, albeit not exclusively, for those phenotypes that change during behavioral maturation.

The gut microbiota accelerates the onset of foraging-like behavior under an automated behavioral tracking system
To determine whether the gut microbiota influences the rate of foraging-like behav iors, we reanalyzed behavioral tracking data from a previous study (7).This experi ment comprised nine pairs of microbiota-depleted (MD) and microbiota-colonized (CL) sub-colonies consisting of ca. 100 age-matched workers.These bees had been manually extracted from nine hives at the pupal stage and incubated under sterile conditions.The newly emerged adult bees were then inoculated (CL), or not (MD), with a gut homoge nate from five nurse bees.Each sub-colony could freely move between a nest-box (30°C, 70% relative humidity [RH] in constant darkness) and a foraging arena subject to cycles of light, temperature, and humidity mirroring the external environment (Fig. 1a; Movie S1).We observed that bees in the foraging arena frequently attempted to fly, while bees in the nest never did, suggesting that they clearly perceived the separation between the nest and foraging environments (Movie S1).The position and orientation of each bee in each sub-colony were tracked by a pair of camera systems using unique matrix barcodes (ARTag library [42]) attached to the bees' thoraces.Bees were tracked for a week, starting 3 days after adult emergence and treatment inoculation so that the gut microbiota would have been fully established (28,43).There was no significant effect of the microbiota status on the number of trips to the foraging arena (Fig. 1b; Wilcoxon matched-pairs signed-rank test: V = 29, P = 0.50) nor the proportion of time spent in the foraging arena (Fig. 1c; Wilcoxon matched-pairs signed-rank test: V = 25, P = 0.82).However, CL bees started to perform trips to the foraging arena on average 15 h earlier than MD bees (when bees were between 5 and 6 days old; Fig. 1d; paired t-test: t = −4.21,df = 8, P = 0.003).This acceleration of the average age at first foraging occurred in all nine sub-colony pairs.

The gut microbiota does not modify the CHC profiles of honeybees
The CHC profile of bees changes during the transition from nursing to foraging (38).To assess the effect of the gut microbiota on the CHC profile of bees, we randomly sampled 8-10 bees from each of the 18 sub-colonies (n = 177) at the end of the automated behavioral tracking experiment (when bees were 10 days old) for CHC analysis.Previ ously performed amplicon sequencing and qPCR analyses targeting the 16S rRNA gene from gut samples of these same bees had confirmed that CL and MD bees differed, as expected, in both gut microbiota composition and total load (see extended Fig. 1 in reference 7).However, in contrast to the previous study (33), we found no significant effect of the gut microbiota on the CHC profile [Fig.2a; Table S1; permutational multivari ate analysis of variance (PERMANOVA) using Bray-Curtis dissimilarities calculated from the centroids of each sub-colony: n = 18, F (1,17) = 0.89, R 2 = 0.04, P = 0.63].
The independence of CHC profile from microbiota status was confirmed by collecting CHC profile data from bees of our previously published RNA-sequencing experiment (7), in which we reared CL and MD bees from 10 different hives.This experiment included two additional treatments where bees were colonized with either (i) a community of 13 strains covering the predominant species of the honeybee gut microbiota (CL_13; see Table S4 in reference 7) or (ii) a single core microbiota member, Bifidobacterium asteroides (CL_Bifi).Bees from 10 different hives were reared in cages of 20 individuals in an incubator for a week after treatment inoculation (one cage per treatment per hive, except for MD bees which were produced in three cages per hive to have a surplus in case of contaminations; see reference 7 for additional details).To assess the effect of the microbiota on body and gut weight (see next section: "The gut microbiota modifies neither body and gut weight nor hypopharyngeal gland size"), we weighed 3-10 bees from 58 cages (548 bees) as well as their guts.We then randomly sampled one to three bees from each of 46 cages (at least one cage per treatment per hive) for CHC analyses (n = 120).The bees of the four different treatments differed both in gut microbiota composition and total bacterial load with the MD bees having lower loads than the other three treatment groups, the CL_Bifi bees being dominated by a Bifidobacterium phylotype, and the other two colonization treatments having more diverse communities as expected (extended Fig. 4 in reference 7).CHC analyses of the bees of these four Boxplots show the median and first and third quartiles.Whiskers show the extremal values within 1.5×, the interquartile ranges above the 75th and below the 25th percentile.**P < 0.01; NS, not significant, as calculated by paired t-test (two sided).
The camera and bee icons in panel a were created with BioRender.com.treatments confirmed our previous results (i.e., there was no significant effect of the gut microbiota on the CHC profile; Fig. 2b; PERMANOVA using Bray-Curtis dissimilarities calculated from the centroids of each cage: n = 46, F (3,45) = 1.13,R 2 = 0.07, P = 0.21).
These results differ from those of Vernier et al. (33) who concluded that the honeybee gut microbiota affects the CHC profile of bees.In our experiments, bees were sampled at two time points in a restricted time window in the life of adult worker bees (8 and 10 days of age for the RNA sequencing and automated behavioral tracking experiments, respectively).To rule out the possibility that the absence of an effect of the microbiota on the CHC profile could be specific to the two selected time points, we conducted two new experiments.We first reared CL and MD bees originating from nine hives in separate groups of 25 bees (18 cages) and tracked the development of the CHC profiles from day 1 to day 11 post-eclosion by sampling one individual per cage every 2 days, starting from the day of treatment inoculation (MD, n = 54; CL, n = 54).We next housed CL and MD bees from a single hive in 20 different cages (10 cages per treatment, which also allowed us to quantify caging effects on the CHC profiles, see below) and sampled them at days 8 (MD, n = 30; CL, n = 30) and 15 post-emergence (MD, n = 29; CL, n = 29).While the CHC profiles changed over time, there was again no significant effect of the gut microbiota on CHC profiles in either follow-up experiment (Fig. 2c and d; treatment effects, PERMANOVA using Bray-Curtis dissimilarities calculated from the centroids of each cage: time-series experiment, n = 18, F (1,17) = 1.22,R 2 = 0.06, P = 0.20; single-colony experiment, n = 20, F (1,19) = 0.85, R 2 = 0.05, P = 0.66; time effects, PERMANOVA using Bray-Curtis dissimilarities: time-series experiment, n = 90, F (1,89) = 30.30,R 2 = 0.21, P = 0.001; single-colony experiment, n = 118, F (1,117) = 12.09, R 2 = 0.08, P = 0.001).

The gut microbiota modifies neither body and gut weight nor hypopharyng eal gland size
Because foragers are lighter than nurses (37) and possess degenerated HP glands (39), we tested whether the microbiota affected these physiological hallmarks of behavioral maturation, using data collected for the RNA-sequencing experiment.There was no significant difference in fresh weight (whole body and gut only) between MD bees and any of the differently colonized bees at 8 days of age (Fig. 3a and b; linear mixed effects models fitted by restricted maximum likelihood (REML) with colony of origin and cage as nested random effects: n = 548, body weight, F (3,41) = 0.94, P = 0.43, gut weight, F (3,41) = 0.22, P = 0.88).We also measured HP gland size from a subset of the bees (n = 28) used in brain and gut RNA sequencing.There was also no significant difference in HP gland size between treatments (Fig. 3c and d; Kruskal-Wallis test: χ² = 2.75, df = 3, P = 0.43).
These findings are inconsistent with Zheng et al. (32), who found that bees inoculated with a gut homogenate exhibit greater weight gain (for both body and gut weight) than microbiota-depleted bees.However, Zheng et al. (32) reported differences in body weight (relative to initial body weight) between CL and MD bees from day 7 onward, while we had assessed the effect of the gut microbiota on weight only in 8-day-old bees.We also used sterilized pollen to feed the bees, while Zheng et al. (32) used a sterilized bee bread diet in their longitudinal experiment.Therefore, we performed two additional experiments to better match the experimental procedure of Zheng et al. (32).We reared CL and MD bees from six hives in groups of 30 (one MD and one CL cage per hive).Individuals were uniquely paint marked and weighed every 2 days for 10 days from the day of treatment inoculation.In a first experiment, bees from three hives were fed bee bread and sugar water (SW) ad libitum, while in a second experiment, bees from three hives were fed pollen instead of bee bread.The microbiota had no significant effect on weight gain in the bee bread experiment, and there was no significant interaction between time and treatment either (Fig. 3e; linear mixed effects model fitted by REML with colony of origin, cage, and bee individual as nested random effects: n = 510, time, F (4, 400) = 54.06,P < 0.0001, treatment, F (1,2) = 0.02, P = 0.91, time × treatment, F (4, 400) = 0.66, P = 0.622).However, there was a statistically significant interaction between time and treatment in the pollen experiment, with MD bees being heavier than CL bees from day 7 onward (Fig. 3f; linear mixed effects model fitted by REML with bee individual, cage, and colony of origin as nested random effects: n = 400, time, F (4, 312) = 126.98,P < 0.0001, treatment, F (1,2) = 2.10, P = 0.28, time × treatment, F (4, 312) = 2.94, P = 0.021).This result is in the opposite direction compared to the effect reported by Zheng et al. (32) who concluded that the microbiota promotes weight gain.

Gnotobiotic bees segregate into nurses and foragers with distinct physiology and behavior
While analyzing the CHC profiles of the experiments mentioned above, we realized that bees always clustered into two distinct groups independently of the treatment (Fig. 4a; Fig. S1A, B, and C).These two types of CHC profiles corresponded to the typical nurse and forager CHC profiles described in Kather et al. (38).To confirm that these CHC clusters represented nurses and foragers, we compared the CHC profiles of our 8-day-old gnotobiotic bees from the RNA-sequencing experiment to those of conventional nurses (randomly sampled within hive cells and with pollen-filled guts, n = 51) and foragers (randomly sampled on landing boards, carrying pollen, and with nectar-filled guts, n = 9) from the same 10 hives used in the RNA-sequencing experiment.The CHC profiles of these conventional nurses and foragers perfectly segregated into the two clusters (Fig. 4a).
This CHC-based assignment was further validated by physiological and behavioral data.Consistent with previous studies (37,39), CHC-classified foragers were lighter than nurses (both for whole body and gut weight) and also exhibited a lower number of Actin gene copies in the gut as measured by qPCR on gut DNA extractions, suggesting differences in cell numbers between nurse and forager guts [Fig.4b, c, and d; linear mixed effects models fitted by REML with colony of origin and cage as nested random effects: n = 120, whole body weight, F (1,116) = 12.61, P = 0.0006, gut weight, F (1,118) = 15.68,P = 0.0001, log(Actin copies), F (1,110) = 13.60,P = 0.0004].Forager-like gnotobiotic bees also had more degenerated HP glands than nurses (Fig. 4e; Kruskal-Wallis test: χ² = 8.07, df = 1, P = 0.005).Finally, gnotobiotic bees which had a CHC profile typical of foragers at the end of the automated behavioral tracking experiment (10-day-old bees) had interacted significantly less frequently with nestmates, performed more foraging trips, initiated foraging trips earlier, and spent more time in the foraging arena over the week of tracking than bees that had a nurse-like CHC profile at the end of the experiment (Fig. 4f, g, h, and i; linear mixed effects models fitted by REML with experimental replicate and sub-colony as nested random effects: social interactions, n = 171, F (1,159) = 20.17,P < 0.0001; foraging trips, n = 171, F (1,166) = 9.18, P = 0.003; age at first foraging trip, n = 111, F (1,108) = 6.33,P = 0.013; proportion of time spent in the foraging arena, n = 171, F (1,169) = 53.44,P < 0.0001).
Nurses and foragers are also known to differ substantially in brain gene expression (35).Consistent with this, the comparison of the RNA-sequencing profiles of CHC-clas sified nurses and foragers revealed a differential expression of 894 genes (i.e., 7% of the transcriptome; Fig. 4j; Table S2).To assess whether the gut microbiota affects behavioral maturation-related gene expression, we compared the identity of these genes with those that were differentially expressed as a function of gut microbiota composition (91 genes [7]).The overlap (11 genes) between these gene lists was not greater than expected by chance (Fig. 4j; hypergeometric test: representation factor = 1.67,P = 0.06).Furthermore, differential gene expression by microbiota treatment was most pronounced in the antennal lobe and subesophageal ganglion region (as shown in reference 7 for the same experimental bees), while differential gene expression by behavioral maturation was most pronounced in the mushroom body and central complex region (Fig. 4k).Finally, in the gut, 15 genes were differentially expressed between CHC-classified nurses and foragers, of which only one featured among the 4,988 genes differentially expressed between the gut microbiota treatments (Fig. 4l).The overlap between these DEG lists was again not greater than expected by chance (representation factor = 0.16, P = 0.99).Together these results indicate that, across tissues, the transcriptomic effects of the gut microbiota are not directly related to behavioral maturation.

Co-housing homogenizes CHC profiles and produces skewed distributions of nurses and foragers
Experiments applying treatments to bees (e.g., microbiota, antibiotics, and pesticides) often involve housing bees in shared environments ("cages").Co-housing may influence the variables tested due to non-independence (e.g., social interactions) of bees sharing the same cage.We therefore assessed whether uncontrolled co-housing effects on behavioral maturation, which have been previously reported (44), could provide an explanation for inconsistencies between our study and those reporting an effect of the microbiota on host CHC profile and weight gain (32,33).
We first tested whether co-housing could drive convergence in CHC profiles.To do this, we analyzed the effect of caging on CHC profiles in the experiment where bees from a single hive were placed in 20 different cages (the experimental design involved 10 cages per treatment allowing us to assess the effects of caging and microbiota treatment simultaneously).Bees collected from the same cage (five to six bees per cage) had CHC profiles more similar than bees from different cages (PERMANOVA using Bray-Curtis dissimilarities, n = 118, F (1,117) = 2.67, R 2 = 0.31, P = 0.001).Additionally, we tested whether the proportions of CHC-classified nurses and foragers were more skewed across cages than expected by chance using the CHC data collected in the automated behavioral tracking experiment because we had CHC data from a minimum of 8 bees in each of the 18 cages.This analysis revealed a significant co-housing effect on the proportion of individuals that matured into foragers (range from 0 to 0.6; chi-square test: χ² = 30.78,df = 17, P = 0.02).
Because co-housing can lead to skewed proportions of nurses and foragers, individuals within a given cage should not be treated as independent values to study the role of the gut microbiota on behavioral maturation-related phenotypes.Given that previous studies did not control for such an effect, we tested whether the gut microbiota affected the distribution of nurses and foragers across our experiments.For both the 8-and 10-day-old bees in the RNA-sequencing and automated behavioral tracking experiments, there was no significant difference in the proportion of nurses and foragers (classified based on CHC profiles) between MD bees and bees of the different colonization treatments (Fig. 5a and b; RNA-sequencing experiment: generalized linear mixed model [GLMM] fitted by maximum likelihood using a binomial distribution with colony of origin and cage as nested random effects, n = 120, CL_Bifi, estimate = −0.69,se = 0.67, z = −1.02,P = 0.31; CL_13, estimate = −0.91,se = 0.67, z = −1.35,P = 0.18; CL, estimate = −0.57,se = 0.66, z = −0.86,P = 0.39; automated tracking experiment: GLMM fitted by maximum likelihood using a binomial distribution with experimental replicate and sub-colony as nested random effects, n = 177, estimate = 0.23, se = 0.53, z = 0.44, P = 0.66).This is consistent with the observation that there was no difference between the microbiota treatments in the time bees spent in the foraging arena and the total number of foraging trips performed per bee in the automated tracking experiment.Similarly, there was no significant effect of the gut microbiota on the proportion of foragers at both day 8 and day 15 in the experiment designed to assess the effect of co-housing on CHC profiles (Fig. 5c; GLMM fitted by maximum likelihood using a binomial distribution with cage as random effect: n = 118, time, estimate = −0.87,se = 0.58, z = −1.49,P = 0.14; treatment, estimate = 0.72, se = 0.72, z = 1, P = 0.32; time × treatment, estimate = −0.74,se = 0.84, z = −0.88,P = 0.38).There was also no significant effect of the gut microbiota on the proportion of foragers at the end (day 11) of either our weight gain experiment with a bee bread diet (Fig. 5d; GLMM fitted by maximum likelihood using a binomial distribution with colony of origin and cage as nested random effects: n = 102, estimate = −0.09,se = 1.26, z = −0.07,P = 0.94) or the weight gain experiment with a pollen diet (Fig. 5e; GLMM fitted by maximum likelihood using a binomial distribution with colony of origin as random effect: n = 80, estimate = 1.30, se = 1.62, z = 0.80, P = 0.42).Finally, the longitudinal CHC experiment, for which we had collected CHC data every 2 days from adult emergence until day 11 (Fig. 2c), allowed us to more precisely classify the bees that were transitioning between nurse and forager states, as we could identify intermediate groups in the clustering and ordination analyses (Fig. 2c; Fig. S1D and E).There was again no statistically significant difference in the proportion of foragers between CL and MD treatments (Fig. 5f; cumulative link mixed model with hive as random effect: n = 90, treatment, LR = 2.04, P = 0.15; time × treatment, LR = 0.94, P = 0.33).
To further confirm that there was no effect of the gut microbiota on the proportion of foragers, we performed a global analysis comparing the CL and MD treatments across all data sets (n = 602 individuals classified as either nurses or foragers across 94 cages, 35 hives, and 6 experiments).For this, we also assessed the effect of the number of co-housed bees at the time of sampling.There was a clear effect of time and group size but no effect of gut microbiota treatment on the proportion of foragers nor an interaction between time and treatment (GLMM fitted by maximum likelihood using a binomial distribution with experiment, colony of origin, and cage as nested random effects: n = 602, time, estimate = −0.17,se = 0.07, z = −2.50,P = 0.013; group size, estimate = 0.02, se = 0.01, z = 3.35, P < 0.001; treatment, estimate = 0.92, se = 0.87, z = 1.06,P = 0.29; time × treatment, estimate = −0.06,se = 0.09, z = −0.72,P = 0.47).These results suggest that the gut microbiota has no effect on the proportion of foragers of honeybee colonies.

DISCUSSION
The honeybee is a powerful model to advance evolutionary and mechanistic under standing of host-microbe interactions (26,45).Previous studies have identified several effects of gut microbes on honeybee phenotypes, including weight (32), CHC profile (33), learning and memory (29,30), and frequency and patterning of social interactions (7).All these phenotypes change during behavioral maturation (37,38,46), with, for example, foragers being lighter and having different CHC profiles than nurses.This raises the question of whether the reported effects may be indirect (i.e., a consequence of an effect of the microbiota on behavioral maturation).Our laboratory studies showed that the gut microbiota accelerates the time at which bees make their first trip to a foraging arena by a few hours.However, the total number of foraging trips was not affected, suggesting that the effect of the microbiota is mild and transient, and possibly overshadowed by other mechanisms that regulate behavioral maturation (47)(48)(49)(50).Division of labor among honeybee colony members is highly plastic and known to be regulated by feedback mechanisms dependent on diet and on the behavior of nestmates and their pheromonal secretions (47)(48)(49)(50).Hence, the microbiota seems to play only a minor role, if any, in determining the proportion of foragers in groups of co-housed bees.The role of other factors in our study is also highlighted by our finding that within cages the proportion of foragers is more skewed than expected by chance.Consistent with these behavioral analyses, our data also showed that the microbiota has no significant effect on the proportion of individuals that transition to a forager state and on various physiological hallmarks of behavioral maturation, such as CHC profile, gut or body weight, the expression of behavioral-maturation-related genes, or hypopharyngeal gland development.Therefore, except for an effect on the timing of first foraging event, our study revealed no influence of the microbiota on hallmarks of behavioral maturation.However, we cannot exclude that the microbiota may affect host traits that we did not quantify and that are also associated with behavioral maturation.The effect of the microbiota on the timing of first foraging event may alternatively be due to improved cognitive abilities (29,30) promoting boldness and a more general tendency for exploratory behavior.Before embarking on trips to collect food, honeybees typically perform orientation flights and become more explorative and phototactic (51,52).Our laboratory assay could not discern between these foraging-related behaviors.Whether the tendency of colonized honeybees to embark earlier on trips to the foraging arena in the laboratory indicates an effect of the microbiota on the onset of foraging in the field will therefore require further testing.A recent study monitoring the foraging behavior of gnotobiotic honeybees mono-inoculated with a subset of core gut microbes in the field (53) reported that Bifidobacterium asteroides increased foraging rates, while Bombilactobacillus mellis and Lactobacillus melliventris decreased foraging rates.None of these microbes affected the onset of foraging behavior.Consistent with our study, there was no effect on the proportion of bees performing foraging trips.However, this study is not directly comparable to ours due to the different inoculation treatments (the complete gut microbiota was inoculated in our study, as opposed to mono-inoculations in reference 53) and because our study was performed under controlled laboratory conditions, whereas the other study was conducted in the field.
Our results are in contrast to two previous studies which reported that the honeybee gut microbiota affects CHC profile (33) and promotes weight gain (32).We found that honeybees kept in the same laboratory cage can either take a nurse-like or a forager-like state with correlated changes in physiology and behavior, including CHC profile and body and gut weight.These two types of bees occur in skewed proportions across experimental cages (i.e., individuals within a cage are more similar to each other than individuals between cages).It is likely that previous studies have used only one or few cages per treatment, meaning that the reported associations could stem from co-hous ing effects.For example, Vernier et al. (33) concluded that the honeybee gut microbiota affects the CHC profile of bees.This study involved a series of experiments that identified gut microbiota-associated changes in CHC profile and acceptance behavior of bees.No information is provided on the replication of the cage setup of the laboratory experiments testing the effects of the gut microbiota on CHC profile.Consequently, it is impossible to determine whether the reported differences in gut microbiota composition and CHC profile in multivariate analyses were due to the experimental treatment or co-housing effects (e.g., social interactions among bees sharing a cage reducing variation in CHC profiles and skewed behavioral maturation producing spurious differences between treatments).However, a re-analysis of the CHC profiles from the key experiment comparing bees inoculated with either live or dead bacterial suspensions showed that, as in our experiments, bees had segregated into nurses and foragers (Fig. S2A) and that there were twice as many foragers in the live inoculum than in the heat-killed inoculum treatment (heat killed: 6 foragers and 10 nurses; live: 11 foragers and 5 nurses; Fig. S2A), driving most of the difference in CHC ordination space (Fig. S2B).Whether the increase in foragers in the live inoculum treatment is due to an effect of the microbiota cannot be determined without details of cage replication.In that respect, it should be noted that an effect of the microbiota is unlikely because 16S rRNA gene amplicon sequencing data from this experiment show that bees in both treatment groups had been colonized by core gut microbes and that the microbiota treatments determined a statistically significant difference in the relative proportion of only a few opportunistic bacteria (absolute bacterial loads were not assessed in this study; Fig. S2C; Table S3).This was unlikely to be due to the amplification of the 16S rRNA gene from dead bacteria in the heat-killed treatment because the same pattern was apparent in the experiment in which bees were mono-inoculated with the non-core microbes Acineto bacter and Lonsdalea quercina (Fig. 3b in reference 33): bees in both treatments had been colonized by all core microbiota members, and only Lonsdalea quercina showed statistically significant differences in relative abundance between treatments (Fig. S2D; Table S4).Additionally, there was no bacterium classified as Acinetobacter, suggesting that its inoculation was unsuccessful (Fig. S2D).
Similarly, our results contrast with those of Zheng et al. (32) who reported a higher weight gain in microbiota-colonized than microbiota-depleted bees.We could not find such effect across three independent experiments employing larger sample size and cage replication.If anything, in one of our longitudinal experiments, there was significant effect of the microbiota in the opposite direction, with CL bees exhibiting reduced weight compared to MD bees from day 7 onward.However, this effect, unique to one of our three experiments, may have been due to the fact that there were slightly more foragers in the cages assigned to the CL group in this experiment compared to MD cages (Fig. 5e; this difference was not statistically significant).According to Zheng et al. (32), bees originated from four hives were hosted in different cages.Unfortunately, we could not access the original data, precluding testing for co-housing effects.It is still possible that other factors play a role for the discrepancy of these results, such as host genotype, gut homogenate used, or seasonal differences between bees.
Our study reveals that bees within cages can be at different stages of their behavioral maturation-a fact that has been previously reported in "single-cohort" colonies (i.e., outdoor hives composed of a few hundred or thousand age-matched young bees [54,55]) and in groups of age-matched bees kept in the laboratory (44).This had been neglected in studies of the honeybee gut microbiota.Social effects on behavioral maturation can confound gnotobiotic bee experiments and need to be controlled for by randomly sampling individuals from separate cages and increasing the number of replicate cages beyond what has been used in many previous studies.In conclusion, our study indicates that the gut microbiota induces an earlier onset of foraging-related behavior but does not regulate the proportion of foragers in honeybee colonies.Additionally, our findings suggest that previous reports of an association between the gut microbiota and weight gain and CHC profiles are likely due to bees being more similar within than between cages because of the effect of social interactions or other cage-specific characteristics.

Experimental overview
In the present study, we have reanalyzed behavioral and RNA-sequencing data of the automated behavioral tracking and RNA-sequencing experiments published in reference 7 to investigate foraging-like behaviors (Fig. 1) and to compare brain and gut gene expression between CHC-classified nurses and foragers (Fig. 4j, k, and l).We also used the honeybee collections from these two experiments to obtain CHC profile and hypophar yngeal gland size data (Fig. 2a and b; Fig. 3c and d), as well as weight data (Fig. 3a and b).We then performed four new experiments: two experiments aimed at further investigat ing the effect of the gut microbiota on CHC profiles (a time-series CHC experiment, Fig. 2c, and an additional experiment to assess the effect of co-housing on CHC profiles, Fig. 2d) and two longitudinal experiments to assess the effect of the gut microbiota on bee weight gain (Fig. 3e and f).

Rearing of gnotobiotic bees
Across the experiments, bees were reared as previously described (7,23,56).Using sterile forceps, we extracted melanized dark-eyed pupae from capped brood cells and placed them in groups of 25-30 into sterilized plastic containers lined with moist cotton.We kept these pupae in an incubator at 70% relative humidity and 34.5°C in the dark for 3 days, then transferred newly emerged worker bees into corresponding cup-cages built using a sterile plastic cup placed on top of a 100-mm petri dish.To colonize bees, an aliquot of a gut homogenate was thawed and diluted 10× in 1× PBS and subsequently 1:1 in sugar water.Microbiota-depleted controls were provided only a 1:1 PBS:SW solution.To inoculate bees, three 100 µL droplets of treatment solution were added to the bottom of each cage.Bees were then kept in their cages in an incubator at 70% RH and 30°C in the dark (except for the bees in the automated behavioral tracking experiment, which were kept under the tracking systems in groups of ca. 100 bees to monitor their behavior; see below and reference 7 for additional details) and continuously fed by providing sterile SW and pollen (except for one of the longitudinal weight gain experiments where bee bread was used instead) ad libitum.

Preparation of gut homogenates to inoculate bees
For each experiment, we randomly collected five nurse bees from each of three hives.We anesthetized bees on ice, dissected their guts, and placed them individually in 1 mL 1× PBS containing 0.75-1 mm sterile glass beads.Guts were homogenized at 6 ms −1 for 45 s using a FastPrep-24 5G homogenizer (MP Biomedicals).The five gut homogenates were pooled by hive of origin, and serial dilutions of these pools from 10 −3 to 10 −12 were plated onto BHIA, CBA + blood, and MRSA + 0.1% L-cys + 2% fructose media using the drop method (10 µL droplets).These plates were then incubated under both anaerobic and microaerobic conditions to verify bacterial growth.Additionally, we prepared lysates of the homogenates by mixing 50 µL of each homogenate with 50 µL lysis buffer, 5 µL proteinase K (20 mg mL −1 ), and 5 µL lysozyme (20 mg mL −1 ) and incubating these mixtures for 10 min at 37°C, 20 min at 55°C, and 10 min at 95°C in a PCR machine.Lysates were centrifuged for 5 min at 2,000 g, and the supernatants used as templates for diagnostic PCR.We performed diagnostic PCRs using specific primers (as done in reference 7) to verify the absence of known honeybee pathogens (Nosema apis, Nosema ceranae, trypanosomatids, Serratia marcescens) and fungal growth in bee guts, as well as the presence of bifidobacteria as initial validation that the homogenates contained members of the core gut microbiota.Homogenates with the lowest amplification of pathogen DNA were selected, spiked with glycerol to a final concentration of 20%, and stored at −80°C.Prior to using a selected homogenate in an experiment, we thawed an aliquot and plated it on various media as described above to verify that the homoge nates were viable after storage at −80°C.For the time-series CHC experiment and the two weight gain experiments, we used the same homogenate that had been previously prepared for the RNA-sequencing experiment.The gut homogenates for the automated behavioral tracking experiment and the single-colony CHC experiment were prepared anew.

Measurement of fresh body and gut wet weight
At the end of the RNA-sequencing experiment (see reference 7 for additional details), we measured fresh body and gut wet weight of the 8-day-old bees across the 58 experimen tal cages (reared from10 different hives and randomly assigned to 4 gut microbiota treatment groups).To do this, we anesthetized bees on ice and weighed them using an electric balance sensitive to 0.0001 g.We then dissected their guts as described in reference 7, placed them in previously weighed 2 mL screw-cap tubes, and used the same electric balance to weigh them.The weight of the tube was then subtracted from the total measurement.
Next, we performed two longitudinal weight gain experiments.For each experiment, we reared gnotobiotic bees from three hives in six different cages (one per treatment per hive).Bees from each cage were paint marked with unique combinations of colors, and their body weight was measured every 2 days for 10 days (including the day of adult emergence and treatment inoculation).At each time point, the cages were placed on ice to anesthetize bees, and each bee was weighed using an electric balance sensitive to 0.0001 g.At the end of the experiment (day 10), bees were anesthetized on ice, snap frozen in liquid nitrogen, and stored at −80°C for subsequent CHC analyses.

Hypopharyngeal gland size
During brain dissection for RNA sequencing, we quantified the size of the hypopharyng eal glandular system of 28 bees using a semi-quantitative scale from 1 to 5 (from the most degenerated to the most developed), assigning the score blindly with respect to gut microbiota treatment or CHC group.

Chemical analysis of cuticular hydrocarbons by GC-MS
We collected cuticular hydrocarbon data from bees across multiple experiments.These included the bees at the end of the RNA-sequencing experiment (n = 120) and auto mated behavioral tracking experiment (n = 177), when bees were 8 and 10 days old, respectively.We also collected CHC data across the two longitudinal weight gain experiments described above (11-day-old bees; pollen experiment, n = 80; bee bread experiment, n = 102).We then designed a longitudinal experiment to follow the development of the CHC profile of gnotobiotic bees produced from 9 different hives and kept in 18 different cages (one cage per treatment per hive).We collected one bee per cage every 2 days for CHC analyses, starting from the day of treatment inoculation until bees were 11 days of age (n = 108).Finally, we performed an additional experiment rearing gnotobiotic bees from a single hive in 20 distinct cages and collecting 3 bees per cage after 8 and 15 days (n = 118).All bees were stored at −80°C until CHC analyses were performed.
The thorax and abdomen after gut extraction, or only the abdomen for samples of the automated behavioral tracking experiment (thoraxes had been previously used for hemolymph extraction), were submerged in pure hexane for 10 min.These extracts were evaporated to a residue of approximately 100 µL.The hexane extracts were run with a DB-5 capillary column (0.25 mm × 30 m × 0.25 μm; JW Scientific) on an Agilent 6890-5975 GC-MS at the University of Würzburg (RNA-sequencing experiment) or with an HP-5MS column (0.25 mm × 30 m × 0.25 μm; Agilent) on an Agilent 8890-5977B GC-MS at the University of Lausanne (all other experiments).A temperature program from 60°C to 300°C with 5°C/min and finally 10 min at 300°C was used for the RNAsequencing experiment data, with data collection starting 4 min after injection.The mass spectra were recorded in the electron ionization mode, with an ionization voltage of 70 eV and a source temperature of 230°C.The chromatography protocol at the Univer sity of Lausanne was shortened by ramping the oven from 65°C to 215°C at 25°C/min and then to 300°C at 8°C /min.Data were acquired and processed with the ChemSta tion software v.F.01.03.2357 (Agilent Technologies).Identification of the compounds was accomplished by comparison of library data (NIST 20) with mass spectral data of commercially purchased standards for n-alkanes, diagnostic ions, and retention indices.

CHC data analyses
To calculate the relative abundance of CHC compounds, the area under each compound peak on the GC was quantified through integration using the ChemStation software and divided by the total area under all CHC peaks.The raw data were aligned using the R package GCalignR v.1.0.5 and afterward analyzed using the packages vegan v.2.6-4 and dendextend v.1.17.1.Polar compounds and contaminations were identified using the mass spectral data (all non-hydrocarbons) and removed from the data set.Afterward, we removed compounds that were not present in at least half the samples of one treatment or that were only present in trace amounts (<0.1%) in all samples.Lastly, samples which had a too low concentration of CHC compounds (due to failed extractions) were excluded from the analysis.All analyses were done using RStudio v.1.4.1717 and R v.4.1.0and the package ggplot2 v.3.4.2 for visualization.Area under the peak values was converted to relative proportions, after which we calculated Bray-Curtis dissimilarities between samples and performed non-metric multidimensional scaling (NMDS) ordination analyses and PERMANOVA with 999 permutations to assess differences between experimental groups.To account for sampling multiple individuals from the same cages, we calculated the multivariate centroids from each cage using the Betadisper function (package vegan) and tested for the main treatment effect using the resulting matrices, while within-subject effects were tested separately using the original data sets.For PERMANOVA analyses of the CHC profiles in the time-series experiment, we removed data from the day of treatment inoculation (day 1) as we did not expect the treatment to have produced immediate effects on CHC profiles (repeating the analysis including day 1 did not change the statistical results qualitatively).Because we analyzed the CHC profiles from the RNA-seq experiment in two separate GC-MS runs, we used the removeBatchEffect function in edgeR v. 3.34.1 (57) to remove the batch effect prior to plotting the NMDS ordination.We used the hclust function of the base R package "stats" to perform hierarchical cluster analyses of Euclidean distances between CHC profiles using the Ward's criterion prior to plotting heatmaps.Nurses and foragers were then identified based on the resulting clusters.We used generalized linear mixed models fitted by maximum likelihood using a binomial distribution to assess the effect of gut microbiota treatment on the proportion of these CHC-classified nurses and foragers.We always accounted for sampling multiple individuals from the same cages by adding cage as random effect to the models.Based on the hierarchical clustering and ordination analyses of the CHC data collected in the time-series experiment, we were able to identify intermediate clusters (newly emerged bees, bees transitioning to the nurse cluster, nurses, and bees transitioning to the forager cluster and foragers).To test the effect of gut microbiota treatment on the proportion of these CHC clusters, we used a cumulative link mixed model with treatment, time, time × treatment, and hive as fixed effects and cage as random effect using the clmm function in the package "ordinal" v.2022.11-16.To do this, we again excluded data from day 1 (this did not change the statistical results qualitatively).

Quantification of foraging tendency under the automated behavioral tracking systems
In the automated behavioral tracking experiment (see reference 7 for additional details on experimental procedures and data post-processing), bees were housed in a doublebox setup, meaning that they had access to a nest box (kept in constant darkness) connected via a tube to a foraging box (subject to day-night condition cycles; Movie S1).The nest was kept at 30°C and 70% RH, while the foraging arena had cyclic conditions with light and temperature of 25°C during the day, and dark condition and 18°C during the night.Transitions in the foraging arena initiated at 4:00 and 16:00 and lasted for 4 h, during which the climate system performed a linear interpolation between the two states.Both boxes contained three 2 mL vials filled with SW, which were continuously replaced, and a trough filled with 1 g of pollen.Bees were placed into the nest box at the start of the experiment, allowing us to quantify three metrics for each individual: (i) the time (and hence the age) at which the individual first ventured into the foraging box, (ii) the total proportion of time spent in the foraging arena (i.e., total frames in which an individual was detected in the foraging arena/total number of frames in which the individual was detected in either box), and (iii) the number of box switches (i.e., each time the individual moved from the nest box to the foraging box and vice versa).We performed all statistical analyses in R v.4.1.0.To assess the effect of the gut microbiota on behavioral variables (average values for each sub-colony), we first checked whether the differences between paired values were normally distributed using the Shapiro-Wilk normality test and then ran either paired t-tests or Wilcoxon matched-pairs signed-rank tests.

RNA-sequencing data analyses
We reanalyzed our previously published RNA-sequencing data (7) to identify differentially expressed genes between CHC-classified nurses and foragers and assess the overlap between these DEGs and those that we had previously identified in gut microbiota treatment comparisons from the same bees (gut, n = 38; antennal lobes and subesophageal ganglion, AL, n = 39; mushroom bodies and central complex, MB, n = 39; optic lobes, OL, n = 38).See reference 7 for details on data processing to obtain the raw read counts, which we reanalyzed in the present study, and for the differential expression analyses of gut and brain between gut microbiota treatment groups.For comparisons of gene expression between CHC-classified nurses and foragers, we used the same parameters as done previously for the gut microbiota comparisons in reference 7. Briefly, we filtered out genes not represented by at least 20 reads in a single sample using the filterByExpr function in edgeR (57).Next, we used the Limma Bioconductor package v. 3.48.3 (58) for differential expression analyses.For the gut, we used the formula 0 + CHC classification + batch, whereas, for the brain, we used the formula 0 + group + batch, where "group" represented every possible combination of brain region and nurse or forager group, and "batch" represented the different experimental and RNA-seq library preparation batches.As we had sampled multiple brain regions from the same individuals, we accounted for the individual random effect using the duplicateCorrelation function.For the brain, the contrasts between CHC-classified nurses and foragers were performed overall and within each brain region separately.P values were adjusted for multiple testing with an FDR of 5%.Hypergeometric tests were used to compare the overlap in DEGs by gut microbiota treatment and by CHC classification of nurses and foragers in both the gut and the brain.

FIG 1
FIG 1 The gut microbiota accelerates the onset of foraging-like behavior under an automated behavioral tracking system.(a) Timeline and experimental setup for the automated behavioral tracking experiment.In each experimental replicate, two groups of 100 gnotobiotic bees (either microbiota depleted or colonized) could move freely between two Plexiglass boxes connected by a plastic tube and hosted in separate climate-controlled chambers.(b) Average number of trips between the nest and the foraging arena per bee for each sub-colony in the automated behavioral tracking experiment.(c) Average proportion of time spent in the foraging arena per bee for each sub-colony.(d) Average age at which bees made their first trip to the foraging arena for each sub-colony.Lines connect paired sub-colonies and are colored by experimental replicate.

FIG 2
FIG 2 The gut microbiota does not affect CHC profile.(a) Non-metric multidimensional scaling (NMDS) of Bray-Curtis dissimilarities between CHC profiles in the automated behavioral tracking experiment (MD, n = 88; CL, n = 89).(b) NMDS of Euclidean distances between CHC profiles in the RNA-sequencing experiment, after removal of batch effects from two separate GC-MS runs (MD, n = 31; CL_Bifi, n = 29; CL_13, n = 30; CL, n = 30).(c) NMDS of Bray-Curtis dissimilarities between CHC profiles in the longitudinal experiment (MD, n = 54; CL, n = 54).(d) NMDS of Bray-Curtis dissimilarities between CHC profiles in the single-colony experiment (MD, n = 59; CL, n = 59).Samples are colored by gut microbiota treatment, and shapes indicate nurses and foragers in panels a, b, and d.Samples in panel c are colored by age of the sampled bees, and shapes indicate treatment group.

FIG 3
FIG 3 The gut microbiota does not affect weight gain and hypopharyngeal gland size.Boxplots reporting fresh body weight (a) and gut wet weight (b) by gut microbiota treatment group in the RNA-sequencing experiment (8-day-old bees).Boxplots show the median and first and third quartiles.Whiskers show the extremal values within 1.5×, the interquartile ranges above the 75th and below the 25th percentile.(c) Photographs showing examples of maximally developed and degenerated hypopharyngeal glands and (d) proportion of hypopharyngeal gland sizes across gut microbiota treatment groups in the RNA-sequencing experiment.NS, not significant.(e, f ) Fresh body weight growth curves of individual bees colored by gut microbiota treatment group, shown separately for the bee bread (e) and pollen (f ) experiments.Values are proportions of initial body weight at the time of gut microbiota colonization (1-day-old bees).Thicker lines represent mean values, and bars indicate SD.

FIG 4
FIG 4 Gnotobiotic bees reared in cages diverge into nurses and foragers showing differences in physiology and behavior.(a) Heatmap of relative abundance of detected CHCs on the cuticle of gnotobiotic bees in the RNA-sequencing experiment (8-day-old bees; n = 120; shown in gray in the annotation column toward the left) and conventional nurses (n = 51) and foragers (n = 9) randomly collected from the same hives in blue and pink, respectively.The dendrogram toward the left shows clustering of CHC profiles based on Euclidean distances using Ward's criterion.(b) Fresh body weight, (c) gut wet weight, (d) number of Actin copies in the gut, and (e) hypopharyngeal gland size of CHC-classified nurse and forager gnotobiotic bees in the RNA-sequencing experiment (8-day-old bees).(f-i) Boxplots reporting the number of head-to-head interactions (normalized by group size) (f ), trips to the foraging arena (g), the age at first foraging trip (h), and the proportion of time spent in the foraging arena (i) over the week of tracking for the gnotobiotic nurses and foragers classified based on their CHC profiles at the end of the automated behavioral tracking experiment (10-day-old bees).***P < 0.001; **P < 0.01; *P < 0.05; NS, not significant.Numbers at the bottom of boxplots and stacked bars in panels b to i indicate sample sizes.(j) Venn diagram reporting overlap in the brain between the differentially expressed genes (DEGs) associated with the gut microbiota (as identified in reference 7) and those associated with behavioral maturation (CHC-classified gnotobiotic nurses vs foragers).(k) Venn diagram reporting overlap in DEGs in brain region-specific comparisons of CHC-classified gnotobiotic nurses vs foragers.(l) Venn diagram reporting overlap in the gut between the DEGs associated with the gut microbiota (as identified in reference 7) and those associated with behavioral maturation.The brain icons were created with BioRender.com.

FIG 5
FIG 5 Proportions of nurses and foragers across the experiments.Stacked bars report the percentage of CHC-classified gnotobiotic nurses and foragers (based on clustering of Euclidean distances in CHC profiles using the Ward's criterion) in the RNA sequencing (a), automated behavioral tracking (b), single-colony CHC (c), longitudinal weight gain with either bee bread (d) or pollen diet (e), and time-series CHC (f ) experiments.Numbers at the bottom of stacked bars indicate sample sizes.NS, not significant.