BMC Genomics RESEARCH ARTICLE Open Access Honey bee (Apis mellifera)larval pheromones may regulate gene expression related to foraging task specialization Rong Ma Juliana Rangel and Christina M.Grozinger Abstract Background:Foraging behavior in honey bees (Apis mellifera)is a complex phenotype that is regulated by signals.How toewngeratleithemoescgarleeliomnodulbtreforagng nbrain the queen and larvae.Larval pheromones can also stimulate foragers to leave the colony to collect pollen.However. the mechanisms und nng this rapid behaviorl plasticityin foragers that specialize pollen over ectar,and erent ,re Here edooaalpheomonsboodheooneB时hdeoakee0menghone hypothesized that both pheromones would alter expression of genes in the brain related to foraging and would differentially impact brain gene expression depending on foraging specialization. Results:Combining data reduction,clustering,and network analysis methods,we found that foraging preference ectar vs.pollen)and pheromone exposure are each associated with specific brain gene expression profiles ore,pneromone e e nas a strong enect on genes that are tudies revealed significant overlaps for both pheromone communicationand foring task of foraging-relatec nsights into how social signals and task specialization are potentially integrated at the molecular level,and s the bothat brnmy playnoybee behavrme Keywords:Animal behavior,Behavioral plasticity,Communication,Differential gene expression,Gene networks, Larval pheromone signals,Task specialization One of the hallmarks of insect sociality is division of ral studies have deary labor,whereby group members specialize on different demonstrated that complex animal behaviors,induding so tasks that tial to group survival and cial interactions,are regulated by transcriptional,ne eural,and reprod on emn vever,the n gy.Hu mediating more rapid NBMC data made
R E S EAR CH A R TIC L E Open Access Honey bee (Apis mellifera) larval pheromones may regulate gene expression related to foraging task specialization Rong Ma1* , Juliana Rangel2 and Christina M. Grozinger1 Abstract Background: Foraging behavior in honey bees (Apis mellifera) is a complex phenotype that is regulated by physiological state and social signals. How these factors are integrated at the molecular level to modulate foraging behavior has not been well characterized. The transition of worker bees from nursing to foraging behaviors is mediated by large-scale changes in brain gene expression, which are influenced by pheromones produced by the queen and larvae. Larval pheromones can also stimulate foragers to leave the colony to collect pollen. However, the mechanisms underpinning this rapid behavioral plasticity in foragers that specialize in collecting pollen over nectar, and how larval pheromones impact these different behavioral states, remains to be determined. Here, we investigated the patterns of gene expression related to rapid behavioral plasticity and task allocation among honey bee foragers exposed to two larval pheromones, brood pheromone (BP) and (E)-beta-ocimene (EBO). We hypothesized that both pheromones would alter expression of genes in the brain related to foraging and would differentially impact brain gene expression depending on foraging specialization. Results: Combining data reduction, clustering, and network analysis methods, we found that foraging preference (nectar vs. pollen) and pheromone exposure are each associated with specific brain gene expression profiles. Furthermore, pheromone exposure has a strong transcriptional effect on genes that are preferentially expressed in nectar foragers. Representation factor analysis between our study and previous landmark honey bee transcriptome studies revealed significant overlaps for both pheromone communication and foraging task specialization. Conclusions: Our results suggest that, as social signals, pheromones alter expression patterns of foraging-related genes in the bee’s brain to increase pollen foraging at both long and short time scales. These results provide new insights into how social signals and task specialization are potentially integrated at the molecular level, and highlights the possible role that brain gene expression may play in honey bee behavioral plasticity across time scales. Keywords: Animal behavior, Behavioral plasticity, Communication, Differential gene expression, Gene networks, Larval pheromone signals, Task specialization Background One of the hallmarks of insect sociality is division of labor, whereby group members specialize on different tasks that are essential to group survival and reproduction [1, 2]. Understanding the proximate and ultimate mechanisms mediating social behavior, division of labor, and task specialization is a major focus of behavioral sociobiology [3–8]. Several studies have clearly demonstrated that complex animal behaviors, including social interactions, are regulated by transcriptional, neural, and physiological networks [9–12]. Moreover, several studies have suggested that behavioral ontogeny is mediated by differential regulation of core, well-conserved transcriptional or physiological “toolkits” that regulate behavioral modules [4, 13–19]. However, the mechanisms mediating more rapid shifts in behavior and task specialization have not been examined as thoroughly [20–22]. © The Author(s). 2019 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. * Correspondence: rongma.10@gmail.com 1 Department of Entomology, Center for Pollinator Research, Center for Chemical Ecology, Huck Institutes of the Life Sciences, Pennsylvania State University, University Park, PA, USA Full list of author information is available at the end of the article Ma et al. BMC Genomics (2019) 20:592 https://doi.org/10.1186/s12864-019-5923-7
(201920592 Page 2 of 15 As in many social insects,honey bee (Apis mellifera) pollen forag ging within an hour of exposure and lasting for 3 a form of age- based ta prod oping shift over the coure fetime 231.This hat FRo is duced early in laryal。 while br called age-based polyethism- -is regulated iust before pupation Both larv enetically and vironmer ally,and stem rs.In fact,bro the first weeks of their lives performing tasks within the including modulation of sucros relative safety of the hive,including tending to the needs ovary dev opment,foraging ontoe developing larvae (i.e oning and hypopharyngea forage,workers may further specialize by collectin ically.brood pheromones cause an increase in the number of ominantly one flora resource type (either pollen aging trip and size of pollen loads【24 and pr ty fo e xhibit dis behavioral,physiological,and trans criptional traits.For the transition of bees from performing within-hive roles t example, he colony,nectar foragers 6. some componen EBO is d b nto hon comb themselves 28 291 Nectar and pollen gers produce ethyl oleate.a co of BP both oragers lso differ ior [48,49.Queen 30]an 32 an comp od i 9.33361.ph are ime scal which they induc beh avioral changes in re arvae, and nent pheromone blends often have ca pheromon hav eral Prime rate long (ie hours)and the are also in d in regulating the size of the foraging labor force in the long especially in the brai 361 er m (i.e weeks】 queen hat。oni of h bee to phe ne ession of large numbers of genes in worker brains 33. 6.In contrast,releaser pheromones elicit rapid behavioral and chromatin remodeling 50).However,it is unclear if nges eith by activa n pheromone ng pa the ala m be (BD honey bees elicits aggressive behaviors against intruders late fora y activating the expression of imm diate early gen nes in the How the behavioral tran sitions acros s different temporal 13 ponent of que or heir unde ing genetic,epigen orkers by binding to an olfact or in the anten tobe dete ined nae,activating dopamine receptors in the brain,and regu In previous studies.the effects of bP on gene expression lating brain gene expr ession[33,40,41. brain expre arv r an 360 age, which provides a fascinating rtunity to unders regulation of behavior acros time scales Two larvae ing their foraging prefer ence.Consequently,we seek to more produced pheromon brood pheron e (BP nd (E)-beta precisely ch e the tr ces assoc cimene (EBO),have been sho toelicit rapi increases in with rapid changes in honey
As in many social insects, honey bee (Apis mellifera) workers exhibit a form of age-based task allocation in which behavioral repertoires incrementally expand or shift over the course of an individual’s lifetime [23]. This phenomenon—called age-based polyethism—is regulated both genetically and environmentally, and provides a tractable system in which to investigate temporal dimensions of behavioral plasticity [24, 25]. Honey bees spend the first weeks of their lives performing tasks within the relative safety of the hive, including tending to the needs of developing larvae (i.e., nursing), before transitioning to increasingly dangerous tasks near the nest entrance and beyond, including foraging [26]. Once they begin to forage, workers may further specialize by collecting predominantly one floral resource type (either pollen or nectar [27]), and their proclivity for pollen vs. nectar foraging can persist throughout their lives. Bees that specialize on nectar vs. pollen foraging exhibit distinct behavioral, physiological, and transcriptional traits. For example, upon returning to the colony, nectar foragers regurgitate collected nectar to nestmates waiting to process it, while pollen foragers pack their pollen loads into honeycomb themselves [28, 29]. Nectar and pollen foragers also differ in their neural and sensory responses to sugar [30] and pheromones [31, 32]. Pheromone communication in honey bees plays a key role in mediating behavioral transitions across time scales [9, 33–36]. Pheromones are typically categorized by the time scale at which they induce behavioral changes in receivers: primer pheromones cause slow, enduring changes in physiology, while releaser pheromones cause rapid, ephemeral responses. Primer pheromones generate longterm changes in behavior and physiology by altering patterns in gene expression, especially in the brain [9, 33–36]. For example, brood and queen pheromones delay the behavioral transition from nurses to foragers by altering the expression of large numbers of genes in worker brains [33, 36]. In contrast, releaser pheromones elicit rapid behavioral changes either by activating or modulating neural circuits, triggering molecular signaling pathways, or regulating gene expression [34, 37–39]. For example, the alarm pheromone in honey bees elicits aggressive behaviors against intruders by activating the expression of immediate early genes in the brain [34], while one component of queen pheromone, homovanillyl alcohol, elicits grooming behavior from workers by binding to an olfactory receptor in the antennae, activating dopamine receptors in the brain, and regulating brain gene expression [33, 40, 41]. Honey bee larval pheromones cause primer and releaser effects that blur the distinction between these categories, which provides a fascinating opportunity to understand regulation of behavior across time scales. Two larvaeproduced pheromones, brood pheromone (BP) and (E)-betaocimene (EBO), have been shown to elicit rapid increases in pollen foraging within an hour of exposure and lasting for 3 hours [42]. Both pheromones are produced by developing larvae but differ in the timing of their peak production, such that EBO is produced early in larval development while BP is produced later on, just before pupation [42]. Both larval pheromones cause additional behavioral and physiological effects in honey bee workers. In fact, brood pheromone induces the greatest number of known primer responses in honey bees, including modulation of sucrose response thresholds, ovary development, foraging ontogeny, foraging choice behavior, and hypopharyngeal gland development [43]. The effect of brood pheromones on forager behavior seems to be driven by an increase in pollen foraging. Specifically, brood pheromones cause an increase in the number of foraging trips and the size of pollen loads [42, 44], and this effect is not driven by task-switching from nectar to pollen foraging [42]. Both pheromones also increase the size of the foraging force of the colony in the long term, accelerating the transition of bees from performing within-hive roles to foraging [44–46]. Interestingly, some components of EBO and BP are also produced by honey bee adults as well. For example, EBO is also produced by mated queens [47], and foragers produce ethyl oleate, a component of BP [48]; both impact the ontogeny of foraging behavior [48, 49]. Queens and larvae both produce another BP component, ethyl palmitate, which inhibits ovarian development [37]. Although BP components are also produced in adults, the full blend of BP and EBO has only been described in honey bee larvae, and multi-component pheromone blends often have synergistic effects [37]. Overall, larval pheromones have a strong effect on pollen foraging but not nectar foraging in the short term (i.e., hours), and they are also involved in regulating the size of the foraging labor force in the long term (i.e., weeks). Chronic exposure of honey bee adults to pheromones that cause primer effects, including BP, have been shown to affect the expression of genes involved in methylation and chromatin remodeling [50]. However, it is unclear if similar epigenetic effects are observed when pheromones act at the short-term, releaser time scale. This is a fascinating system because both pheromones (BP and EBO) regulate foraging behavior, but at different temporal scales. How these behavioral transitions across different temporal scales are related, or how their underlying genetic, epigenetic, and physiological mechanisms interact to regulate foraging behavior, remains to be determined. In previous studies, the effects of BP on gene expression were evaluated on whole brain expression patterns from bees collected at five and fifteen days of age, after life-long exposure to brood pheromone [36]. However, in that study, the bees were collected without regard to their behavior, including their foraging preference. Consequently, we seek to more precisely characterize the transcriptional differences associated with rapid, pheromonally-regulated changes in honey Ma et al. BMC Genomics (2019) 20:592 Page 2 of 15
Ma et aL BMC Genomics 201920:592 Page 3 of 15 per sample gration centers of the brain igned to generate transcrin abundance for each anno (ie,mushroom bodies).Given that foragers have similar be avioral responses to BP an EBO 51],we hypot 1 tha annotation mel HAv3. ional file e S1) senting%of th EBO have more ounced effects on pollen foraging than 12 332 annotated honey bee genes pare nalysis was performed to te of nl type,and the interaction between pheromone and forage predictions:1)foragers specializing on pollen vs. There were 533 diffe expressed gene nBP erns or gene FDR n con nd EBO files n the brair of for ger h in the related to fo (Table me behavior at different time scales (ie transition to Additional file2:Tabl S2).Additionally,there were 13 aging)util ize similar mo DEGs that showed a s me pollen for. dad in rsthan nectar foragers Of the 269 DEGs related to pheromone treatment Combining differential gene expression, omone-related DEGs),there ere 58 DEGs between netwo analyse al l tween EBO an that arval herop ulate e of ion was almos enes involved in foraging tasks specifically.nectar and ,there were 14 genes that showed difference between BP and EBO samples.Because there were many xpres genes th wer ly expr regulated by similar sets of more shared DEGs than those The results of the study did not support the hypothesis from random expectation among pheromone treatments that larval pheromones affect gene expression more and between pheromone treat ents and fo ager type ongly in pol age cta re wer EBO rcgulate expression of a common subset of genes or sion profiles that significantly overlapped v with those of netic pathways(Table 2). ctar fo gers not pollen foragers.Our study DEGs were then mpared to n in bor ecular highli ted DEGs).While nd roe that brain gene expression plays in behavioral plasti overlaps betweer foraging-related and pheromone city across time scale It also probes the interface be related DEGs(Table 3),it is important to note that ne ween ephe more 60 tar vs pollen I agingwas a binary trait, so genes th bavioral and complexity across time mlated in the site fo ontext For example,ge nes that w upregulated in pollen foragers Results were downregulated in nectar fo ragers, and vic To ed in this stud小y ed DEG from mushroom bodies of pollen and in pollen foragers (and thus downregulated in nectar for posed to one of three phen none treatments paraffin oi agers)and those that were upregulated in nectar foragers ntrol,br ood pheron (BP),or E-beta-ocimene (EBO) do nregulate d in pollen forag DEGs fom ea nt caps
bee foraging, and to juxtapose these rapid changes with more stable differences in gene expression associated with task specialization, specifically in integration centers of the brain (i.e., mushroom bodies). Given that foragers have similar behavioral responses to BP and EBO [51], we hypothesized that these two pheromones regulate a common set of foraging genes in the brain (i.e., a foraging “toolkit”). Because BP and EBO have more pronounced effects on pollen foraging than nectar foraging [42, 45], we further hypothesized that larval pheromones affect foragers differentially depending on foraging task specialization. We thus compared the effects of EBO and BP exposure on foragers previously found to specialize on nectar or pollen to test the following four predictions: 1) foragers specializing on pollen vs. nectar foraging exhibit distinct patterns of gene expression in the brain, 2) BP and EBO stimulate the same transcriptional profiles in the brains of forager bees, 3) changes in the same behavior at different time scales (i.e., transition to and/or stimulation of pollen foraging) utilize similar molecular mechanisms, and 4) both larval pheromones have more pronounced effects on gene expression in pollen foragers than nectar foragers. Combining differential gene expression, clustering, and network analyses, our study presents several lines of evidence that support the predictions of the hypothesis that larval pheromones regulate a common suite of genes involved in foraging tasks. Specifically, nectar and pollen foragers showed distinct patterns of brain gene expression, BP and EBO do regulate a common set of genes, and changes in short-term and long-term shifts in foraging behavior are regulated by similar sets of genes. The results of the study did not support the hypothesis that larval pheromones affect gene expression more strongly in pollen foragers than nectar foragers, however. Contrary to our prediction, the data showed that exposure to larval pheromones produced gene expression profiles that significantly overlapped with those of nectar foragers but not pollen foragers. Our study provides insights into the molecular mechanisms underlying task allocation in honey bees, and highlights the possible role that brain gene expression plays in behavioral plasticity across time scales. It also probes the interface between ephemeral and more consistent changes in behavior to gain insight into mechanisms that permit behavioral plasticity and complexity across time. Results Transcript quantification The RNA samples collected in this study were extracted from mushroom bodies of pollen and nectar foragers exposed to one of three pheromone treatments: paraffin oil control, brood pheromone (BP), or E-beta-ocimene (EBO) (Fig. 1). The number of RNA-seq reads per sample ranged from 41 to 94 million, with an average of 65 million reads per sample. After quality filtering and adapter trimming, an average of 69% of the reads per sample were pseudoaligned to generate transcript abundance for each annotated transcript in the recently updated honey bee genome annotation (Amel_HAv3.1; Additional file 1: Table S1). Overall, 9179 genes were detected in all samples and were included in subsequent analyses, representing 74% of the 12,332 annotated honey bee genes. Differential gene expression Differential gene expression analysis was performed to characterize the effects of pheromone treatment, foragertype, and the interaction between pheromone and forager type. There were 533 differentially expressed genes (DEGs) whose expression varied in at least one contrast (FDR < 0.05), including 269 DEGs related to pheromone treatment and 326 DEGs related to forager type (Table 1; Additional file 2: Table S2). Additionally, there were 131 DEGs that showed a statistically significant interaction between forager type and pheromone treatment. The lists of all DEGs are provided in Additional file 2: Table S2. Of the 269 DEGs related to pheromone treatment (pheromone-related DEGs), there were 58 DEGs between BP and control samples, and 152 DEGs between EBO and control samples, indicating that EBO’s effect on gene expression was almost three times greater than that of BP. In addition, there were 148 genes that showed differences between BP and EBO samples. Because there were many genes that were differentially expressed in more than one contrast, we performed hypergeometric tests to further determine if there were more shared DEGs than those from random expectation among pheromone treatments, and between pheromone treatments and forager type. There were significant overlaps between all pairwise comparisons of pheromone treatment, indicating that BP and EBO regulate expression of a common subset of genes or genetic pathways (Table 2). Pheromone-related DEGs were then compared to DEGs that differed between nectar and pollen foragers (foraging-related DEGs). While we found significant overlaps between foraging-related and pheromonerelated DEGs (Table 3), it is important to note that nectar vs. pollen foraging was a binary trait, so genes that were upregulated in one foraging context were necessarily downregulated in the opposite foraging context. For example, genes that were upregulated in pollen foragers were also downregulated in nectar foragers, and viceversa. To further explore these results, we split the foraging-related DEGs into those that were upregulated in pollen foragers (and thus downregulated in nectar foragers) and those that were upregulated in nectar foragers (and thus downregulated in pollen foragers), and again looked for overlaps with DEGs from each pheromone treatment. Interestingly, there were significant overlaps Ma et al. BMC Genomics (2019) 20:592 Page 3 of 15
Ma et al BMC Genomics (201920592 Page 4of 15 Colonies exposed Pollen and nectar RNAseq libraries to pheromone roragers collectec generate Control ·Pollen师师 898738 X 。Nectar师辰辰 (E)-beta-ocimene ·Pollen师乘辰 (Oci】 ·Nectar辰辰 ·Pollen师师元 (BP) ·Nectar师乘元 Fig.1 Overview of e ng numbers of reads per sample and between pheromone-related DEGs and DEGs upregu lated in but not ol phe om thth ding that the result and ERO alata DEC with RD had fo effect was driven primarily by genes upregulated in nec- ipid biosynthesis and integral com onents of the membrane tar foragers relative to pollen foragers. (FDR<0 DEGs associated to EBO exposure were 10 und r integral treatment.we performed pentose phosphate pathway.There was a significant verp analysis for DEGs associated with pheromone treatment,for- of 39 genes between BP and EBO exposed foragers ager type,and their intera and these DEGs were type en tike y en DEGs related top ne treat t were enriched for in tegral components of membrane,fatty acid metabolism,and Hierarchical clustering and principal components analysis ipid biosynthe (FDR<0.05).Finally, EGs related to the PCA) chical clust ring analysis and PCA nd h The dEgs associated with either ebo or bp were also an n all variance-stabilized gene expression values of alvzed separately.Because there were few upregulated genes DEGs,hierarchical clustering grouped samples with Table 1 Numbers of DEG in all pairwise comparisons Upregulated Downregulated Pheromone Main Effect BP ys Control 12 46 EBO ys Contro 14 138 Food Main Effect Pollen ys Nectar 79 24 BP y Control and Food 30 EBOvCont and Food Up-and down
between pheromone-related DEGs and DEGs upregulated in nectar foragers (Table 4; hypergeometric tests, p < 0.01), but not between pheromone-related DEGs and DEGs upregulated in pollen foragers. In summary, BP and EBO both regulated foraging-related genes, and this effect was driven primarily by genes upregulated in nectar foragers relative to pollen foragers. To better understand the function of differentially expressed genes associated with forager type and pheromone treatment, we performed gene ontology (GO) enrichment analysis for DEGs associated with pheromone treatment, forager type, and their interaction. DEGs associated with forager type were significantly enriched for GO terms related to lipid metabolism and trypsin-like serine proteases (FDR < 0.05). DEGs related to pheromone treatment were enriched for integral components of membrane, fatty acid metabolism, and lipid biosynthesis (FDR < 0.05). Finally, DEGs related to the interaction of pheromone treatment and forager type were enriched for lipid biosynthesis and metabolism (FDR < 0.05). The DEGs associated with either EBO or BP were also analyzed separately. Because there were few upregulated genes associated with either pheromone, up- and down-regulated genes for each pheromone were pooled during pathway enrichment analysis, with the understanding that the results for pheromone could potentially be driven by down-regulated genes. DEGs associated with BP exposure were enriched for lipid biosynthesis and integral components of the membrane (FDR < 0.05). DEGs associated to EBO exposure were enriched for integral components of membrane, fatty acid biosynthetic processes, fatty acid metabolism, and the pentose phosphate pathway. There was a significant overlap of 39 genes between BP and EBO exposed foragers compared to controls (P < 0.05), and these DEGs were significantly enriched for metabolic pathways and fatty acid metabolism (FDR < 0.05). Hierarchical clustering and principal components analysis (PCA) Hierarchical clustering analysis and PCA were used to better understand broad patterns across all DEGs. Based on all variance-stabilized gene expression values of DEGs, hierarchical clustering grouped samples with Fig. 1 Overview of experimental design and sequencing. RNA-seq libraries were generated from nectar and pollen foragers exposed to three pheromone treatments. Three pooled pollen forager samples and three pooled nectar forager samples were collected for each pheromone treatment. Each bee diagram represents a sample, though two brains were used for each sample. Resulting numbers of reads per sample and percentages of those reads that mapped to the honey bee genome are presented in a table to the right Table 1 Numbers of DEG in all pairwise comparisons Upregulated Downregulated Pheromone Main Effect BP vs Control 12 46 EBO vs Control 14 138 Food Main Effect Pollen vs Nectar 79 246 Interaction Effect BP v Control and Food 29 39 EBO v Control and Food 55 32 Genes whose expression differed between groups were considered differentially expressed when they had a false discovery rate (FDR) of <0.05. Up- and downregulation of significantly differentially expressed genes was determined by whether log fold change was above or below zero, respectively Ma et al. BMC Genomics (2019) 20:592 Page 4 of 15
Ma et aL BMC Genomics 201920:592 Page 5 of 15 Table 2Overlaps between pheromone-related DEG tric tes and for ager type (Fig.2)significantly more often than random than nectar foragers in both princ expectation based on 10,000 ite ampling (P 上gur Overlaps with landmark studies d to EBO dustered with pollen foragers To explore the relationship between the results shown above and those of previous similar studies,we performed repre Whi EBO had than BP field et al 52]identified DEGs related to fo Pollen fora sed to BP or EBO were more similar to while Alaux et al.3).identificd DEGs related each other d of their e nd thos s emerged metric test p<005:Table 6)Thus eenes that wer To better understand the contributions of pheromone differentially expressed in the brains of nectar and pollen for treatment and forager the expre gers significantly wit genes inal com Similarly.we found a si nificant de sed of a linear combination of many genes. metrictest P)between DEGs associated with BPx a and the PC were use 15 days of continuous expos ing tasks o nectar and pollen foragers,indicating that the greatest lap significantly with short-term changes in brain expression axis of variation in gene regulation was rela to forager ated with the stimulatio n of foraging behavior type This i t with resu more DEGs associated with forager with pheromone exposure. the variance in the DEGs and Weighted gene co-expression network analysis(WGCNA) egan GCN to cons struct networks of genes based eir express PC2 seemed to separate bees exposed to control phero These module were independe of tra mone treatment from those exposed to BP,while sam- les fro be expose O were more y we looked ragers, est:pollen ys.nectar fo Table 3 Overlaps between pheromone-and foraging-related were significantly coreated to forager type,exposuret pheromone genes raging Genes Overlap BP 005:Fig or a comb orrelated with only one trait.Module 10 was the only 152 module that was ass ciated with all traits,while Module was a si ant overlap between pheror one-related DEGs and DEGs b was associated with forager type and EBO exposu verlap of genes than expected by chance:Pc0.001;
identical combinations of pheromone treatment and forager type (Fig. 2) significantly more often than random expectation based on 10,000 iterations of multiscale bootstrap resampling (P < 0.05; Additional file 4: Figure S1). Nectar foragers exposed to either BP or control pheromone treatments clustered together. However, nectar foragers exposed to EBO clustered with pollen foragers, suggesting that EBO exposure resulted in gene expression patterns of nectar foragers that were more similar to those of pollen foragers. This is consistent with the observation that EBO had a greater effect on overall gene expression than BP. Pollen foragers exposed to BP or EBO were more similar to each other than either group was to pollen foragers exposed to control treatments. Genes were also clustered based on the similarity of their expression, and several large clusters of genes emerged. To better understand the contributions of pheromone treatment and forager type on patterns of gene expression, we performed PCA on all DEGs with samples grouped by treatment. Each principal component (PC) was composed of a linear combination of many genes. Together, the first two PCs explained 63% of variance in the data, and the PCs were useful in separating samples by both pheromone treatment and forager type (Fig. 3). The first PC explained 46% of variance and separated nectar and pollen foragers, indicating that the greatest axis of variation in gene regulation was related to forager type. This is consistent with results from the differential gene expression analysis, which showed that there were more DEGs associated with forager type than with pheromone exposure. The second PC explained 17% of the variance in the DEGs and began to separate pheromone treatment from each other, although the separation was less distinct than for forager type. Specifically, PC2 seemed to separate bees exposed to control pheromone treatment from those exposed to BP, while samples from bees exposed to EBO were more intermediate. Pollen foragers, especially those exposed to EBO and control treatments, seemed to have a lower variance than nectar foragers in both principal components. PC3 and PC4 explained 14% and 5% of the variance in DEGs, respectively (Additional file 6: Figure S3). Overlaps with landmark studies To explore the relationship between the results shown above and those of previous similar studies, we performed representation factor analysis between our results and landmark honey bee transcriptome studies (Tables 5, 6) [36, 52]. Whitfield et al. [52] identified DEGs related to foraging ontogeny, while Alaux et al. [36]. identified DEGs related to long-term exposure to BP (i.e., primer pheromone effects). We found a significant overlap between the foraging-related DEGs identified in our study and those identified by [52] (hypergeometric test, P < 0.05; Table 6). Thus, genes that were differentially expressed in the brains of nectar and pollen foragers (our study) overlapped significantly with genes that were differentially expressed in nurses and foragers [52]. Similarly, we found a significant degree of overlap (hypergeometric test, P < 0.05) between DEGs associated with BP exposure in our study and BP-related DEGs identified in [36] after 15 days of continuous exposure. Thus, long-term changes in gene expression associated with impacts of BP exposure on the transition from nursing to foraging tasks overlap significantly with short-term changes in brain expression patterns associated with the stimulation of foraging behavior by BP. This ultimately suggests that behavioral plasticity utilizes common suites of genes at vastly different time scales. Weighted gene co-expression network analysis (WGCNA) We used WGCNA to construct networks of genes based solely on the similarity of their expression patterns to organize co-expressed genes into groups, called modules. These modules were constructed independently of trait information and were then correlated to traits using a generalized linear model. Specifically, we looked at relationships between each module and three traits of interest: pollen vs. nectar foraging, BP vs. control, and EBO vs. control. The WGCNA identified 16 modules that were significantly correlated to forager type, exposure to BP, exposure to EBO, or a combination thereof (GLM, P < 0.05; Fig. 4). Fourteen modules were significantly correlated with only one trait. Module 10 was the only module that was associated with all traits, while Module 16 was associated with forager type and EBO exposure, but not BP exposure. For each module, the most highly connected gene in the network was identified (Table 7), Table 2 Overlaps between pheromone-related DEG First Contrast Second Contrast DEGs in First Contrast DEGs in Second Contrast Overlap BP vs Control EBO vs Control 58 152 39* There was a significant overlap between BP-related DEGs and EBO-related DEGs in a hypergeometric test *significantly greater overlap of genes than expected by chance; P < 0.001; hypergeometric test Table 3 Overlaps between pheromone- and foraging-related DEG Pheromone genes Foraging Genes Overlap BP vs Control 58 386 41* EBO vs Control 152 386 71* There was a significant overlap between pheromone-related DEGs and DEGs related to foraging *significantly greater overlap of genes than expected by chance; P < 0.001; hypergeometric test Ma et al. BMC Genomics (2019) 20:592 Page 5 of 15
Ma et al BMC Genomics (201920:592 Page of 15 Table 4 Overlaps between pheromone-and foraging-related genes,separated by foraging preference Pheromone genes Pollen Upregulated Nectar Upregulated Overlap Pollen Overlap Nectar BP vs Control 58 246 40 EBO vs Control 46 0 ersac linear model c polen forgers and ownreged n re ulated in po antly greater werlap of genes than expected by chance:P<0.001:hypergeometric tes mo nways,ca To better understand the functions of the gene modules 05)for ger identified in this analysis,we performed phospholipid metabolism,neuroactiv od pathway n M path relation with food and like FoxO and AGE-RAGEian both brood pheromones,and modules 3 and 7 were se nent pathw ys like wnt signaling,and immune pathways lected based on their strong correlations with BP and like Toll and Imd signaling pathways(Wilcoxon,P<0.05). ood 10 Fig.2 Heatmap for the hie nectar were exposed to pheromone amples.Food an nple are represente ate pheromone treatments
providing a list of candidate genes. The top five most connected genes for each module can be found in the Additional file 7: Table S4. To better understand the functions of the gene modules identified in this analysis, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis on three modules (Table 8). Module 10 was chosen based on its significant correlation with food and both brood pheromones, and modules 3 and 7 were selected based on their strong correlations with BP and EBO, respectively. Module 10 was enriched for KEGG pathways related to metabolic pathways, carbon metabolism, fatty acid metabolism, and peroxisomes (Wilcoxon, P < 0.05). Module 7 was significantly enriched for glycerophospholipid metabolism, neuroactive ligand-receptor interaction, and hippo signaling pathway (Wilcoxon, P < 0.05). Module 3 was enriched for metabolic pathways like FoxO and AGE-RAGE signaling pathways, development pathways like wnt signaling, and immune pathways like Toll and lmd signaling pathways (Wilcoxon, P < 0.05). Table 4 Overlaps between pheromone- and foraging-related genes, separated by foraging preference Pheromone genes Pollen Upregulated Nectar Upregulated Overlap Pollen Overlap Nectar BP vs Control 58 79 246 1 40* EBO vs Control 152 79 246 0 71* Because foraging preference was a binary trait in the generalized linear model (i.e. either pollen or nectar), DEGs that were up-regulated in nectar foragers were by definition down-regulated in pollen foragers, and vice versa. Foraging-related DEGs were upregulated in pollen foragers (and downregulated in nectar foragers) when log fold change was greater than 0 and upregulated in nectar foragers (and downregulated in pollen foragers) when log fold change was less than zero. There was a significant overlap between pheromone-related genes and genes that were upregulated in nectar foragers *significantly greater overlap of genes than expected by chance; P<0.001; hypergeometric test Fig. 2 Heatmap for the hierarchical clustering of brain gene profiles. Honey bees foraging on pollen or nectar were exposed to pheromone treatments: Brood pheromone (BP), E-beta-ocimene (EBO), or a control. Rows correspond to differentially expressed genes, and columns represent samples. Food and pheromone treatments for each sample are represented between sample dendrogram and heatmap. The scale bar indicates variance stabilized gene expression values, with highly expressed genes in lighter colors and lower expression in darker colors. Clustering of samples shows two branches main branches, which correspond broadly to nectar foraging (left) and pollen foraging (right); however, nectar foragers exposed to EBO have expression profiles more similar to pollen foragers. Within pollen and nectar branches, there is also a split in pheromone treatments Ma et al. BMC Genomics (2019) 20:592 Page 6 of 15
Ma et aL BMC Genomics 2019)20:592 Page7 of 15 △A Pheromone 0 Food 0 10 PC1:46%variance or nectar forager-type.The percentageo Discussion pheromone communication and n the present study,weinv alization and s tim that d that (BP)and larval pheromones,brood pheromone and pollen foragers in the mushroom bodies of ocimene (EBO),woul a regu com era tha quantitative ion diff ntially ariation in the】 epending on task specialization.We found that nectar variation in the sion sugar concentration of nectar collected and the and larval ph en an brought back to the hive for necta Moreover,comparisons with previous studies suggest vs.pollen foraging 27,53,55-57.In our study,foraging that similar genes regulate the ontogeny of foraging be pecialization on nectar vs pollen foraging was a or and in gen xpre (wi and pollen foragers,which ounted for 46% of the pheromones affected transcriptic onal pathways more overall variation in DEGs (Fig 3).To elucidate transcrip trongly in nectar foragers than poll orage onal pathw respond t larval pher mone to pre that lar we (WGCNA)to ene nds ciated with traits of interest [58,59] specialization.Howeve ntified 16 genetic modules that were signifi esult that larva phero hich ith osur raging ialization (Fi Short exposure to both BP and EBO significantly al study begins to elucidate between tered gen expression protiles in the brains of foragers Table5 Overlaps be wee pheromon related genes and those of Alaux et al enes repre d in both after 5 days Alaux et al.B et al.136l used a micr ray to characterze br -pp d ment in Ala
Discussion In the present study, we investigated the genes and transcriptional pathways underlying rapid behavioral responses to pheromone signals in honey bee foragers specialized in pollen or nectar foraging. We hypothesized that two larval pheromones, brood pheromone (BP) and E-beta-ocimene (EBO), would regulate a common set of foraging genes in the brain, and that these pheromones would affect gene expression differentially depending on task specialization. We found that nectar and pollen foragers have distinguishable gene expression profiles, and that both larval pheromones do indeed regulate a shared set of genes and transcriptional pathways, supporting predictions 1 and 2, respectively. Moreover, comparisons with previous studies suggest that similar genes regulate the ontogeny of foraging behavior and foraging task specialization, and a common set of genes mediate both short- and long-term responses to BP, supporting prediction 3. However, larval pheromones affected transcriptional pathways more strongly in nectar foragers than pollen foragers, contrary to prediction 4. Therefore, we found support for the hypothesis that larval pheromones regulate a shared set of foraging genes in the brain, and that their effect depends on foraging preference or specialization. However, our results did not support our prediction that larval pheromone have a greater effect on pollen foragers. Instead, the data revealed that larval pheromones regulated genes are positively correlated with nectar foraging. Thus, our study begins to elucidate mechanistic links between larval pheromone communication and foraging specialization and suggests that common transcriptional pathways may regulate behavior across time scales. The present study demonstrates for the first time that there are transcriptional differences between nectar and pollen foragers in the mushroom bodies of honey bees (prediction 1). Several quantitative trait loci have been identified which underlie colony-level variation in the propensity to collect pollen vs. nectar, and these loci are associated with variation in the sugar concentration of nectar collected and the amount of pollen and nectar brought back to the hive [53, 54]. Previous studies have examined the genetic and behavioral differences associated with preference for nectar vs. pollen foraging [27, 53, 55–57]. In our study, foraging specialization on nectar vs. pollen foraging was associated with substantial differences in gene expression profiles (with almost 400 DEGs; Table 1), and with variation among nectar and pollen foragers, which accounted for 46% percent of the overall variation in DEGs (Fig. 3). To elucidate transcriptional pathways that respond to larval pheromones, we utilized weighted gene correlation network analysis (WGCNA) to provide a more detailed view of the molecular processes associated with traits of interest [58, 59]. WGCNA identified 16 genetic modules that were significantly correlated with foraging or pheromone exposure (Fig. 4), most of which were associated with foraging specialization (Fig. 4). Short exposure to both BP and EBO significantly altered gene expression profiles in the brains of foragers, Fig. 3 Principal component analysis of all DEG. The first two principal components (PCs) are displayed, together representing 63% of the total variation. Each point represents a single sample. PC1 separates samples based on food preference, whereas PC2 separates pheromone treatment, particularly for nectar foragers. Shape represents pheromone treatment. Color represents pollen or nectar forager-type. The percentage of variation in transcript expression patterns explained by each PC is shown in the axes Table 5 Overlaps between pheromone-related genes and those of Alaux et al Genes represented in both BPgenes Alaux et al., BP after 5 days Alaux et al., BP after 15 days Overlap BP5 Overlap BP15 BP vs Control 6039 49 104 85 1 2* Alaux et al. [36] used a microarray to characterize brain gene expression differences related to long-term exposure to BP (i.e. primer pheromone effects) at two time points, 5 days (BP5) and 15 days (BP15). Shown here are the results of a hypergeometric test between the DEGs related to BP in the present study and in Alaux et al. at each time point. There was a significant overlap between DEGs in the present study and the 15-day treatment in Alaux et al *significantly greater overlap of genes than expected by chance; P < 0.05; hypergeometric test Ma et al. BMC Genomics (2019) 20:592 Page 7 of 15
Ma et al BMC Genomics (201920:592 Table6Overlaps between foraging-related genes and Whitfield et al Genes represented in both Foraging-related BP vs Control 603g 839 tly gre than expected by and both phero sure had stron ene with 169 DEGs, which profiles.Furthe rmore,WGCNA revealed that module 10, than the numbe of DEGs epres 9 genes with cor expression patte rulated by BP an that RP a EBO(Table 2),and the overlapping genes were enriched EBO regulate overlapping genetic modules and pathways for fatty acid metabolism. Hierarchical clustering and that are enriched for energy metabolism.Decreasing whole 09 ME3 0.4 ME10 0.2 ME5 0 ME25 -0.2 ME28 ME22 P=04到 ME15 ME7 P-0.05 ME19 ME21 ME8 ME24 ME16 0.29 ME13 Fig.4 Weighted gene co t he colorized according correlation coefficient,varying from high values in yellow to low values in purple
and both pheromones regulated overlapping sets of genes (prediction 2). Exposure to EBO was associated with 169 DEGs, which was nearly three times greater than the number of DEGs regulated by BP (Table 1). Yet, even in this limited gene set, there was a statistically significant overlap in the DEGs regulated by BP and EBO (Table 2), and the overlapping genes were enriched for fatty acid metabolism. Hierarchical clustering and principal component analyses confirmed that pheromone exposure had strong and consistent effects on gene expression profiles. Furthermore, WGCNA revealed that module 10, representing 239 genes with correlated expression patterns, was significantly downregulated in samples exposed to either pheromone. Together, these results suggest that BP and EBO regulate overlapping genetic modules and pathways that are enriched for energy metabolism. Decreasing wholeTable 6 Overlaps between foraging-related genes and Whitfield et al Genes represented in both Foraging-related Whitfield et al Overlapping genes BP vs Control 6039 264 839 48* Whitfield et al. [52] identified DEGs related to foraging ontogeny during the transition between nurses and foragers, controlling for the effect of age. Presented here are the results of a hypergeometric test between the foraging-related DEGs in the present study and the DEGs identified in Whitfield et al., which show a significant overlap *significantly greater overlap of genes than expected by chance; P < 0. 05; hypergeometric test Fig. 4 Weighted gene co-expression network analysis. Rows represent gene modules. Columns represent sample traits. Each cell contains two values: a correlation coefficient between the module and sample trait and the associated p-value in parentheses. Significant correlations are colorized according correlation coefficient, varying from high values in yellow to low values in purple Ma et al. BMC Genomics (2019) 20:592 Page 8 of 15
Ma et aL BMC Genomics (2019)20:592 Page 9 of 15 Table7 WGCNA Module Hub genes Hub Gene Description A Forager-type.BP& B45943 BP Onhy GB42728 Sodium channel protein paralytic EBO Only 13 217 G4542 transmembrane protein 900 GB52595 zinc finger and BTB domain-containing protein 20 GB45063 LIMhomeobox protein Lhx9 24 267 GB19920 phosphopantothenoylcysteine decarboxylase 540 GB44289 ataxin3 Forager-type Only 19 127 GB50823 rotein kinase ATM 21 CB49517 DENN G851059 omains pro 1 GB45147 18 inding factor GB40539 405 ribosomal protein S20 25 89 GB51029 band 4.1-ike protein 5 that were o iffe ntially expressed in at least one contrast .that of fatty-acids s as e scale tasks to foragin asks 601.sugg sting that larval pherd nonal mones regulate for ging behavior by specifically activating this intriguing after hose of two andmark honey bee transcriptome stuc Table 8 KEGG analysis of selected WGCNA modules y5P<00 (EASE <005) on meta Fatty acid metabolisn BP alone Glycerophospholipid metabolism lon channel Neuroactive ligand-receptor interaction Hipoo signaling pathway EBO alone Pentose and clucuronate interconversion Intearal components of membrane Metabolic pathwavs FOXO signaling pathway Neuroactive ligand-rec Notch signaling pathway Toll and Imd signaling pathway Modle 10oncant corelation h ood and bo brood pheromones and modesandere selected aed on ther strong
brain energy metabolism, including that of fatty-acids, is associated with long-term behavioral transition from in-hive tasks to foraging tasks [60], suggesting that larval pheromones regulate foraging behavior by specifically activating pathways involved in the natural ontogeny of foraging behavior. Changes in the same behavior at different time scales, such as the ontogeny of pollen foraging and the pheromonal upregulation of pollen foraging, may utilize similar molecular mechanisms (prediction 3). We reached this intriguing conclusion after comparing our results to those of two landmark honey bee transcriptome studies Table 7 WGCNA Module Hub genes Regulation Pattern Module Size Hub Gene Hub Gene Description EBO & Forager-type 16 166 GB52658 Transcription factor All: Forager-type, BP, & EBO 10 239 GB45943a Collagen alpha-5 chain BP Only 7 560 GB42728 Sodium channel protein paralytic EBO Only 13 217 GB45423 transmembrane protein 3 900 GB52595 zinc finger and BTB domain-containing protein 20 8 90 GB45063a LIM/homeobox protein Lhx9 24 267 GB19920 phosphopantothenoylcysteine decarboxylase 6 540 GB44289 ataxin-3 Forager-type Only 19 127 GB50923 serine-protein kinase ATM 21 145 GB49517 DENN domain-containing protein 4C 22 121 GB51059 four and a half LIM domains protein 2 15 168 GB45147 a clavesin-2 26 82 GB41641a mitochondrial cardiolipin hydrolase 28 58 GB50931 box A-binding factor 5 594 GB40539 40S ribosomal protein S20 25 89 GB51029 band 4.1-like protein 5 a hub genes that were also differentially expressed in at least one contrast Table 8 KEGG analysis of selected WGCNA modules Module Trait association Significantly enriched KEGG pathways (P < 0.05) Significantly enriched GO categories (EASE <0.05) 10 BP, EBO, Foragertype Metabolic pathways Integral components of membrane, Fatty acid biosynthetic process Carbon metabolism Fatty acid metabolism Peroxisome 7 BP alone Glycerophospholipid metabolism Ion channel Neuroactive ligand-receptor interaction Hippo signaling pathway 3 EBO alone Pentose and glucuronate interconversions Integral components of membrane Metabolic pathways FOXO signaling pathway Neuroactive ligand-receptor interaction Lysosome Wnt signaling pathway Dorso-ventral axis formation Notch signaling pathway Toll and Imd signaling pathway AGE-RAGE signaling pathway in diabetic complications Module 10 was chosen based on its significant correlation with food and both brood pheromones, and modules 3 and 7 were selected based on their strong correlations with BP and EBO, respectively Ma et al. BMC Genomics (2019) 20:592 Page 9 of 15
Ma et al BMC Genomic (2019)20592 Page 10 of 15 [36,52].Whitfield et al.[52]compared nurses and for- ontogeny of foraging behavior [60].Pankiw et al. ers,cont or8eand exposure to increased poller on ge and found more ging to nollen forag which the authors found than 200 DEGs betweer age-matcl ed bees that were ex explanation may posed to usly for multiple days that short exposur an se tha t were notx are ake the ora tudies we comr ed 1)the DEGs betw e which may be a way to buffer aga inst ephemeral swing oragers in our sti ly with those identified by Whit t al emands of dev loping larvae 52.an the I A mo ated to b We e of Whitfield et al (P<001)and of Alau et al (P. but not epige netic pathways.which sug sts that metabol )The een our st and the sign may play that for s time scales (see 59)),and sup tore and changes in insulin signaling [61].Thes ports the idea hat ogical changes during the tion from in-hIv lar of fo references for nectar s genes foragin spec bu olism lipid signa ing path study tha 41 bel es been shown to elicit astrated the imp rtance of brain oli ollen agers.For exam ple, exp on influe individual variation in be old 42 ony-l enr and the individual e ollen for [44]How ng nath way in our ever,prior to this study,there were no docu nted im orts the role of sulin signaling pathways in mediating pacts of exposure to brood phe mo on ne inct and or in insects 64,65].For ex ommo t of DEC necialization (Table 3)which was driven primarily nd is related to FOXO signalin DEGs in nectar foragers but not pollen forag ers (Tabl Module3 was enriched for FOXO signaling and significantiy )Hierarchical clust ring analysis showe that,for the EBO treatment,so its hul genes may serve P s were nto po nt s nt he im on phe ectar foragers exp sed to EBO had expre ofiles ion and fora that were more like those of pollen foragers (Table 4) Although the results of this study are consistent with PC that ne rage exp erpre that ph tion may other n tar fora 3)The on which this in vely small revealed that two modules were associated with The consistency of the expression differences between our pheromon treat ent and foraging, one of which a tudy and previous studies 36 52],the patterns obtained ane compe ent rgy m he upe A g.hi stering ism hy which larval nheromone modulate colon indicate that our data ma reveal biologically meaningfu pollen foraging behavior could be by downregulating patterns despite the small sample size.However,future metabolic path role that energy metabolism plays in variation among foragers
[36, 52]. Whitfield et al. [52] compared nurses and foragers, controlling for age, and found over 1000 DEGs. Alaux et al. [36] were the first to study the effects of brood pheromone on gene expression, and found more than 200 DEGs between age-matched bees that were exposed to BP continuously for multiple days (i.e., five or 15 days) and those that were not exposed. To test the degree of overlap between our results and those from previous studies, we compared 1) the DEGs between nectar and pollen foragers in our study with those identified by Whitfield et al. [52], and 2) the DEGs between pheromone treatments in our study with those identified by Alaux et al. [36] . We found significant overlaps between DEGs identified in our results and those of Whitfield et al. (P < 0.001) and of Alaux et al. (P < 0.001). The significant overlap between our study and the two microarray studies, which validate the expression patterns related to foraging specialization and brood pheromone exposure, suggests that foraging-related gene expression shows a degree of consistency across time scales (see [59]), and supports the idea that pheromones regulate the transcriptional pathways underlying foraging specialization. Our data supported the hypothesis that exposure to larval pheromones alters expression of foraging related genes depending on foraging task specialization, but contrary to our prediction, the pheromones had more pronounced effects on gene expression in nectar foragers than pollen foragers (prediction 4). Larval pheromones have been shown to elicit specific responses in pollen foragers. For example, exposure to brood pheromone (BP) increased colony-level pollen foraging 2.5 fold [42], the ratio of pollen to non-pollen foragers [44], and the individual effort of pollen foragers [44]. However, prior to this study, there were no documented impacts of exposure to brood pheromones on nectar foraging. There was a common set of DEGs that were associated with both pheromone treatment and foraging specialization (Table 3), which was driven primarily by DEGs in nectar foragers but not pollen foragers (Table 4). Hierarchical clustering analysis showed that, for the most part, samples were clustered into pollen and nectar foraging “branches,” with the intriguing exception that nectar foragers exposed to EBO had expression profiles that were more like those of pollen foragers (Table 4). Similarly, PCA showed that nectar foragers exposed to EBO clustered more closely with pollen foragers than other nectar foragers (Fig. 3). The gene network analysis revealed that two modules were associated with both pheromone treatment and foraging, one of which was enriched for membrane components and energy metabolism (Table 8). These results suggest that one mechanism by which larval pheromones modulate colony-level pollen foraging behavior could be by downregulating metabolic pathways in the nectar forager brain, which is consistent with the role that energy metabolism plays in the ontogeny of foraging behavior [60]. Pankiw et al. [44] found that short exposure to BP increased pollen foraging, but did not observe task-switching of nectar foragers to pollen foraging, which the authors found puzzling. Our results indicate that one explanation may be that even after short exposures to larval pheromones, nectar foragers are primed to switch to pollen foraging even before they actually make the behavioral transition, which may be a way to buffer against ephemeral swings in the nutritional demands of developing larvae. DEGs and WGCNA modules related to both pheromone treatment and foraging specialization were enriched for several metabolic pathways, including fatty-acid metabolism, but not epigenetic pathways, which suggests that metabolic processes and lipid signaling in integration centers of the honey bee brain may play a role in behavioral plasticity. The transition from nursing to foraging involves large-scale changes in metabolic pathways, including reductions in lipid stores and changes in insulin signaling [61]. These physiological changes during the transition from in-hive tasks to foraging are associated with changes in energy metabolism (including insulin signaling), gustatory response, and foraging preferences for nectar vs pollen [62, 63]. Therefore, the prominence of energy metabolism, lipid signaling pathways, and related metabolic pathways in our study’s brain transcriptome data supports the idea that these pathways in the brain play a role in insect behavior [64, 65]. Other studies have demonstrated the importance of brain metabolic processes on influencing individual variation in behavior, particularly aggression [65–67]. The enrichment of metabolic pathways in DEGs and the prominence of the FOXO signaling pathway in our gene co-expression networks further supports the role of insulin signaling pathways in mediating neuronal function and behavior in insects [64, 65]. For example, an insulin binding protein, Queen brain-selective protein-1 (Qbp-1), was differentially expressed in response to pheromone treatment and is related to FOXO signaling. Module 3 was enriched for FOXO signaling and significantly correlated with EBO treatment, so its hub genes may serve as useful candidate genes for subsequent studies investigating the impact of insulin signaling on pheromone communication and foraging. Although the results of this study are consistent with the interpretation that pheromone communication may possibly regulate foraging task specialization, the sample sizes on which this interpretation rests are relatively small. The consistency of the expression differences between our study and previous studies [36, 52], the patterns obtained in the PCA, and the results of the unsupervised clustering strategies (WGCNA & hierarchical clustering) all serve to indicate that our data may reveal biologically meaningful patterns despite the small sample size. However, future studies will be required to assess whether any confounding factors—such as individual variation among foragers, Ma et al. BMC Genomics (2019) 20:592 Page 10 of 15
