The genomic signature of social interactions regulating honey bee caste development

Social evolution theory posits the existence of genes expressed in one individual that affect the traits and fitness of social partners. The archetypal example of reproductive altruism, honey bee reproductive caste, involves strict social regulation of larval caste fate by care-giving nurses. However, the contribution of nurse-expressed genes, which are prime socially-acting candidate genes, to the caste developmental program and to caste evolution remains mostly unknown. We experimentally induced new queen production by removing the current colony queen, and we used RNA sequencing to study the gene expression profiles of both developing larvae and their care-giving nurses before and after queen removal. By comparing the gene expression profiles between both queen-destined larvae and their nurses to worker-destined larvae and their nurses in queen-present and queen-absent conditions, we identified larval and nurse genes associated with larval caste development and with queen presence. Of 950 differentially-expressed genes associated with larval caste development, 82% were expressed in larvae and 18% were expressed in nurses. Behavioral and physiological evidence suggests that nurses may specialize in the short term feeding queen- versus worker-destined larvae. Estimated selection coefficients indicated that both nurse and larval genes associated with caste are rapidly evolving, especially those genes associated with worker development. Of the 1863 differentially-expressed genes associated with queen presence, 90% were expressed in nurses. Altogether, our results suggest that socially-acting genes play important roles in both the expression and evolution of socially-influenced traits like caste.


21
queen, and we used RNA sequencing to study the gene expression profiles of both developing larvae and their 22 care-giving nurses before and after queen removal. By comparing the gene expression profiles between both 23 queen-destined larvae and their nurses to worker-destined larvae and their nurses in queen-present and queen-24 absent conditions, we identified larval and nurse genes associated with larval caste development and with queen 25 presence. Of 950 differentially-expressed genes associated with larval caste development, 82% were expressed in 26 larvae and 18% were expressed in nurses. Behavioral and physiological evidence suggests that nurses may 27 specialize in the short term feeding queen-versus worker-destined larvae. Estimated selection coefficients 28 indicated that both nurse and larval genes associated with caste are rapidly evolving, especially those genes 29 associated with worker development. Of the 1863 differentially-expressed genes associated with queen presence,

42
The social insect sterile worker caste is the archetypal example of reproductive altruism that initially puzzled 43 Darwin (Darwin 1859) and spurred Hamilton (Hamilton 1964) to develop kin selection theory. Kin selection theory 44 presupposes the existence of genes that are expressed in one individual but have fitness effects on 45 relatives (Hamilton 1964). Despite this clear focus of social evolution theory on socially-acting genes, empirical 46 studies of the genetic basis of social insect traits, including caste, have widely overlooked the potential effects of 47 genes that are expressed in one individual but affect the traits of social partners.

48
Honey bee female caste is considered to be an exemplar polyphenism, whereby the expression of 49 alternate queen and worker morphs is controlled by environmental cues (Evans and Wheeler 1999). Unlike some 50 other well-studied polyphenisms that are controlled by simple abiotic factors such as temperature or 51 photoperiod (Nijhout 2003), honey bee queen-worker dimorphism critically depends on social control of larval 52 development by adult nestmates (Linksvayer et al. 2011). In vitro rearing studies demonstrate that in the absence 53 of social control, queen-worker dimorphism disappears and a continuous range of phenotypes are 54 produced (Linksvayer et al. 2011).

55
Honey bee colonies only rear new queens during specific life history stages, for example in the spring 56 when the colony is large enough to split in half, or upon the death of the current queen. Queen rearing is an 57 emergent, colony-level process involving the coordinated activities of hundreds or thousands of adult workers.

58
Necessary steps include the construction of special queen cells by nurse bees (Fig. 1), distinct provisioning behavior 59 of nurses coupled with distinct qualitative and quantitative differences in the nutrition fed to queen-and worker-60 destined larvae (colloquially known as "royal jelly" vs. "worker jelly") (Haydak 1970;Brouwers et al. 1987), the 61 larval developmental response to these environmental signals, and finally, selection by nurses of a subset of larvae 62 in queen cells to be reared to adulthood (Hatch et al. 1999

73
As a result, the contribution of genes expressed in adult nestmates (e.g., nurses and foragers) to the 74 genetic basis and evolution of the honey bee caste developmental program has received relatively little attention.

85
We use a novel extension of the conventional genetic approach to begin to characterize the full set of 86 molecular interactions underlying social interactions that regulate reproductive caste development (Linksvayer et 87 al. 2012), or the "social interactome" of caste. Instead of searching only for associations between an individual's 88 own patterns of gene expression and its traits, we also search for associations between social partners' patterns of 89 gene expression and the traits of focal individuals. Specifically, we used RNA sequencing of queen-and worker-90 destined larvae as well as "royal nurses" and "worker nurses", nurses collected in the act of feeding queen-and 91 worker-destined larvae, respectively. Subsequently, we used a new honey bee population genomic dataset (Harpur 92 et al. 2014) to study how rates of molecular evolution vary at the genes we identified as being associated with the 93 caste developmental program. We also determined whether there was evidence for behavioral and physiological 94 specialization of nurses to feed queen-versus worker-destined larvae, as such specialization is expected to 95 strengthen the transcriptional signature of social effects on caste development.

99
Basic setup, honey bee behavioral observations, and sampling

108
The first study used two replicate observation hives. Every three days beginning 24 days before the start 109 of the study, we individually marked 400 newly emerged adult workers with a unique combination of numbered 110 tag glued onto the mesosoma and an age-specific abdomen paint mark, and we added 200 individually marked 111 workers to each observation hive. Frames of known-aged brood were produced by caging queens on empty frames 112 for 24 hours and then checking for the presence of eggs. Four days later, one frame with only similarly-aged 1 st 113 instar larvae was placed into each observation hive, and the queen was removed to initiate emergency queen 114 rearing. These frames were the source of young focal larvae, a fraction of which were reared as new queens, and 115 the rest as workers. Within the first two days of queen removal, nurse workers build wax queen cells over young 116 focal brood and begin provisioning these queen-destined larvae differentially than worker-destined larvae in 117 worker cells (Figure 1). We continually observed areas of the frame with focal brood that contained both queen 118 cells and worker cells and recorded the date, time, and identity of nurses observed provisioning queen or worker 119 cells (i.e. "royal nurses" or "worker nurses"). Feeding behavior was defined when workers had their head 120 positioned deep enough into the worker or queen cell to be in contact with the larva, and remained motionless 121 except for a rhythmic motion of the abdomen for at least 5 seconds.

122
The second study used three replicate observation hives. The setup followed the first study, except that 123 we collected samples of focal brood under both queen-present and queen-removed conditions. First, on the fourth 124 day after introducing focal brood, samples of five 4 th instar worker-destined larvae and 20 nurses observed feeding 125 4 th instar worker-destined larvae were collected. Two days later, a new frame of same-aged 1 st instar larvae was 126 added to each of the three observation hives, and each colony queen was removed in order to initiate emergency 127 queen rearing. On the fourth day after introducing focal brood and removing the queen, we collected five 4 th instar 128 worker-destined larvae from the frame of focal brood, and we collected 20 worker nurses in the act of feeding

137
(3) queen larvae from queenless colonies, (4) worker nurses from colonies with a queen, (5) worker nurses from 138 queenless colonies, and (6) royal nurses from queenless colonies. Thus, for both larvae and nurses, there were 139 three conditions that were associated with: the production of new workers in hives with a queen; the production 140 of new workers in queenless hives; and the production of new queens in queenless hives. We extracted RNA from 141 whole larvae, but from nurses we dissected the two main glandular sources of proteinaceous brood food, the (maxround=40) to ensure convergence. With DESeq2 we used default settings and ran two separate analyses to 179 identify genes with differential expression associated with queen-vs. worker-development and genes with 180 differential expression associated with queen presence vs. absence. We focus on the EBSeq results for subsequent 181 analyses because EBSeq is most appropriate for our study, but we also report DESeq2 results because the DEseq2 182 analysis was more conservative for identifying genes associated with caste ( Fig. S1), but not for genes associated 183 with queen presence (Fig. S2). Subsequent analyses were qualitatively similar following either EBSeq and DESeq2

218
Differential expression analysis 219 220 We identified 950 differentially expressed genes associated with whether larvae developed into new queens or 221 workers (Table 1; Table S2). As expected, the majority of these genes (82%; 779/950) were differentially expressed 222 in the larvae themselves, depending on whether the larvae were queen-or worker-destined larvae. 18% (171/950) 223 were differentially expressed in nurses collected while feeding queen-destined larvae compared to nurses 224 collected while feeding worker-destined larvae (3 expressed in MG, 105 H, and 63 HPG) ( Table S2). Overlap of 225 differentially-expressed genes associated with caste development is shown by tissue type in Figure S3.

226
We also identified 2069 genes that were differentially expressed depending on queen presence, i.e. 227 whether the mother queen was present or removed, irrespective of larval caste fate or nurse behavior (Table 2;   228   Table S3). 90% (1863/2069) were expressed in nurse tissues, especially MG (1744 MG, 105 H, 15 HPG), and 206 229 were expressed in larval tissue. Overlap of differentially-expressed genes associated with queen presence is shown 230 by tissue type in Figure S4.

231
Considering the top 25 most highly expressed genes for each tissue (Table S1) (Table S1). Approximately 234 one third of each set of most highly-expressed genes was unique to each nurse tissue, whereas ~90% (22/25) of 235 the most highly expressed larval genes were unique to larvae (Fig. S5).

236
GO enrichment analysis for differentially-expressed genes associated with caste or queen presence are 237 shown by tissue type in Tables S4 and S5, respectively. Among genes associated with caste: genes differentially 238 expressed in nurse HPG tissue were enriched for GO terms associated with translation and several categories 239 associated with immune function; nurse head tissue genes showed a weaker signal of enrichment for a range of 240 GO terms, including signaling; and larval-expressed genes were enriched for terms such as metabolic processes 241 and chromatin assembly. Among genes that were differentially expressed depending on queen presence: nurse 242 MG genes were enriched for a range of terms including translation and transcription, macromolecular 243 biosynthesis, signal transduction, metabolism, and immune response; nurse head tissue genes were enriched for 244 immune system function, brain development, and chromatin assembly; and larval genes were enriched for terms 245 such as response to oxidative stress and metabolism.

249
Differentially-expressed genes, whether associated with caste development or queen presence had higher average 250 selection coefficients (γ) than non-differentially expressed genes ( Fig. 4; glm on log-transformed gamma estimates, 251 6 all P < 10 -8 ), and furthermore, genes with expression associated with caste or both caste and queen presence had 252 higher γ than genes with expression only associated with queen presence (Fig. 4; Tukey contrasts, both P < 10 -4 ).

253
Next we focused on genes with caste-associated expression. To further compare patterns of molecular 254 evolution at genes associated with queen vs. worker production, we defined genes up-regulated in queen larvae or 255 royal nurse tissues as "queen-associated genes" and genes up-regulated in worker larvae or worker nurse tissues 256 as "worker-associated genes". Mean γ for worker-associated genes was higher than queen-associated genes (glm 257 with log-link on γ + 2 values, t = 2.47, df = 824, P = 0.014), and did not depend on whether the genes were 258 expressed in larval or nurse tissues (P = 0.33) (Fig. 5). When only considering nurse-expressed genes, γ was higher 259 for queen-associated vs. worker-associated genes (t = 3.71, df = 135, P = 0.0076), but γ was not significantly 260 different when only considering larval-expressed genes (t = 1.78, df = 688, P = 0.076).

264
We simultaneously studied the gene expression profiles of two classes of socially-interacting individuals --

303
Genes with caste-associated expression had higher estimated selection coefficients than non-304 differentially-expressed genes and genes with expression dependent on queen presence (Fig. 4). Furthermore, 305 among caste-associated genes, genes up-regulated in worker larvae and worker nurses had higher selection 306 coefficients than genes up-regulated in queen larvae and royal nurses (Fig. 5). Altogether, these results indicate 307 that both genes with direct effects and putative indirect effects on larval development -especially those 308 associated with worker development -may have experienced positive selection and contributed to the evolution 309 of the honey bee caste system. Our results fit with two recent honey bee studies showing that genes associated 310 with adult worker traits are also rapidly evolving. The first study shows that genes encoding proteins that are more

325
As expected, genes whose proteins make up the primary components of royal jelly, including 8 of the 9 326 major royal jelly proteins (MRJPs), were among the most highly-expressed genes in nurse tissues (Table S1) and 327 were also differentially expressed (Tables 1 and 2). However, of the mrjp genes, only the expression of mrjp3, 328 which has previously been implicated as promoting queen development (Huang et al. 2012), depended on nurse 329 behavior: it was up-regulated in the head tissue of royal nurses (Table 1). All eight differentially-expressed mrjp 330 genes, including mrjp1, also implicated as promoting queen development(Kamakura 2011), were differentially 331 expressed in nurse mandibular glands or head tissue, depending on queen presence. Most were up-regulated in 332 the queen-removed condition (Table 2), presumably related to colony-level changes associated with the rapid shift 333 to emergency queen rearing.

334
Notably, 4 of the 6 described honey bee antimicrobial peptides (Evans et al. 2006) (defensin 1, abaecin, 335 hymenoptaecin, and apisimin) were up-regulated in the HPG and / or the head tissues of nurses feeding queen-336 destined larvae (Table 1) and caste-associated nurse-expressed genes were enriched for Gene Ontology terms for 337 immune function (Table S4). Interestingly, hymenoptaecin and another antimicrobial peptide, apidaecin were also 338 upregulated in queen-destined larvae. Altogether, these results suggest that queen-and worker-destined larvae 339 may require different levels of antimicrobial peptides, some of which may be produced by nurse workers and 340 transferred to larvae through royal jelly (Schonleben et al. 2007;Furusawa et al. 2008;Zhang et al. 2014).

341
Besides broad differences in gene expression, royal nurses also had larger HPG acini and were 1.5 days 342 younger than worker nurses ( Fig. 2 and 3). Previous studies have shown that nurse HPG gland size and 343 activity (Ohashi et al. 2000), as well as the composition of nurse glandular secretions (Haydak 1970), and patterns of 344 nurse brain gene expression (Whitfield et al. 2006) all vary with nurse age and social environment. While it is not 345 clear how exactly these differences are related to the observed differences in nurse provisioning behavior, 346 individually-marked nurses did tend to specialize on feeding either queen or worker cells within a day, but not 347 across multiple days. Thus, individual nurses may be physiologically primed to contribute to queen rearing only on 348 the short term. Longer-term tracking of individuals during queen rearing will be necessary to definitively 349 demonstrate the degree to which nurse specialization occurs. The key point for this study of colony-level caste 350 regulation is that queen-vs. worker-destined larvae interact with nurses that are on average transcriptionally and 351 physiologically distinct, resulting in distinct rearing environments and alternate caste developmental trajectories.