Division of labor in honey bees is associated with transcriptional regulatory plasticity in the brain

A strong reciprocal relationship exists between behavioral plasticity and changes in brain gene expression (Cardoso et al., 2015; Zayed and Robinson, 2012), yet the regulatory mechanisms connecting the two are not well understood. Shifts in the expression of specific transcription factors (TFs) can initiate and maintain changes in brain transcriptomic state and lead to corresponding alterations in behavior (Zayed and Robinson, 2012), implicating the involvement of brain transcriptional regulatory networks (TRNs) in the regulation of neural and behavioral state. Although it is well known that TFs influence behavior, it is not yet known whether the functional connections between TF and target genes represented in brain TRNs change quantitatively (i.e. changes in connection strength) or qualitatively (i.e. shifts in the polarity or addition/loss of connections) to give rise to transcriptional regulatory plasticity in association with different behavioral states.

Transcriptional regulatory plasticity (or ‘rewiring’) has been clearly demonstrated at evolutionary time scales, mediated by the addition or removal of TF binding sites in promoter and enhancer regions (Sorrells and Johnson, 2015). In addition, the relationship between TFs and target genes during development can be transiently influenced by epigenetic factors that affect the accessibility of TF binding sites (Araya et al., 2014) or changes in the expression of cofactors and partner TFs that influence the affinity of a TF for specific regulatory sequences (Sorrells and Johnson, 2015). These forms of transcriptional regulatory plasticity have been observed, for example, in the establishment of cell identity leading to cell type-specific patterns of TF–target gene expression (Araya et al., 2014; Lorberbaum et al., 2016).

We hypothesized that context-dependent shifts in nervous system function associated with behavior may also be linked to changes in TF–target relationships over physiological time scales. It is now known that many of the processes that induce regulatory plasticity during development are capable of acting in terminally differentiated cells such as neurons (Sweatt, 2013) at relatively short time scales (Meadows et al., 2016). This suggests that transcriptional regulatory plasticity may also occur in the context of behavior, but it is not yet known whether this is the case.

The nervous system is highly dynamic; the strength of synaptic connections can shift within seconds of neurotransmitter release, and dramatic changes in neural connectivity have been reported in the context of behavior in many organisms, including the honey bee (Fahrbach et al., 1998). Thus, although transcriptional regulatory plasticity has been demonstrated in other contexts, it is not known whether it also can occur on the faster time scales associated with neural and behavioral plasticity. Further, determining whether transcriptional regulatory plasticity in the brain is linked to social behavior is critical for understanding how the nervous system generates adaptive behavioral responses to complex social stimuli.

Social insects such as honey bees are excellent models to study the mechanisms underlying transcriptomic and behavioral plasticity. Colonies of the Western honey bee (Apis mellifera) exhibit a division of labor that is based on the highly stereotyped behavioral maturation of the worker bees, which perform different tasks necessary for their colony's growth and development as they age. During the first approximately 3 weeks of adult life, a worker bee specializes for several days on brood care and other in-hive tasks. Later, it switches to nest defense and/or foraging for nectar and pollen for the remainder of its 5–7 week life. In addition, bees are able to speed up, slow down or reverse behavioral maturation in response to the needs of their colony. The duration of their stable and relatively long-term behavioral states thus results from a complex interplay of genotypic, endocrine, neurochemical, developmental, social and abiotic factors during development and adulthood (Robinson, 1992). The fact that these strong behavioral states are reliably associated with large-scale changes in brain gene expression (Zayed and Robinson, 2012) makes honey bees particularly useful for this study.

A honey bee brain TRN was previously modeled based on extensive brain transcriptomic data that encompassed multiple behaviors and social contexts (Chandrasekaran et al., 2011). Computational inferences from this TRN unexpectedly implicated the possibility of transcriptional regulatory plasticity as a function of behavior. To investigate the concept of transcriptional regulatory plasticity in the context of the brain and social behavior, we performed the following analyses. We first used bioinformatic analyses of the previously published brain transcriptome data (Chandrasekaran et al., 2011) to test whether individual TFs exhibit signatures of regulatory plasticity in the adult bee. We then performed experimental tests by perturbing two TFs identified by the bioinformatic analyses with RNA interference (RNAi) and endocrine treatments. We determined the effects of these perturbations on brain gene expression, behavior and transcriptional regulatory plasticity using field and laboratory assays for behavioral maturation (age at onset of foraging and brood care), and a laboratory assay for hive defense (aggression toward a foreign bee paired with exposure to alarm pheromone). We also explored the idea that pleiotropic TFs can regulate discrete behavioral states by recruiting downstream regulators with more specific functions.

Bioinformatic analyses were conducted on transcriptomic results from our previously published studies (Alaux et al., 2009; Chandrasekaran et al., 2011), which used custom honey bee microarrays (Array Express accession numbers E-TABM-604, E-TABM-605, E-TABM-606 and E-TABM-607). These experiments assessed the interaction between genotype and environmental factors using cross-fostered bees in three different behavioral categories: aggressive bees (soldiers, guards and responders to alarm pheromone), foragers and hive bees (4 day old bees that were not assayed for behavior but were most likely performing brood care and related tasks within the colony; Seeley, 1982). Although the experiments of Alaux et al. (2009) were not designed to assay transcriptomic differences between aggressive bees, foragers and hive bees, they were performed by the same individual, used bees of the same genotypes and colonies of origin, and were performed during the same 3 month period, thereby minimizing the impact of many variables that could interfere with our assessment.

To further ensure that the effect of these variables on brain gene expression was minimized, we processed mean-centered expression values from the microarrays using RUVcorr (Freytag et al., 2015), an implementation of the Remove Unwanted Variation (RUV) algorithm (Gagnon-Bartsch and Speed, 2012) that is designed for use in gene co-expression analyses (File S1, figshare). This algorithm uses internal reference genes that are unlikely to vary across experimental conditions to correct for technical variability in microarray data. We used a set of 1000 probes identified as having the lowest variability across all samples (excluding the genes previously identified as putative targets for the TFs in question). Ridge Regression (Freytag et al., 2015) was then performed to adjust the expression of the remaining genes (Fig. S1A,B), and these values were used in downstream analyses (File S1, figshare).

We modeled changes in putative TF–target relationships (first identified in a honey bee brain TRN; Chandrasekaran et al., 2011) across the behavioral categories using analysis of covariance (ANCOVA). TF expression was treated as a continuous linear predictor of putative target gene expression, and the behavioral category was assessed as a discrete variable. A significant interaction term (false discovery rate, FDR<0.1) between these two variables indicates a change in the regulatory relationship between the TF and putative target genes between categories, i.e. it is ‘context dependent’ (see Fig. 1A for one TF–target gene example). The assumption of a linear relationship between a TF and its target genes is an oversimplification, as genes are regulated by multiple TFs and cofactors. However, linear techniques, together with algorithms designed to ensure that only direct TF–target relationships were inferred, were previously used by Chandrasekaran et al. (2011) to accurately predict TF–target gene relationships for about 25% of the genes expressed in the honey bee brain. This implies that the assumption of linearity also can be used as a heuristic tool to probe transcriptional regulatory plasticity.

We performed these analyses on 20 TFs (Table 1) that were previously predicted to be major regulators of the honey bee brain transcriptome; combined, they were predicted to regulate ca. 75% of the genes in the honey bee brain TRN (Chandrasekaran et al., 2011). Although we corrected for technical variation in the microarray data as detailed above, we cannot fully rule out the possibility that some individual instances of plasticity between the TFs and their putative targets might be false positives. We therefore related TF expression (File S2, figshare) to a subset of microarray probes that we curated to consist of both the TFs’ putative targets and the transcriptomic background (the rest of the microarray probes; 9544 probes in total; Files S2–S4, figshare) from the microarrays used in the original experiments reported in Chandrasekaran et al. (2011). PCA plots based on the putative target and background genes exhibited similar sample distributions, suggesting neither set of genes was biased by technical variation (Fig. S1C,D), and that group-wide comparisons between them should be a valid method to detect whether a TF demonstrated a signature of transcriptional regulatory plasticity.

Signatures of transcriptional regulatory plasticity were inferred for each pair of categories (aggressive bees, foragers and hive bees) by comparing the proportion of putative target genes (Files S3 and S4, figshare) identified in Chandrasekaran et al. (2011) with significant (FDR<0.10) differences in expression for the pairwise contrast with the proportion of the remaining probes (i.e. the transcriptomic background) for the same contrast (one-sided Fisher’s exact test, FDR<0.05). Only if at least 15% of a TF's putative target genes exhibited a significant contrast between the relevant pair of behavioral categories was it considered for this analysis. Two TFs (Broad and NF-κB) were represented by multiple probes on the array. To account for this, the analysis was run independently for each individual TF probe, and significant effects were compiled across the putative target gene and background sets for Fisher’s exact test. We also determined the effect size for each comparison by comparing the odds ratio (OR), calculated by dividing the proportion of putative target genes with significant TF×category interaction effects by the proportion of the transcriptomic background that had significant interaction effects.

All experiments took place at the University of Illinois Bee Research Facility, Urbana, IL, USA, between the months of June and October in 2013 and 2014. To minimize the effects of genetic variation on behavior and molecular analyses, all experiments used adult worker bees, Apis mellifera Linnaeus 1758, from colonies headed by queens each instrumentally inseminated with semen from a single drone; because of haplodiploidy, this results in an average coefficient of relatedness of 0.75. Honeycomb frames containing pupae were removed from colonies 1–4 days prior to the beginning of each experiment and maintained in a dark incubator at 34°C and 50% relative humidity. Adult workers 0–24 h old were obtained from these frames, placed into Plexiglas cages in groups of 50, and fed 50% (w/v) sucrose solution and pollen paste (45% honey/45% pollen/10% water) ad libitum. The bees were kept in these cages in the incubator for 4 days prior to RNAi treatment, and an additional 2 days post-treatment (except for the bees used in the foraging experiments, see below).

We focused on the transcription factor Broad to experimentally test for a direct connection between network plasticity and behavior, for four reasons. First, Broad was one of only two TFs to exhibit a signature of plasticity between more than one pair of behavioral categories in the bioinformatic analysis described above (see Results). Second, Broad is known to be a pleiotropic integrator of diverse endocrine cascades in multiple insect species (Bonneton and Laudet, 2012). It mediates insulin, juvenile hormone (JH) and vitellogenin (Bonneton and Laudet, 2012; Hamilton et al., 2017) signaling, all of which are known to play a causal role in honey bee behavioral maturation (Hamilton et al., 2017; Zayed and Robinson, 2012). Third, Broad also has consistently been predicted to be an important regulator of honey bee behavior in studies employing diverse bioinformatic methods, and its predicted targets are associated with neuronal function and activity (Ament et al., 2012a; Chandrasekaran et al., 2011). Fourth, an increase in the number of Broad cis-regulatory elements across the genome has previously been linked to the evolution of social behavior in multiple species of bees (Kapheim et al., 2015), implying that its regulatory relationships have been ‘rewired’ over evolutionary time scales.

We also targeted a TF related to Broad that did not exhibit a strong signature of transcriptional regulatory plasticity (see Results), to complement and contrast with Broad: Fushi tarazu transcription factor 1 (Ftz-F1). Ftz-F1 is a downstream target of Broad in highly conserved endocrine cascades that are critical to insect development (Bonneton and Laudet, 2012) and, like Broad, it is predicted to be a key regulator of several behaviors in honey bees (Ament et al., 2012a; Chandrasekaran et al., 2011; Faragalla et al., 2018). Additionally, like Broad, Ftz-F1's putative targets are enriched for gene ontology terms related to neuronal plasticity, cognition and learning and memory (Chandrasekaran et al., 2011).

Dicer-substrate (Integrated DNA Technologies, Coralville, IA, USA) RNAi constructs (D-siRNA) were designed against exons shared between all predicted isoforms of broad complex or fushi tarazu transcription factor 1 (ftz-f1). A cocktail of three such constructs (Table S1) in equivalent concentrations was used to target each gene. Perfluorocarbon nanoparticles were fabricated as previously described (Kaneda et al., 2010). To create nanoparticle–D-siRNA complexes, 1 μmol l⁻¹ of the D-siRNA cocktail or an exogenous control (targeted to green fluorescent protein, gfp) and 2 nmol l⁻¹ nanoparticles were combined in 1× phosphate buffered saline (Sigma-Aldrich, St Louis, MO, USA). This solution was incubated at room temperature for 30 min with light agitation before use. These nanoparticles have been used previously to successfully administer siRNA to honey bees (Li-Byarlay et al., 2013).

In the 2013 experiments, 20% high molecular weight dextran (Sigma-Aldrich) was added to the nanoparticle complexes post-incubation to increase the viscosity of the solution and prevent it from diffusing away from the head capsule during injection. However, dextran use resulted in a substantially higher mortality, so this practice was discontinued in 2014.

We manipulated Broad and Ftz-F1 expression by applying broad and ftz-f1 RNAi to the brain via direct injection, a method that has previously been shown to reduce gene expression in the honey bee brain (Farooqui et al., 2003; Rein et al., 2013). Bees were cold-anesthetized and mounted on a Plasticine (Flair Leisure Products PLC, Cheam, UK) pedestal; dental floss was used to stabilize the head. Pedestals were placed in a tray of dry ice to ensure that the bees did not wake prior to the end of the injection. Prior to injection, two paint marks (Testors, Rockford, IL, USA) were applied to the individual's thorax and a preliminary incision was made to the posterior part of the median ocellus using a 28-gauge needle. A 34-gauge Nanofil needle (World Precision Instruments, Sarasota, FL, USA) was inserted 500 μm into the incision, and 500 nl of nanoparticle complexes was administered using a UMP3 Microinjector run by a Micro4 controller (World Precision Instruments) at a rate of 5 nl s⁻¹. Bees that exhibited backflow of the injection mixture were discarded. After injection, bees were kept anesthetized for ca. 10 min before being returned to their cages to prevent the injected solution from circulating away from the brain. Mortality as a result of the injections was ca. 15–20%, with no statistically significant differences as a function of treatment type (Kruskal–Wallis test; Fig. S2). The efficacy of these injections, measured with qPCR on broad and ftz-f1 expression resulted in significant knockdowns of 17–37% in 80% of the trials (Table 2; Fig. S3).

Broad and Ftz-F1 mediate the transcription-dependent effects of JH on cell physiology during development in other insect species (Bonneton and Laudet, 2012) and (at least in the case of Ftz-F1) in pharate adult bees (Mello et al., 2018). JH plays a strong role in regulating honey bee behavioral maturation (Robinson, 1992). We co-administered RNAi and the JH analog methoprene (JHA) and observed effects on brood care behavior. Methoprene was used instead of JH-III because of its longer half-life and the fact that it is known to induce many of the same transcriptional changes as behavioral maturation (Whitfield et al., 2006). In 2014, bees in half of the cages were administered methoprene using a method adapted from established protocols (Ament et al., 2012b). Methoprene was mixed with the pollen paste (20 mg g⁻¹) until the day prior to RNAi treatment, at which point it was replaced with standard pollen paste. After RNAi treatment, methoprene-treated bees and control diet bees were combined in the same cage and co-housed for the remainder of the experiment. Previous research has shown that similar doses of methoprene cause precocious foraging and forager-like gene expression patterns involving hundreds of genes in both the brain (Whitfield et al., 2006) and fat bodies (Ament et al., 2012b).

Single-cohort colonies initially composed of ca. 2000, 1 day old worker bees were established according to previous protocols (Tenczar et al., 2014). Each colony was given two honeycomb frames of honey, one frame of pollen, one empty frame and a mated queen unrelated to the worker bees. Each colony was formed 2 days prior to the eclosion of the treated bees. Instead of being housed in an incubator, Plexiglas cages containing these bees (provisioned as noted previously and given a piece of comb from the parent colony) were housed within the colony prior to RNAi treatment. Immediately after RNAi treatment, treated and control bees were each tagged with a pair of RFID (p-chip) transponders (Pharmaseq, Princeton, NJ, USA) (Tenczar et al., 2014). The bees were then returned to their cages and placed within the colony until their release 12 h later. Together with the presence of common pieces of comb, the co-housing procedure ensured that the other bees in the colony were not aggressive toward the treated bees. Up to three cohorts of bees were added to each colony on consecutive days. The flight activity of the treated and control bees was continuously monitored by a pair of barcode readers (Pharmaseq) mounted at the hive entrance that captured all outgoing and incoming trips (Tenczar et al., 2014). Data were collected continuously for 10 days post-RNAi treatment, and age at onset of foraging information was extracted as described previously (Tenczar et al., 2014) and summarized below.

Two days after RNAi treatment, bees were transferred to Petri dishes (100×20 mm) provisioned with honey, water, a pollen ball and a piece of wax comb from their original cage. In 2013, groups of six individuals (three control RNAi and three treatment RNAi) were used; in 2014, groups of eight individuals were used, consisting of two bees from each JHA×RNAi treatment combination, yielding a 2×2 factorial. Each group member was marked uniquely to enable behavioral analyses at the individual level. The bees were then placed in a climate-controlled testing room (28°C) with normal fluorescent lighting and allowed to acclimate for at least 60 min prior to testing.

Aggression assays were adapted from a previously developed intruder assay (Rittschof et al., 2014). An intruder collected from the entrance of an unrelated colony immediately prior to the assay was introduced to a dish of bees. The bees were observed for 5 min, after which the intruder was removed from the dish, regardless of whether it was still alive. The bees were then exposed to 5 μl of isopentyl acetate (Sigma-Aldrich), the active ingredient in honey bee alarm pheromone, diluted 1:10 with mineral oil as described previously (Alaux et al., 2009). This additional treatment was added to simulate a high-level threat to the group, even if none of the residents responded to the intruder in an aggressive manner. Previous studies have shown that this intruder assay produces behavior that is consistent with aggressive defense behavior in colonies in the field (Rittschof et al., 2014).

In 2013, treatment and control bees were collected 1 h after the presentation of the stimulus by flash-freezing them in liquid nitrogen. Unexposed control groups served as a baseline for gene expression assays. These groups were acclimated to the same chamber and collected during the same time frame, but were not administered an intruder, isopentyl acetate or control stimulus. In 2014, bees were instead collected either immediately after the end of the assay (before the assay could elicit changes in transcription) or 60–120 min later. Those groups collected immediately after the end of the assay served as the baseline for gene expression studies.

Brood care assays were performed as described previously (Shpigler and Robinson, 2015). A waxen naturally built cell containing a 3 or 4 day old queen larva was introduced into a Petri dish of bees, and the interactions between the residents and the queen were scan-sampled at 10 s intervals. The queen cell was removed after 10 min of observation. Control dishes were each given an empty queen cell to determine whether the RNAi treatment caused any differences in general sensitivity to stimuli. Previous studies have shown that this assay produces behavior that is consistent with brood care behavior in colonies in the field (Shpigler and Robinson, 2015). Collections followed the same protocol as the aggression assays performed in 2014, described above.

Sample preparation and quantitative PCR (qPCR) were performed using established protocols (Ament et al., 2011). Whole brains were dissected into RNA-later ICE (Thermo Fisher Scientific, Waltham, MA, USA) at −80°C, homogenized and extracted using RNeasy spin columns (Qiagen, Valencia, CA, USA) with DNase treatment (Invitrogen, Carlsbad, CA, USA). RNA was quantified using a Nanodrop spectrometer (Thermo Fisher Scientific) and Qubit fluorimeter (Thermo Fisher Scientific), and cDNA syntheses were performed using Arrayscript (Thermo Fisher Scientific) reverse transcriptase. An exogenous RNA spike-in (Root Cap Protein 1 from Arabidopsis thaliana) was used to assess cDNA synthesis efficiency. qPCR was performed using SYBR Green dye (Thermo Fisher Scientific) on an ABI Prism 7900ht (for samples analyzed in 2013) or ABI QuantStudio 6 (for samples analyzed in 2014). Sample–probe combinations with an inter-triplicate coefficient of variation (CV) higher than 30% were discarded. For qPCR probes, see Table S1.

The age at onset of foraging was determined with an algorithm that relates the proportion of RFID p-chip reads during peak foraging hours to the proportion during the time when pre-foraging orientation flights are likely to occur (Tenczar et al., 2014). Foragers were classified as bees with more than eight reads for at least 2 consecutive days, with greater than half of their reads collected during peak foraging hours (12:00 h–15:00 h central standard time) (File S5, figshare). A previous study using similar, but less stringent, thresholds found an excellent correspondence between this metric and human observations (Tenczar et al., 2014). Up to three cohorts of bees were treated with RNAi and/or JHA on consecutive days and added to each background colony. The age at onset of foraging data were analyzed with a Cox Proportional Hazards (coxph in R; File S1, figshare) model (Ben-Shahar et al., 2002), stratified by the colony the bees were added to and the cohort they belonged to within this colony. Bees that died before the end of the experiment and did not become foragers were removed from the analysis, as the time of their death could not be accurately determined.

Behavioral data from the aggression assay were scored using an index adapted from Rittschof et al. (2014). Bees were assayed using an approach that involves scan sampling for 5 min, with an observer blind to treatment noting all interactions between residents and the intruder during consecutive 10 s scans (Table 3). In addition, these values were summed for each 10 s window that an individual resident interacted with the intruder, and the total score for each bee was normalized by the duration of the potential interaction time (i.e. if a resident or intruder died prior to the end of the assay, additional bins would not be factored into the normalized score; File S5, figshare). The average score for each treatment group of bees was used as the experimental unit, with individual bees regarded as subsamples. Because of the invasive nature of the RNAi treatments and the potential for treatment-related injuries to interfere with normal behavior, individuals that did not respond to the intruder (30% across all assays, with no significant difference between treatments) were not scored and were discarded from the study. Brood care behavior was measured and scored in the same manner as aggression, using a behavioral index adapted from Shpigler and Robinson (2015) (Table 3; File S5, figshare).

The distribution of aggression and brood care scores across samples deviated significantly from normality, so a Box–Cox transformation was applied as a correction. Post-transformation, all data from 2013 met requirements of normality and homoscedasticity. Behavioral data from 2014 met all assumptions post-transformation, except for skewness and, in some instances, kurtosis (likely due to the lower number of pseudoreplicates in each sample); however, mixed models have been shown to be robust to violations of skewness and kurtosis when total sample sizes larger than 45 are used (Arnau et al., 2013). Therefore, we analyzed the transformed data using a linear mixed effects model (the lme4 package in R; Bates et al., 2015) with Type III Sums of Squares and the Kenward–Roger approximation of degrees of freedom (reported to the nearest degree of freedom; File S1, figshare). To account for group-wise differences in variables such as the stimulus presented (which cannot be easily quantified or controlled for across groups), we included a random variable that accounted for the dish that each set of bees belonged to (nested within the date of testing and the bees’ colony of origin). We also initially blocked for the experimental observer and included the time of testing as a linear covariate. However, the influence of observer and the time of testing were not found to be significant, so these variables were removed from the final model.

The effects of the RNAi (broad and ftz-f1 RNAi) were analyzed independently, as the injections were performed on distinct sets of bees. In 2014, JHA treatment was added as a response variable of interest, and the interaction between JHA and broad/ftz-f1 RNAi was modeled as well.

Expression values for each gene were normalized to the geometric mean expression of four reference genes (S8, RP49, GAPDH and eIF1-a). Each reference gene was assessed individually to ensure that the CV was less than 40% in each experiment/colony combination, and there were no significant differences in reference gene expression between experimental groups (File S6, figshare).

Six predicted target genes of Broad and six predicted target genes of Ftz-F1 were selected based on their strong relationship (a regression coefficient in the 10th percentile) with their parent TF in previous computational analyses (Chandrasekaran et al., 2011). The influence of the TFs on these putative target genes was determined in two ways: (1) using a multifactorial ANOVA and looking at the effect of RNAi on target gene expression or (2) modeling TF expression levels as a continuous variable and using an ANCOVA to determine the relationship between TF and target gene expression. In both cases, a fixed effects model was used with Type III Sums of Squares (File S1, figshare). Sample colony of origin was included as a blocking factor; time and date of testing were analyzed but did not have a significant effect on gene expression and were dropped from the final models.

Transcriptional regulatory plasticity between Broad, Ftz-F1 and their target genes was determined by conducting several different behavioral and endocrine manipulations and observing their effect on the TF–target gene relationships. We tested for context-dependent effects of behavioral state (how the bee responded to a social stimulus, i.e. an intruder or a queen cell), social environment (i.e. the mere presence of an intruder or queen cell, independent of behavioral response) and endocrine state (whether the bees were treated with JHA or not). This was accomplished by examining the two-way interactions of each of these response variables with TF expression (via ANCOVA) as in the bioinformatic analyses detailed above.

During the 2013 field season, the baseline for social context was provided by matched control groups that were not exposed to an intruder or alarm pheromone but were otherwise treated in the same manner. The effect of behavioral state on gene expression was therefore nested within stimulus exposure (as behavioral state could not be ascertained for the unexposed group). In 2014, by contrast, the data were organized in a time series that was analyzed as a split-plot in time with three levels (see information on collections above for differences between 2013 and 2014), allowing for a direct comparison of individuals before and after gene expression changes were elicited by the social stimulus. Therefore, both behavioral state and the response to stimulus could be resolved, allowing each to be treated as a main effect, with their interaction term delineating the difference between the way individuals were influenced by the social stimulus. Additionally, the effect of methoprene on gene expression and the TF’s regulatory relationships was also assessed. Working with a small set of genes made it more feasible to use gene expression analysis to explore a variety of factors that might affect network plasticity.

In the behavioral maturation experiment, the efficacy of RNAi treatments was confirmed by sampling treated bees on the day they were released into the colony. The influence of RNAi was analyzed as above but blocking only for colony of origin.

To ensure that broad knockdowns did not induce global changes in transcription, we also assessed the expression of a panel of six predicted target genes of a third TF (CG17912) that were not predicted targets of Broad and Ftz-F1. The fact that only one of these genes was differentially expressed in response to RNAi indicates that the effect of these treatments was largely specific to the putative targets (Fig. S4; File S5, figshare).

Signatures of transcriptional regulatory plasticity were found for 5 out of 20 TFs: Broad, Yl-1, Trithorax-like (Trl), Creb and MTA1-like (Fig. 1B, Table 1). Broad and Yl-1 each exhibited evidence of plasticity across two pairs of behavioral categories. Broad's putative regulatory relationships varied between aggressive bees and foragers (OR=1.57, FDR=0.014) as well as between aggressive bees and hive bees (OR=1.89, FDR=2.94e⁻⁴), and Yl-1's putative regulatory relationships varied between aggressive bees and hive bees (OR=2.02, FDR=0.037) as well as between foragers and hive bees (OR=2.14, FDR=0.046). Trl's putative regulatory relationships varied between aggressive bees and hive bees (OR=1.98, FDR=2.38e⁻³), and Creb's and MTA1-like's putative regulatory relationships varied between aggressive bees and foragers (OR=1.52, FDR=0.030 and OR=1.63, FDR=0.048, respectively). Broad and Yl-1 thus demonstrated the strongest evidence of plasticity, as the only TFs to exhibit a signature of plasticity between more than one pair of behavioral categories.

In addition to these results, ANCOVA revealed that the putative targets of 14 out of the 20 TFs we studied were overrepresented relative to the transcriptomic background when testing for a main effect of TF expression (Table 1). These results are largely concordant with the predictions of the bee brain TRN (Chandrasekaran et al., 2011) despite using different statistical methodology and only a subset of the samples that the original TRN was trained on. The fact that we observed extensive significant main effects for the TFs on their putative target genes also serves as a contrast to the signatures (interaction effects) of transcriptional regulatory plasticity detected for Broad, Yl-1, Trl, Creb and MTA1-like.

As predicted, we found that broad RNAi influenced both behavioral maturation (age at onset of foraging and brood care) and aggression, in automated field behavior analyses using RFID tags (Tenczar et al., 2014) and laboratory behavioral assays. broad RNAi caused a significantly earlier foraging onset age (Cox Proportional Hazards: z=4.67, P<1e⁻⁵; Fig. 2A). However, bees in different treatment groups showed no foraging differences 3 days after RNAi treatment (Cox Proportional Hazards, P=0.34), suggesting that the effect of RNAi on the age at onset of foraging is delayed. broad RNAi also caused a significant decrease in the intensity of aggression toward intruder bees (linear mixed effects model: F_1,45=19.95, P<1e⁻⁴ for data collected in 2013 and F_1,245=10.27, P<1e⁻² for data collected in 2014; Fig. 2B, Table 2) and a significant increase in the intensity of brood care (linear mixed effects model: F_1,213=26.98, P<1e⁻⁶; Fig. 2D, Table 2).

As with broad RNAi, ftz-f1 RNAi also significantly decreased the age at onset of foraging (Cox Proportional Hazards: z=3.77, P<1e⁻³; Fig. 2A) and increased the intensity of brood care (linear mixed effects model: F_1,173=20.10, P<1e⁻⁴; Fig. 2E, Table 2). In addition, as predicted, ftz-f1 RNAi did not alter overall levels of aggression (linear mixed effects model: F_1,42=0.67, P>0.4; Fig. 2C, Table 2).

To explore the hypothesis that JH alters behavioral maturation at least in part through Broad- and Ftz-F1-dependent pathways, we co-administered broad or ftz-f1 RNAi and JHA and observed the effects on brood care behavior. Consistent with this hypothesis, JHA partially rescued the effect of RNAi on brood care (Fig. 2D,E, Table 2; linear mixed effects model, broad RNAi: F_1,220=13.50, P<0.01; ftz-f1 RNAi: F_1,174=6.31, P<0.05).

RNAi treatment altered the expression of 11 out of 12 predicted target genes, confirming their status as targets and providing additional support for the conclusions of our bioinformatics and behavioral analyses. Importantly, broad knockdown also reduced the expression of ftz-f1 and many of the ftz-f1 target genes (Fig. 3A) but not vice versa (Fig. 3B). This result indicates that the previously reported hierarchical relationship of these two TFs in Drosophila melanogaster (Bonneton and Laudet, 2012) also exists in the adult honey bee brain. It also suggests that pleiotropic TFs involved in behavior, such as Broad, induce state-specific transcriptional programs via the recruitment of more specialized downstream regulators (e.g. Ftz-F1).

Using ANCOVA of brain gene expression, and consistent with the results of the bioinformatic analyses presented above, we found that several of the TF–target gene relationships varied as a function of behavioral state (Fig. 4A), providing empirical evidence for transcriptional regulatory plasticity. Additionally, we detected transcriptional regulatory plasticity as a function of social context independent of behavioral response (the presence of an intruder and alarm pheromone that provokes aggression or a queen larva that provokes brood care; Fig. 4B) and as a function of endocrine state (treatment with JHA independent of behavioral effect; Fig. 4C).

Our findings extend the paradigm of transcriptional regulatory plasticity, which is well established in evolutionary biology (Sorrells and Johnson, 2015) and developmental biology (Araya et al., 2014; Lorberbaum et al., 2016; Sorrells and Johnson, 2015), to neurobiology and to the shorter time scales associated with neural and behavioral plasticity. We demonstrated that the precise relationships between a TF and its target genes depend on the behavioral and endocrine state of each individual bee as well as its social context (Figs 1 and 4).

Through bioinformatic analyses, we showed the existence of brain transcriptional regulatory plasticity by examining changes in putative TF–target gene relationships associated with behavior. Although we lack evidence that these TFs bind directly to the cis-regulatory modules of their predicted targets, TF–target gene relationships were inferred by algorithms designed to find direct TF–target gene connections (Chandrasekaran et al., 2011). Moreover, the results of the TRN that produced these inferences (Chandrasekaran et al., 2011) are generally consistent with cis-regulatory analyses of genes differentially expressed across behavioral contexts (Ament et al., 2012a; Khamis et al., 2015). Together with the results of our RNAi experiment (Fig. 3), this suggests that a significant number of the predicted TF–target gene relationships represent direct interactions.

Signatures of transcriptional regulatory plasticity were found for only five out of 20 TFs analyzed: Broad, Creb, MTA1-like, Trl and Yl-1, suggesting that the phenomenon of transcriptional regulatory plasticity may be preferentially associated with specific subsets of TFs. Orthologs of these five TFs in Drosophila are either highly pleiotropic integrators of cell signaling pathways (Broad and Creb) or epigenetic regulators of gene expression (Yl-1, Trl and MTA1-like). Broad is known to be a critical signal integrator in diverse endocrine cascades (Bonneton and Laudet, 2012), and was one of only four TFs predicted to regulate all three different behavioral categories (aggression, maturation and foraging) in the bee brain TRN (Chandrasekaran et al., 2011). Creb is essential for the transduction of kinase signaling cascades related to neuronal activation, plasticity and learning and memory, and has been implicated in gating behavioral responses to sensory stimuli (Gehring et al., 2016) and as a critical regulator of foraging behavior in honey bees (Khamis et al., 2015). Yl-1 is a histone acyltransferase (Kusch et al., 2004), while MTA1-like is a histone deacetylase (Sen et al., 2014; Swenson et al., 2016) and both are components of nucleosome remodeling complexes (Sen et al., 2014; Swenson et al., 2016; Weber et al., 2014). Trl (the insect ortholog of GAGA factor) is a chromatin remodeler that inhibits Polycomb Group activity (Ringrose and Paro, 2007). Epigenetic factors (including DNA methylation and histone modifications) have been implicated in behavioral plasticity in honey bees and other social species (Herb et al., 2018, 2012; Simola et al., 2016), and are also capable of inducing transcriptional regulatory plasticity (Araya et al., 2014). We therefore speculate that transcriptional regulatory plasticity may rely especially on pleiotropic signal integrators or epigenetic regulators, which can initiate and maintain contextually specific transcriptional programs that depend on the internal state of the individual organism.

RNAi treatment confirmed that Broad is a pleiotropic regulator of honey bee behavior. While broad RNAi influenced both behavioral maturation (i.e. foraging and brood care) and aggression, ftz-f1 RNAi influenced only behavioral maturation (Fig. 2). It may seem puzzling that both brood care and the likelihood of foraging increased as a result of RNAi, because they are generally seen as mutually exclusive behavioral states. Brood care is generally performed when the bee is young, and foraging is performed when older. However, there is evidence to suggest that JH affects behavioral maturation in a dose-dependent manner. There is a dose-dependent effect of the JHA methoprene on the age at onset of foraging (Robinson, 1987). In addition, low doses of JH treatment increase the development of the brood food-producing hypopharyngeal glands, while high doses prematurely decrease their development (Rutz et al., 1976). Based on these results, we suggest that broad and ftz-f1 RNAi increased the likelihood of brood care at young ages and increased the likelihood of foraging at older ages. This is consistent with the finding that bees in different treatment groups showed no foraging differences 3 days after RNAi treatment, which was the age at which brood care was assayed. By contrast, co-administration of RNAi with JHA advanced maturation past the point where brood care is more likely to occur, causing a reversion to baseline levels of this behavior.

By examining changes in the TF–target gene relationship as a function of aggression and brood care, we found empirical evidence for transcriptional regulatory plasticity in the context of behavioral state, endocrine state and social environment for both Broad and Ftz-F1. Because these results were derived from the knockdown of individual TFs, they suggest that the observed changes in behavior and brain gene expression are due to changes in the strength of specific TF–target gene relationships, rather than more indirect effects. The fact that transcriptional regulatory plasticity was observed for these TFs even though only one of them was identified in our bioinformatic screen reinforces the idea that quantitative (rather than qualitative) differences in the prevalence of regulatory plasticity exist, potentially related to a TF's properties.

broad RNAi reduced the expression of Ftz-F1 and its target genes in the brain but not vice versa (Fig. 3), confirming the previously shown hierarchical relationship between Broad and Ftz-F1, with Broad upstream of Ftz-F1 (Bonneton and Laudet, 2012). As the two TFs modulated brood care and the age at onset of foraging behavior in a similar fashion, it is likely that Broad's influence on these behaviors also involves Ftz-F1. Broad's effect on aggression, however, is probably independent of Ftz-F1, as ftz-f1 RNAi did not affect aggressive behavior (though Ftz-F1 may play some role in aggression in honey bees, as it has been found to be differentially expressed in previous aggression-related brain transcriptomic studies: Alaux et al., 2009; Rittschof et al., 2014). These results lead us to speculate that one possible route for a pleiotropic behavioral regulator like Broad to influence multiple behavioral states is via discrete regulatory subnetworks with branching context-specific transcriptional cascades. This hypothesis is consistent with the importance of network hierarchy in determining the contribution of a TF to establishing phenotypic plasticity (Bhardwaj et al., 2010).

The possibility of off-target and deleterious effects always exists when using RNAi, so we note that the broad and ftz-f1 RNAi results reported here can be seen as ‘gain of function’ changes. Causing an early onset of foraging (the most physically and cognitively demanding task in the honey bee repertoire; Robinson, 1992) and an increase in the intensity of brood care (another complex behavior) is evidence that the treated bees were capable of functioning normally. Moreover, it suggests that the RNAi treatment had specific effects on behavioral state. Early onset of foraging in honey bees also occurs in response to various stressors such as pathogens (Goblirsch et al., 2013), starvation (Schulz et al., 1998) and demographic disruption (Huang and Robinson, 1996). However, we used conservative thresholds for the identification of foragers (requiring multiple foraging trips over more than one day; see Materials and Methods), making it unlikely that these bees were substantially impaired by the treatment.

Future studies should focus on modeling plasticity across the entire transcriptome to generate TRNs that capture shifts in regulatory architecture associated with behavior. In addition, an integration of the contributions of site-specific epigenetic modifications (Hamilton et al., 2018) and motif-binding properties of effectors (Khakharn et al., 2018) will be necessary to fully elucidate the causal connections between brain transcriptional regulatory plasticity and behavior.

The findings presented here are derived from whole-brain transcriptomes, so they represent an aggregate of the individual states of each neuronal and glial cell in the brain. Brains are highly compartmentalized and consist of numerous specialized regions and neuronal subtypes. Only a subset of these neurons is likely to be activated in a given social context, and recently it was shown that even neurons of the same subtype and lineage can exhibit strikingly different transcriptomic profiles (Poulin et al., 2016). However, transcriptome-wide differences in brain expression associated with naturally occurring behavior have been reported for a variety of species (Hughes et al., 2012; Oliveira et al., 2016) ever since they were first discovered in honey bees (Whitfield et al., 2003). Given that many behaviors are known to be orchestrated in specific regions of invertebrate and vertebrate brains, why should there be such robust patterns of behaviorally related gene expression at the whole-brain level in honey bees and other organisms? A technical explanation is that the whole-brain transcriptomic profile largely reflects the profiles of the larger brain regions; in honey bees, this would include the mushroom bodies and optic lobes, which together account for ∼2/3 of the neurons in the adult bee brain. A biological explanation is that the whole-brain transcriptomic profile arises because there are similarities in gene expression in different brain regions. For instance, numerous neuromodulators and hormones are known to directly contribute to long-term changes in honey bee behavior (Hamilton et al., 2017), and many are known to have receptors in multiple regions of the adult insect brain (Baumann et al., 2017; Perry and Barron, 2013). This makes it likely that at least some neuromodulators and hormones are capable of coordinating gene expression across brain regions. Determining the role that either or both of these explanations play in contributing to the dynamics of brain transcriptional regulatory plasticity, and neural and behavioral plasticity in general, will require integrating information at the cellular level. Recent advances in single-cell RNA sequencing make this feasible for the first time, and future studies should use these exciting new developments to examine transcriptional regulatory plasticity at the level of individual neurons.

In summary, transcriptional regulatory plasticity in the brain occurs in the context of social behavior. The transcriptomic architecture underlying complex behavior thus appears to be more dynamic than previously appreciated. This makes it important to consider brain transcriptional regulatory plasticity along with well-known neurochemical, neurophysiological and neuroanatomical mechanisms to achieve a more complete understanding of neural and behavioral plasticity.

We thank Daniel C. Nye for assistance with honey bee rearing and colony maintenance in the field, Emma E. Murdoch for assistance with data collection, Michael J. Scott for synthesizing the perfluorocarbon nanoparticles, and Thomas C. Newman and Amy Cash Ahmed for their help and advice with molecular biology. We also thank Claudia C. Lutz, Lisa C. Stubbs, Rhanor Gillette, Justin Rhodes and the members of the Robinson lab for reviewing the manuscript and offering critical feedback. The methoprene used in this study was generously donated by Doug Vangundy of Central Life Sciences (Schamburg, IL, USA).