Differential expression analysis
We analyzed differential gene expression as a function of aggression on a per-tissue basis. 85, 1571, and 312 genes were differentially expressed in the brain, fat body, and midgut tissues, respectively (Tables S1-S3). Genes in the brain were significantly biased towards up-regulation in low aggression bees (81%, binomial test, P<0.0001), while direction of expression was not significantly biased in the fat body (49% upregulated, binomial test, P=0.27) or midgut (55%, binomial test, P=0.07).
To describe the function of genes differentially expressed as a function of aggression, we performed a Gene Ontology (GO) analysis followed by a REViGO analysis of significant GO terms (Benjamini-Hochberg corrected P < 0.05). REViGO clusters GO terms on the basis of semantic similarity to identify major patterns in long GO term lists (Supek et al. 2011). Differentially expressed genes in the brain were significantly enriched for 23 GO terms (Table S4). The REViGO clustering analysis showed clusters of processes and functions related to chaeta morphogenesis, disaccharide transport, and RNA polymerase II regulatory region sequence-specific DNA binding. These results suggest strong roles for transcriptional regulation, sensory development, and carbohydrate metabolism in differentiating brain gene expression profiles for high versus low aggression bees. Differentially expressed fat body genes were significantly enriched for 188 terms (Table S5), including processes and functions associated with nucleotide and energy metabolism, and transporter activity. Only one GO category, toxin activity, was significantly enriched among differentially expressed midgut genes.
All pairwise tissue comparisons showed some overlap in terms of genes differentially expressed as a function of aggression, with the strongest similarities between the midgut and fat body. Eight genes were differentially expressed in both the fat body and brain (enrichment test for significant overlap, P=0.79), and seven of eight genes showed the same direction of change as a function of aggression (binomial test, P=0.07). For the brain and midgut, six genes overlapped (P=0.006) with five of six genes showing the same direction of change (binomial test, P=0.22). 76 genes overlapped between the fat body and midgut (hypergeometric test, P<0.0001), with 71 showing the same direction of regulation across these two tissues (binomial test, P<0.0001), suggesting robust tissue expression similarity across these tissues. Only a single gene, a homeobox transcription factor (GB51409) was differentially expressed across all three tissues.
Evidence that low aggression is a sickness behavior
Are low aggression bees infected with a pathogen?
We detected five bacterial pathogens, four fungal pathogens, deformed wing virus, and acute bee paralysis virus in all three tissues in at least one individual in our study (Table 1). No pathogen was detected in every individual in any tissue, but most pathogens were present in at least one tissue in one individual. No pathogen was significantly more abundant or more likely to be present in low aggression samples (Table S6-S8), suggesting molecular differences as a function of aggression were not caused by acute pathogen infection.
Do low aggression bees show evidence of heightened immune activity?
To evaluate whether the molecular patterns associated with low aggression resemble a diseased state, we compared our differentially expressed gene lists with a recently published meta-analysis that identified genes for which expression changed in response to immune challenge across a variety of tissue types and combinations of tissues, including the whole bee, whole abdomen, fat body, midgut, and brain (Doublet et al. 2017). This meta-analysis identified 57 genes consistently upregulated and 110 genes consistently downregulated in response to immune challenge (across a range of tissue types, see METHODS), whether the source was parasitic mite feeding, viral or fungal infection, or some combination. We performed two enrichment tests per tissue type in our study, evaluating significance in overlap between our differentially expressed gene lists and the up and down-regulated genes from Doublet et al. (2017). We also evaluated directional concordance, with the hypothesis that genes upregulated with infection would be upregulated in low aggression bees if it is a phenotype associated with disease.
In the brain, only one differentially expressed gene overlapped with the Doublet et al. (2017) upregulated gene list, significant overlap due to the relatively small number of differentially expressed genes in this tissue (particularly after list conversion, see METHODS, hypergeometric test, P=0.03). This single gene, GB42523 (an uncharacterized non-coding RNA), was upregulated in low aggression bees, consistent with a hypothesis of a sickness phenotype. Two genes overlapped with the downregulated Doublet et al. list (P=0.01). GB45913 (lethal (2) essential for life, related to heat-shock proteins) was downregulated in low aggression bees, while the second, GB50116 (chymotrypsin inhibitor) was upregulated in low aggression bees.
In the fat body, 13 genes overlapped with the 56 upregulated genes in the Doublet et al. list (Table 2). This overlap was statistically significant (hypergeometric test, P=0.04). Moreover, 10 of the 13 genes were upregulated in low aggression bees, 77% directional concordance with the hypothesis that the fat body molecular signature of low aggression resembles sickness (a significant directional bias, binomial test, P<0.05). 17 genes overlapped with the downregulated Doublet et al. list (out of 110), but this was not statistically significant (P=0.39), nor was the degree of directional concordance (Table 3, 64%, P=0.17). Notably, one gene, hymenoptaecin, was listed on both the up and down-regulated gene lists in Doublet et al. (2017).
In the midgut, 3 genes overlapped with the 56 upregulated Doublet et al. (2017) genes (hypergeometric test, P=0.06). These were GB42523 (uncharacterized), GB48134 (L-lactate dehydrogenase), and GB44112 (melittin); all three were up-regulated in low aggression bees. 7 genes overlapped with the downregulated Doublet et al. (2017) genes (hypergeometric test, P=0.007). These were GB59710 (protein scarlet), GB42053 (NPC intracellular cholesterol transporter 2), GB47279 (cytochrome P450 6k1), GB40976 (HSP90), GB52023 (cytochrome P450 6AQ1), GB49854 (alpha-amylase), GB44549 (glucose oxidase). Five of seven showed concordance with the hypothesis that low aggression resembles sickness (a non-significant result, P=0.23).
Do low aggression bees show differences in predator-responsive genes?
The pre-adult developmental environment could cause low aggression by modulating the baseline expression of genes that are responsive to alarm and predator cues. To test this possibility, we compared our list of genes differentially expressed in the brain as a function of aggression to genes differentially expressed following alarm pheromone exposure (Alaux et al. 2009), which induces a rapid, aggressive anti-predator response. Two genes (GB40074, hormone-like receptor in 38 and GB45913, protein lethal(2)essential for life) overlapped, a non-significant result (P=0.09).
Do pre-adult and adult colony environment effects on aggression share a molecular signature?
High-aggression "Africanized" honey bees are genetically distinct from relatively docile honey bees of European origin. Using a series of experiments that involved housing adult worker bees of Africanized or European origin in either Africanized or European colonies for several weeks, Alaux et al. (2009) found that certain genes in the brain are differentially expressed as a consequence of colony environment, irrespective of individual genotype. This social treatment also affected expression of aggression (Hunt et al. 2003; Alaux et al. 2009). We compared genes differentially expressed as a function of adult colony environment to those differentially expressed as a function of aggression in our study to determine if similar genes are regulated by the adult and pre-adult social environment. Four genes were shared across these lists (GB54316, cardioacceleratory peptide receptor, GB43805, membrane metallo-endopeptidase-like 1, GB41643, blue sensitive opsin, GB54675, uncharacterized), but this degree of overlap was not significant (P=0.19).
Do low aggression bees show a change in rate of adult behavioral maturation?
Adult workers shift tasks as they age, a process called behavioral maturation. This process is influenced by social and environmental cues (Schulz et al. 1999; Huang & Wang 2015), genotype (Giray et al. 2000), and various stressors (Woyciechowski & Moroń 2009; Goblirsch et al. 2013). Older workers performing foraging tasks are typically more aggressive than younger hive bees, and accelerated transition to foraging is associated with higher aggression (Giray et al. 2000). Juvenile hormone regulates this behavioral progression as well as larval development, suggesting these processes could be mechanistically linked. Alternatively, stress tends to accelerate the transition to foraging, which may manifest as accelerated development in low aggression individuals in our study. To assess whether bees show differences in adult behavioral maturation as a function of pre-adult environment, we compared differentially expressed genes in the brain to those differentially expressed between foragers (older adult workers) and nurses (younger adult workers)(Alaux et al. 2009). We found that seven genes (Table 4) overlapped between these lists, a statistically significant result (P=0.01). Five out of seven genes showed directional concordance between low aggression bees and younger nurse bees, suggesting low aggression bees are developmentally delayed, however this relationship was not statistically significant (P=0.23).