Honey bee (Apis mellifera) larval pheromones regulate gene expression related to foraging task specialization

Background Foraging behavior in honey bees (Apis mellifera) is a complex phenotype which 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 behavior 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, but the mechanisms underpinning this rapid behavioral plasticity are unknown. Furthermore, the mechanisms through which foragers specialize on collecting nectar or pollen, 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). Results We hypothesized that both pheromones would alter expression of genes in the brain related to foraging and would differentially impact expression of genes in the brains of pollen compared to nectar foragers. 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 Social signals (i.e. pheromones) may invoke foraging-related genes to upregulate pollen foraging at both long and short time scales. These results provide new insights into how social signals integrate with task specialization at the molecular level and highlights the important role that brain gene expression plays in behavioral plasticity across time scales.

were enriched for lipid biosynthesis and metabolism (FDR <0.05). 160 The DEGs associated with either EBO or BP were also analyzed separately. Because 161 there were few upregulated genes associated with either pheromone, up-and down-regulated 162 genes for each pheromone were pooled during pathway enrichment analysis; however, it should 163 be noted that the results for pheromone could potentially be driven by down-regulated genes. 164 DEGs associated with BP exposure were enriched for lipid biosynthesis and integral components 165 of the membrane (FDR < 0.05). DEGs associated to EBO exposure were enriched for integral 166 components of membrane, fatty acid biosynthetic processes, fatty acid metabolism, and pentose 167 phosphate pathway. There was a significant overlap of 39 genes between BP and EBO exposed 168 foragers compared to controls (P<0.05), and these DEGs were significantly enriched for 169 metabolic pathways and fatty acid metabolism (FDR<0.05). significantly more often than random expectation based on 10,000 iterations of multiscale 175 bootstrap resampling (P<0.05; Supplementary Fig. 1). Nectar foragers exposed to either BP or 176 control pheromone treatments clustered together; however, nectar foragers exposed to EBO 177 clustered with pollen foragers, suggesting that EBO exposure resulted in gene expression 178 patterns of nectar foragers that were more similar to those of pollen foragers. This is consistent 179 with the observation that EBO had a greater effect on gene expression than BP. Pollen foragers 180 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 182 expression, and several large clusters of genes emerged. 183 To better understand the contributions of pheromone treatment and forager type on 184 patterns of gene expression, we performed PCA on all differentially expressed genes with 185 samples grouped by treatment. Each principal component (PC) was composed of a linear 186 combination of many genes. Together, the first two PCs explained 63% of variance in the data, 187 and the PCs were useful in separating samples by both pheromone treatment and forager type 188 (Fig. 3). The first PC explained 46 % of variance and separates nectar and pollen foragers, 189 indicating that the greatest axis of variation in gene regulation related to forager type. This is 190 consistent with results from the differential gene expression analysis, which showed that there 191 were more DEGs associated with forager type than with pheromone exposure. Nectar foragers and began to separate pheromone treatment from each other, although the separation was less 195 distinct than for forager type. Specifically, PC2 seemed to separate bees exposed to control 196 pheromone treatment from those exposed to BP, while samples from bees exposed to EBO were 197 more intermediate. Pollen foragers, especially those exposed to EBO and control treatments, 198 seemed to have a lower variance than nectar foragers in both principal components. 199 Overlaps with landmark studies. To explore the relationship between the results shown 200 above and those of previous similar studies, we performed representation factor analysis between 201 our results and landmark honey bee transcriptome studies (Tables 5, 6 nursing to foraging tasks overlap significantly with short-term changes in brain expression 212 patterns associated with the stimulation of foraging behavior by BP. This ultimately suggests 213 that behavioral plasticity to utilize common suites of genes at vastly different time scales seems. 214 Weighted gene co-expression network analysis. We used WGCNA to construct networks 215 of genes based solely on the similarity of their expression patterns to organize co-expressed 216 genes into groups, called modules. These modules were constructed independently of trait 217 information and were then correlated to traits using a generalized linear model. Specifically, we 218 looked at relationships between each module and three traits of interest: pollen vs. nectar 219 foraging, BP vs. control, and EBO vs. control. In this way, the WGCNA identified 16 modules 220 that were significantly correlated to forager type, exposure to BP, exposure to EBO, or a 221 combination thereof (GLM, P<0.05; Fig. 4). Fourteen modules were significantly correlated with 222 only one trait. Module 10 was the only module that was associated with all traits, while Module 223 light cyan was associated with forager type and EBO exposure, but not BP exposure. For each 224 module, the most highly connected gene in the network was identified (Table 7), providing a list 225 of candidate genes. The top five most connected genes for each module can be found in the 226

Supplementary Materials. 227
To better understand the functions of the gene modules identified in this analysis, we 228 performed KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis on three 229 modules (Table 8). Module 10 was chosen based on its significant correlation with food and both 230 brood pheromones, and modules 3 and 7 were selected based on their strong correlations with BP 231 and EBO, respectively. Module 10 was enriched for KEGG pathways related to metabolic 232 pathways, carbon metabolism, fatty acid metabolism, and peroxisomes (Wilcoxon, P<0.05). 233 In the present study, we investigated the genes and transcriptional pathways underlying 240 rapid behavioral responses to pheromone signals in honey bee foragers specializing on pollen or 241 nectar foraging. We hypothesized that two larval pheromones, brood pheromone (BP) and E-242 beta-ocimene (EBO), would regulate a common set of foraging genes in the brain, and that they 243 would affect gene expression more strongly in pollen foragers than nectar foragers. We found 244 that nectar and pollen foragers have distinguishable gene expression profiles, and that both larval 245 pheromones do indeed regulate a shared set of genes and transcriptional pathways; however, 246 these transcriptional pathways more strongly affect nectar foragers than pollen foragers.  The present study demonstrates for the first time that there are also transcriptional 264 differences between nectar and pollen foragers in the mushroom bodies of honey bees. Several 265 quantitative trait loci have been identified which underlie colony-level variation in the propensity 266 to collect pollen versus nectar, and these loci are associated with variation in concentration of 267 nectar collected and the amount of pollen and nectar brought back to the hive [54, 56]. In our 268 study, foraging specialization on nectar versus pollen foraging was associated with substantial 269 differences in gene expression profiles (with almost 400 DEGs; Table 1), and with variation 270 among nectar and pollen foragers, which accounted for 46% percent of the overall variation in 271 DEGs (Fig. 3). To elucidate transcriptional pathways that respond to larval pheromones, we 272 utilized weighted gene correlation network analysis (WGNCA) to provide a more detailed view of molecular processes associated with traits of interest [57,58]. WGNCA identified 16 genetic 274 modules that were significantly correlated with foraging or pheromone exposure (Fig. 4), most of 275 which were associated with foraging specialization (Fig. 4). 276 Short exposure to both BP and EBO significantly altered gene expression profiles in the 277 brains of foragers, and both pheromones regulated overlapping sets of genes. Exposure to EBO 278 was associated with 169 DEGs, which was nearly three times greater than the number of DEGs 279 regulated by BP (Table 1). Yet, even in this limited gene set, there was a statistically significant 280 overlap in the DEGs regulated by BP and EBO (Table 2), and the overlapping genes were 281 enriched for fatty acid metabolism. Hierarchical clustering and principal component analyses 282 confirmed that pheromone exposure had strong and consistent effects on gene expression 283 profiles. Furthermore, WGCNA revealed that module 10, representing 239 genes with correlated 284 expression patterns, was significantly downregulated in samples exposed to either pheromone. 285 Together, these results suggest that BP and EBO regulate overlapping genetic modules and 286 pathways that are enriched for energy metabolism. Decreasing whole-brain energy metabolism, 287 including that of fatty-acids, is associated with long-term behavioral transition from in-hive tasks 288 to foraging tasks, suggesting that larval pheromones regulate foraging behavior by specifically 289 activating pathways involved the natural ontogeny for foraging behavior [59]. 290 Our data supported our hypothesis that exposure to larval pheromones alters expression 291 of foraging related genes, but contrary to our predictions, the pheromones had more pronounced 292 effects on gene expression in nectar foragers than pollen foragers. There was a common set of 293 DEGs that were associated with both pheromone treatment and foraging specialization (Table 3), 294 which was driven primarily by DEGs in nectar foragers but not pollen foragers (Table 4). 295 Hierarchical clustering analysis show that, for the most part, samples were clustered into pollen and nectar foraging "branches," with the intriguing exception that nectar foragers exposed to 297 EBO had expression profiles that were more like those of pollen foragers (Table 4). Similarly, 298 principal components analysis show that nectar foragers exposed to EBO cluster more closely 299 with pollen foragers than other nectar foragers (Fig 3). The gene network analysis revealed that 300 two modules were associated with both pheromone treatment and foraging, one of which was 301 enriched for membrane components and energy metabolism (Table 8). These results suggest that 302 one mechanism by which larval pheromones modulate colony-level pollen foraging behavior is 303 by downregulating metabolic pathways in the brains of nectar foragers, which is consistent with 304 the role that energy metabolism plays in the ontogeny of foraging behavior [59]. Pankiw et al. 305 [44] found that short exposure to BP increased pollen foraging, but did not observe task-306 switching of nectar foragers to pollen foraging, which the authors found puzzling. Our results 307 indicate that one explanation may be that even after short exposures to larval pheromones, nectar 308 foragers may be primed to switch to pollen foraging even before they actually make the 309 behavioral transition, which may be a way to buffer against ephemeral swings in the nutritional 310 demands of developing larvae. 311 Changes in the same behavior at different time scales, such as ontogeny of pollen 312 foraging and pheromonal upregulation of pollen foraging, may utilize similar molecular 313 mechanisms. We reached this intriguing conclusion after comparing our results to those of two 314 the effects of brood pheromone on gene expression, and found more than 200 DEGs between 317 age-matched bees that were exposed to BP continuously for multiple days (i.e., 5 or 15 days) and 318 those that were not exposed. To test the degree of overlap between our results and those from 319 previous studies, we compared 1) the number of DEG between nectar and pollen foragers in our 320 study with those identified by Whitfield  The results of our study lay the groundwork for several intriguing lines of inquiry for 329 future studies. First, exposing foragers to a short pulse of BP, which stimulates immediate 330 foraging, regulates a similar set of genes as bees as exposing nurses to BP for 5 days, which 331 modulates the transition from in-hive tasks to foraging. BP potentially regulates foraging 332 behavior, at least in part, by priming their to receptivity to foraging-related or social stimuli, even 333 before nurse bees have made the physiological transition to foraging tasks. This could 334 conceivably involve genes implicated in both foraging and division of labor (e.g., Malvolio, a 335 manganese transporter) [60], or neurochemical regulatory pathways involving octopomine, 336 which has shown to modulate responsiveness to both foraging-stimuli and to BP [61, 62]. 337 Furthermore, our results suggest the hypothesis that social pheromones upregulate pollen 338 foraging by decreasing the expression of nectar foraging genes in the brain, and this would also 339 be productive line of inquiry for future studies. Lastly, our data suggest that rapid changes in 340 brain gene expression in nectar foragers may happen prior to task-switching to buffer against 341 ephemeral environmental conditions. Short-term exposure to larval pheromones may "prime" protein, Queen brain-selective protein-1 (Qbp-1), was differentially expressed in response to 362 pheromone treatment and is related to FOXO signaling. Module 3 was enriched for FOXO 363 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 365 pheromone communication and foraging. 366

Conclusions 367
The neural circuits and molecular pathways underlying behavioral plasticity and task 368 specialization are complex, and our study demonstrates that foraging behavior is regulated in part 369 by common suites of genes across time scales, from long-term behavioral plasticity (nurse to 370 forager) to individual variation in task specialization (pollen vs nectar). Our study further 371 confirms that pheromone communication has a profound effect on gene expression within hours 372 of exposure, and more importantly, that social signals (i.e. pheromones) may invoke foraging-373 related transcriptional pathways to upregulate pollen foraging at both long and short-time time 374 scales. Moreover, there is clear interaction between individual variation in task specialization and 375 responses to social signals (i.e. pheromones), and these social signals seem to be manipulating 376 brain energy metabolism to elicit foraging behavior. Because the mechanisms underlying 377 foraging behavior are deeply conserved in animals, a detailed mechanistic understanding of the 378 foraging in honey bees may provide insights into mechanisms involved in division of labor in 379 insect societies, foraging preference in animals, and complex behavioral phenotypes in general. 380

Animals and Experimental Design 382
We created single-cohort colonies (using same-aged workers) from a common source colony 383 with a naturally mated queen to avoid differences in behavior and gene expression due to 384 variation in age and genetic background; thus, all bees used in the study were half-sisters. 385 Because workers performing any given task in a natural colony can vary widely in age, we 386 constructed single cohort colonies using workers that emerged as adults within a 48-hr period, minimizing differences in age and experience among individuals. After one week, some of the 388 bees in single-cohort colonies transition quickly to foraging (Robinson et al. 1989). feeders. Only bees with multiple paint marks of a single color were used in this study because 398 multiple same-color marks demonstrated consistent preference for pollen or nectar, respectively. 399 Each colony was also provided a strip of queen pheromone (PseudoQueen, Contech Industries, 400 Victoria, BC, Canada) to prevent colonies from developing a "queenless" physiological state and 401 to control for the variation in queen quality that would inherently occur when using live queens. 402 Colonies did not receive any frames containing brood. Although broodless colonies are not the 403 default state of a colony, it is nevertheless a natural occurrence when queens die for any number 404 of reasons [52]. The absence of brood controlled for the natural variation in brood pheromones 405 that may have occurred with the presence of real brood, and minimized the amount of 406 beekeeping interference necessary to maintain identical colony conditions. 407 After two weeks, colonies were exposed to field-relevant dosages (5,000 larval 408 equivalents) of EBO, BP, or a paraffin oil control. We used a BP blend characteristic of older 409 larvae that has been shown to strongly upregulate pollen foraging, as done previously by Ma et 410 al. [45]. Hive entrances were blocked during the one-hour period during which the pheromone 411 treatment was applied, and any foragers outside the colony during that time were removed from 412 the experiment when they landed on the blocked entrance. When the entrances were opened, 413 forager bees (previously marked as nectar or pollen foragers, as described above) were collected 414 as they landed on a pollen or nectar feeder, but before they initiated feeding. Six pollen foragers 415 and six nectar foragers were collected from each colony and placed immediately in dry ice for 416 later brain dissection. Subsequently, we pooled pairs of bees to generate the RNA samples. 417 Thus, in total, there were three pollen forager and three nectar forager samples for each of the 418 three colonies representing the control, BP and EBO treatments ( Fig. 1; total number of samples 419 = 18). Sampled individuals were stored at -80°C until they were dissected. 420

Brain Dissection 421
Insect mushroom bodies are considered important integration centers of the brain because 422 of their role in multimodal information processing and their association with learning and 423 memory [70-72]. These factors make mushroom bodies an ideal candidate brain region to 424 investigate temporal dynamics of communication and behavior. Therefore, mushroom bodies of 425 the brain were dissected from sampled individuals by placing them on dry ice to prevent thawing 426 and degradation of transcripts, as in [37,40]. However, in our study, the brains were not freeze-427 dried to facilitate dissection of the mushroom bodies. For each sample, RNA from two brains 428 were extracted using the RNeasy Mini Kit (Qiagen, Valencia, CA) following the manufacturer's 429 protocol, and pooled RNA quantity and quality were assayed using a Qubit Fluorometer 430 (Invitrogen, Carlsbad, CA). cDNA library preparation and sequencing were performed by the 431 Genomic Sequencing and Analysis Facility at the University of Texas at Austin using an 432 split across four lanes to control for sequencing bias. A total of 18 RNA-seq single-end 50bp 434 libraries were generated, with three libraries for each treatment group from each colony (Fig 1). 435

Data Analysis 436
Reads were trimmed using Trimmomatic [73] to remove adapter sequences, low-quality 437 reads, and short reads (<36 bp). Kallisto software was used to build a transcriptome index based 438 on a recently updated honey bee reference genome annotation (Amel_HAv3.1), and 439 subsequently to quantify the abundance of transcripts represented in each sample [74]. The R 440 package tximport was then used to import transcript abundances and generate a gene-level count 441 matrix that was scaled to both library size and transcript length [75]. The transcript-gene 442 correspondence was derived from the genome annotation using the R package rtracklayer [76]. 443 The R package DESeq2 [77] was used to collapse technical replicates (i.e. identical samples 444 across multiple sequencing lanes) and perform differential gene expression analysis with 445 pheromone treatment, forager type, and the interaction of pheromone treatment and forager type 446 as fixed effects in a generalized linear model. Only genes with an abundance of at least one 447 transcript per million in all samples were used in the analysis. Genes whose expression differed 448 between groups were considered differentially expressed when they had a false discovery rate 449 (FDR) of <0.05. Gene ontology analysis was performed using the Database for Annotation, 450 Visualization and Integrated Discovery (DAVID v6.8) to better understand the biological 451 relevance of differentially expressed genes (DEGs) [78]. 452 The expression patterns of DEGs were further analyzed by performing unsupervised 453 hierarchical clustering (Fig. 2; Supplementary Fig 1) and PCA (Fig. 3)  Variance stabilized gene expression data were grouped into modules based on similarity 471 of expression patterns. Because genes within each module showed very highly correlated 472 patterns, the first principal component of the genes within a module was used to represent the 473 entire module (module 'eigengene'). Then, these module representatives were correlated with 474 sample traits using a generalized linear model, with forager-type and pheromone as fixed effects 475 (Fig. 4). Minimum module size was set to 30, and deep split was set to 2. Modules were built 476 with a standardized connectivity score of -2.5, and module definition was based on "hybrid" 477 branch cutting. A signed gene co-expression network was constructed with a soft threshold of 10. 478 sample traits using a generalized linear model with forager-type, pheromone exposure treatment, generated from nectar and pollen foragers exposed to three pheromone treatments. Three pooled 539 pollen forager samples and three pooled nectar forager samples were collected for each 540 pheromone treatment. Each bee diagram represents a sample, though two brains per used for 541 each sample. Resulting numbers of reads per sample and percentages of those reads that mapped 542 to honey bee genome are presented in a table to the right. 543 foraging on pollen or nectar were exposed to pheromone treatments: Brood pheromone (BP), E-545 beta-ocimene (EBO), or a control. Rows correspond to differentially expressed genes, and 546 columns represent samples. Food and pheromone treatments for each sample are represented 547 between sample dendrogram and heatmap. The scale bar indicates z-scores of variance stabilized 548 gene expression values, with highly expressed genes in lighter colors and lower expression in 549 darker colors. Clustering of samples shows two branches main branches, which correspond 550 broadly to nectar foraging (left) and pollen foraging(right); however, nectar foragers exposed to 551 EBO have expression profiles more similar to pollen foragers. Within pollen and nectar 552 branches, there is also a split in pheromone treatments. 553