Skip to main content
  • Research article
  • Open access
  • Published:

Defense against territorial intrusion is associated with DNA methylation changes in the honey bee brain

Abstract

Background

Aggression is influenced by individual variation in temperament as well as behavioral plasticity in response to adversity. DNA methylation is stably maintained over time, but also reversible in response to specific environmental conditions, and may thus be a neuromolecular regulator of both of these processes. A previous study reported DNA methylation differences between aggressive Africanized and gentle European honey bees. We investigated whether threat-induced aggression altered DNA methylation profiles in the honey bee brain in response to a behavioral stimulus (aggression-provoking intruder bee or inert control). We sampled five minutes and two hours after stimulus exposure to examine the effect of time on epigenetic profiles of aggression.

Results

There were DNA methylation differences between aggressive and control bees for individual cytosine-guanine dinucleotides (CpGs) across the genome. Eighteen individual CpG sites showed significant difference between aggressive and control bees 120 min post stimulus. For clusters of CpGs, we report four genomic regions differentially methylated between aggressive and control bees at the 5-min time point, and 50 regions differentially methylated at the120-minute time point following intruder exposure. Differential methylation occurred at genes involved in neural plasticity, chromatin remodeling and hormone signaling. Additionally, there was a significant overlap of differential methylation with previously published epigenetic differences that distinguish aggressive Africanized and gentle European honey bees, suggesting an evolutionarily conserved use of brain DNA methylation in the regulation of aggression. Lastly, we identified individually statistically suggestive CpGs that as a group were significantly associated with differentially expressed genes underlying aggressive behavior and also co-localize with binding sites of transcription factors involved in neuroplasticity or neurodevelopment.

Conclusions

There were DNA methylation differences in the brain associated with response to an intruder. These differences increased in number a few hours after the initial exposure and overlap with previously reported aggression-associated genes and neurobiologically relevant transcription factor binding sites. Many DNA methylation differences that occurred in association with the expression of aggression in real time also exist between Africanized bees and European bees, suggesting an evolutionarily conserved role for epigenetic regulation in aggressive behavior.

Background

Experiences produce molecular responses in the brain that exert long-term influences on behavior, affecting the way individuals or groups respond to future circumstances [1]. For example, acute stress, including exposure to social stressors, can lead to a transient state of increased aggression [2]. This reaction to stress exists in many animal model systems, including honey bees. When faced with a territorial intruder, bees protect their colony with defensive behavior such as biting and stinging. Once disturbed, they develop a state of aggressive vigilance against other potential threats and exhibit heightened responses to future intrusions [3, 4].

Multiple observations motivate investigation into the epigenetics of this response. First, there are extensive differences in gene expression in the honey bee brain that are related to aggression [4,5,6,7,8,9,10], including between the highly aggressive Africanized honey bee (AHB) subspecies Apis mellifera scutellata and the less aggressive European honey bee subspecies (EHB) [5]. In addition, whole-genome profiling of 5-methylcytosine (5mC) has revealed differences in DNA methylation between AHB and EHB [9]. Although these results suggest the possibility that aggression is epigenetically regulated in honey bees, AHB and EHB differ in other traits besides aggression. The question of whether brain epigenetic profiles are associated specifically with aggressive behavior in bees therefore remains open.

There is a well-established relationship between behavioral plasticity and dynamic DNA methylation in several model systems. In mammals, modulation of epigenetic marks in the brain has been linked with exposure to social stressors, including separation of offspring from their mothers [11, 12] and chronic social defeat [13, 14]. Recent evidence indicates that DNA methylation underlies the maintenance of long-term memory [15,16,17,18], allowing experiences to stably alter behavior. The de novo addition and active removal of DNA methylation facilitate transitions between epigenetic states in the brain [19,20,21,22], suggesting a compelling molecular mechanism underlying behavioral plasticity.

Knowledge of genetic and social effectors of aggression in honey bees, coupled with a fully functional epigenetic toolkit [23], make honey bees an excellent organism to study the molecular basis of aggression. DNA methylation changes in the honey bee brain are associated with behavioral plasticity, most notably in switching between nursing and foraging roles in the hive division of labor [24]. Additionally, differences in brain DNA methylation have been reported between queen and worker bees during development and in adulthood when they display distinct behaviors [25, 26]. Early evidence for DNA methylation’s impact on queen / worker development was illustrated by DNMT knockdown during the larval stage that produced queens despite a feeding regime that typically produces workers [27].

Division of labor is a hallmark of eusocial species, and mounting evidence shows that DNA methylation plays a role in this process, possibly by affecting vitellogenin levels and altering lifespan [28]. In the brain, pharmacological perturbation of DNA methylation affects learning ability, and long-term memory formation appears to require upregulation of DNA methylation modifying enzymes DNMT and Tet [29]. This evidence, in addition to DNA methylation’s predicted role in regulating gene expression and alternative splicing events, makes DNA methylation an attractive candidate for regulating aggressive behavior in honeybees [30].

We used a well-established aggression assay to investigate the role of DNA methylation in mediating aggressive responses to an intruder [4]. Although other stimuli such as alarm pheromone also elicit defensive behavior, we used the intruder assay because it provokes strong aggressive responses that are easily quantified, including biting and stinging. To profile DNA methylation signatures associated with aggression in the brain, we performed whole genome bisulfite sequencing (WGBS) of DNA from the brains of honey bees that responded aggressively to territorial intrusion, and compared the DNA methylation profiles with those of control bees who experienced an inert stimulus. DNA methylation was profiled at early and late time points (5 and 120 min, 6 individuals per phenotype) following the intrusion to investigate epigenetic differences between aggressive and control bees.

Here we show DNA methylation profile signatures of aggressive behavior in honey bees co-localize with gene regulatory elements in cis and are reflected in differences between aggressive Africanized bees and gentle European bees. DNA methylation differences increase over time from an early state immediately after interaction with an intruder to two hours post exposure, which is when memory consolidation typically occurs in other experimental systems [31]. These epigenetic profiles thus correlate with other findings that relate to both evolutionary and physiological time scales and provide insight into the molecular mechanisms of aggression.

Results

Effects of stimulus on methylation levels of individual cytosines over time

Out of 3,897,088 CpG sites with an average coverage of at least 5X, we identified 119,557 methylated CpG sites where at least 3 out of 24 samples had DNA methylation levels greater than 10%. As previously reported in honey bees, most of the DNA methylation resided within exons [25]. We found that 97,245 of the methylated CpG sites resided within gene bodies, with 69,831 overlapping exons. Additionally, 25,190 CpG sites were within promoters 2 kb upstream of transcriptional start sites, some of which also overlapped gene bodies.

Individual t-tests were performed on arcsine transformed methylation percentages and subject to multiple testing correction (FDR < 0.05) to compare aggressive and control bees within each time point and across time. While no individual CpGs were significantly different between aggressive and control bees 5 min after stimulus, 18 CpGs distinguished aggressive and control bees 120 min post stimulus (FDR ≤ 0.05). Most (16/18) of these CpGs were clustered on genomic scaffold 2.11. This cluster included two TRP channel genes, GB50805 and GB50806, that could play a role in sensing environmental stimuli [32]. Additionally, there were two CpGs differentially methylated between 5-min and 120-min control samples, and no differences between 5-min and 120-min aggressive samples. These results provide evidence for methylation changes in the brain in response to intruder exposure with clear differences 120 min post stimulus.

Regional analysis of DNA methylation reveals differences resulting from intruder exposure over time

Since the majority of differentially methylated CpGs co-localized to a region in scaffold 2.11, we next looked for additional regional methylation differences among clusters of CpGs and assessed significance at the regional level. This approach has been used to identify differences in cancer types [33], and it allowed us to directly compare aggressive individuals to control bees within each time point. We identified four differentially methylated regions (DMRs) between aggressive and control bees at the 5-min time point (Additional file 1: Table S1.). These DMRs could either indicate a very rapid epigenetic response to the intruder, or they might have existed prior to intruder exposure; akin to individual differences in social responsiveness previously reported in honey bees [34]. Some of these DMRs were located near or within genes that contain ion channel or receptor domains, with known neuronal functions in insects. One DMR within scaffold 8.9 is located over the 5′ end of NMDR3, an NMDA (n-methyl-D-aspartate) receptor, previously associated with increased aggression [35]. Other DMRs overlap with a homolog of CG42340, a regulator of synaptic plasticity in Drosophila [36].

At 120 min following intruder exposure, there were 50 DMRs between the aggressive and control bees (Additional file 2: Table S2.). This approximate ten-fold increase in the number of DMRs from 5 to 120 min, with no overlap between DMR sets, demonstrates sweeping temporal changes in the brain epigenetic response to an intruder. It is also consistent with the finding from Shpigler et al. (2017a), which performed a similar aggression study with honey bees and found the largest number of differentially expressed genes (DEGs) at 120 min after intruder exposure [4]. This DEG analysis also revealed no expression bias of DNA methyltransferase enzymes in either phenotype, which is reflected in the fact that we observe both hyper- and hypo- methylation within the 120 min DMRs. 120-min DMRs include most (11/16) CpGs from scaffold 2.11 that were identified independently, and confirm that this is a hotspot of epigenetic change, with 21 out of 50 DMRs residing on this scaffold (Fig. 1a). Increasing numbers of DMRs from 5 to 120 min also strengthens our finding of change over time within individual CpGs and strongly suggests that the DNA methylation changes are a direct response to intruder exposure.

Fig. 1
figure 1

Regional DNA Methylation Differences Distinguish Aggressive And Control Bees. Aggressive bees (red) and control bees (blue) show distinct DNA methylation profiles. a) scaffold 2.11 contains 21 of the 50 min-120 DMRs. The region depicted spans the ~ 90 kb hotspot of DNA methylation change and contains many signaling genes. Top panel shows two examples of DMRs within scaffold 2.11 with smoothed lines representing average methylation levels of individual methylated CpG’s (ignoring CpGs with < 10% methylation levels). Middle panel displays “hotspot” region of scaffold 2.11 containing multiple DMRs. Horizontal bars within yellow DMR areas show average methylation levels for each DMR. b-c) Top panel presents DNA methylation level for a short segment of scaffold 11.16, with dots for individual CpGs and smoothed lines for average methylation levels for each phenotype. The DMR in (b) is located over the DH31-R (GB47217) calcitonin receptor, and the DMR in (c) is over the DH44 (GB48796) diuretic hormone gene, both important for stress response

Genes associated with 120-min DMRs include those involved with RNA processing, chromatin remodeling and hormone signaling, processes that have been previously implicated in brain transcriptomic analyses of the bee aggression [4, 7]. DMR-associated chromatin remodeling genes included the tudor domain containing spoon and SET domain binding factor Sbf. There were also highly significant DMRs associated with DH44 (diuretic hormone 44) and DH31-R (calcitonin receptor), which are evolutionarily related to proteins in the mammalian CRF (corticotropin-releasing factor) and CGRP (calcitonin gene-related peptide) stress hormone signaling pathways (Fig. 1b-c). In vertebrates, CRF and its receptors have been identified as regulators of aggression and responses to social stress [13]. Overall, genes associated with 120-min DMRs may help facilitate large-scale gene expression response to intruder and alter signaling pathways that heighten the response to a new intruder in the future, because the presence of one intruder often presages the presence of more [37].

Correlation analysis utilizing individual cytosines suggests a time-dependent impact on gene expression and an enrichment of transcription factor binding sites

Since other forms of behavioral plasticity that involve neuronal remodeling, such as memory formation, require coordinated gene expression changes, we investigated the potential impact of DNA methylation on gene regulation across the genome. On the cellular level, DNA methylation is known to play a critical role in the late-phase of long-term potentiation (L-LTP, ~ 2 h post induction), when gene expression changes stabilize the strengthening of synapses [27]. To explore whether DNA methylation might work on a similar time scale in the context of aggression, we related the present results to previously published transcriptomic profiles in Shpigler et al. (2017a), wherein gene expression profiles were made at 30, 60, and 120 min after intruder exposure in the mushroom bodies (MB), rather than the whole brain [4] (The MB in honey bees are very prominent and constitute ca. 40% of total brain volume [38]). In this exploratory analysis, we checked the overlap of the DEGs from Shpigler et al. (2017a) with suggestive individual CpGs in the present study that passed a relaxed statistical threshold (t-test p-value < 0.05, but not subject to FDR correction). This relaxed threshold identified 3689 suggestive CpGs that were differentially methylated between aggressive and control bees 5 min post exposure, and 4323 CpGs 120 min post exposure. These suggestive CpGs in aggregate distinguish aggressive and control bees on a genome-wide level, but cannot be thought of in this way when considering any individual CpG. They were used to compare our findings to other genome-wide datasets. All conclusions made using suggestive CpGs were evaluated on the genome-wide level, not on the individual CpG level.

This analysis revealed that only 120-min DEGs have a significant overlap with suggestive CpGs (both P-value < 0.001 based on 1000 permutations, Additional file 3: Figure S1 a-f). Surprisingly, both 120-min and 5-min CpGs have significant overlaps with 120-min DEGs, often co-localizing in the same gene (Fig. 2, full GO annotation results in Additional file 4: Table S3). A total of 289 120-min DEGs overlapped 5-min CpGs and 319 120-min DEGs overlapped 120-min CpGs. We also observed a significant co-localization of 120-min and 5-min suggestive CpGs in 140 of the 120-min DEGs (overlap p = 1.1 × 10− 18 by hypergeometric test). Although 120-min DEGs overlapped both 5-min and 120-min CpGs, these CpG groups were largely distinct, with only 11 common CpGs within 120-min DEGs across time points. Together, these findings reveal a significant number of methylation differences marking 120-min DEGs, both immediately after and 2 h post exposure to an intruder. Existence of differential methylation before gene expression change might poise a set of genes ready to respond, much like bivalent epigenetic marks help poise genes in mammalian stem cells to respond to developmental cues [39]. Here, we found ~ 12% (140/1151) of 120-min DEGs marked by suggestive CpGs at 5 min were replaced by other suggestive CpGs at 120 min, potentially switching from a poised state to locking in the gene expression change post exposure.

Fig. 2
figure 2

Gene expression changes in response to intruder are marked by both early and late epigenetic differences. Gene expression differences between aggressive and control bees that arise 120 min after interaction with an intruder are often marked with differences in DNA methylation. There were 140 genes that had differential methylation both before and after gene expression changes occurred. For both (a) and (b), the top panel depicts boxplots of DNA methylation levels at both time points for aggressive and control bees. Asterisks indicate significant differences between aggressive and control bees within time points. Note that different sets of CpGs have significant differences at different time points, thus marking gene before and after gene expression change at 120 min with distinct CpGs. The bottom panel shows the differentially expressed gene and the location of the CpGs within the gene body. (c) Summary of overlap between 120-min DEGs and suggestive CpG’s. Both 5-min and 120-min suggestive CpGs have significant overlaps with120-minute DEGs, with a high degree of co-localization within a subset of genes. (d) GO enrichments for genes both differentially expressed and differentially methylated in response to intruder exposure

DMRs and suggestive CpGs overlap DNA methylation between Africanized and European honey bees, suggesting an evolutionary conserved role for epigenetics in regulating aggression

We explored whether the methylation differences we found to be associated with a real-time aggressive response were also associated with differences in aggression that occur on an evolutionary time scale. Evolutionary differences in colony-level aggression exist between European and African-derived subspecies of honey bees. Africanized honey bees (AHB) have a lower threshold needed to trigger an aggressive response and a larger proportion of aggressive individuals within a colony than do European honey bees (EHB), and Alaux et al. [5] found genes that differ in brain expression between AHB and EHB. Moreover, there was a statistically significant overlap between the genes that differed in expression between AHB and EHB and those that differed in EHB in response to exposure to alarm pheromone, which is associated with territorial intrusion and aggressive responses. We compared our within-colony (EHB) methylation differences to previously published methylation differences distinguishing AHB and EHB colonies [9]. We found a highly significant overlap between 120-min DMRs (and suggestive CpGs) and AHB / EHB differences (both P-value < 0.001 based on 1000 permutations, Additional file 3: Figure S1 g-h). This overlap suggests a conserved role for DNA methylation in modulating aggressive behavior across genetic backgrounds.

DNA methylation co-localizes with TFBS that regulate neuroplasticity and neurodevelopment

To further understand the molecular role of DNA methylation in gene regulation, we assessed the co-occupation of transcription factor binding sites (TFBS) and suggestive CpGs within promoters across the genome. DNA methylation can directly block TFBS in promoters, as seen in mammalian systems [40]. To identify transcription factors that might interact with DNA methylation to mediate the transcriptional response to an intruder, we calculated hypergeometric tests for the overlap of putative TFBSs and 120-min suggestive CpGs within promoters (2 kb) as defined by the OGS 3.2 gene set (Additional file 5: Table S4). Fourteen TFBSs had a highly significant overlap with suggestive CpGs (FDR of < 1.5 × 10− 10), and eleven of these TFBSs have some previously reported role in either neuroplasticity or neurodevelopment. Among these were Hr51 and foxo, both of which were identified as part of a brain transcriptional regulatory network constructed for honey bee aggressive behavior [4]. Hr51 and foxo, along with other enriched transcription factors lola and jigr1 (via interaction with jing) all contribute to neuron remodeling in Drosophila melanogaster by assisting either axon guidance or dendrite branching [36, 41,42,43]. Other enriched transcription factors help shape brain morphology during development and may play a role in adult behavior, including fkh, vvl/dfr, CG10267/Zif and Optix [44,45,46,47]. Because many genes overlap, we also present the enrichment for the non-overlapping portion of the promoter in Additional file 5: Table S4. Outside of the promoter, we observed strong enrichment of DNA methylation and TFBS within genes, potentially serving a protective role against spurious TF binding in gene bodies, as observed in termites [48]. Across the genome, we also observed that TFBS of Cf2, CG8281, jigr, lola, fkh, vvl, and CG10267/Zif were enriched over 120-min suggestive CpG sites by AME analysis [49]. For these TFBS enriched over suggestive CpGs, DNA methylation may be acting directly and altering the binding affinity for the TF. While largely unexplored in insects, DNA methylation might play a similar role in affecting transcription factor binding as in mammals [50].

Discussion

We present, to our knowledge, the first DNA methylation profiles of aggression in response to territorial intrusion. Our choice of honey bees as a model organism allowed us to directly study known aggressive individuals within a social context and their molecular response to an intruder. We found far more methylation differences between aggressive and control bees at 120 min than at 5 min after intruder exposure, which implies that DNA methylation could be assisting in altering neurogenomic states to maintain an aroused and vigilant state. Intruder-induced increases in arousal over the time scale used in this study are well known in the animal world and have been documented in honey bees [4, 5]. The observed brain epigenetic differences had a significant overlap with both gene expression differences in response to an intruder and with DNA methylation differences across genetic strains of honey bees that differ drastically in response to intruders.

Honey bees exhibit an increased sensitivity to future territorial intrusions after an initial interaction with an intruder in the colony [3, 4]. A change in response threshold to subsequent intruder interactions likely requires neuronal remodeling and has been shown to involve gene expression changes in genes related to steroid hormone receptor activity and chromosome organization [4]. Epigenetic mechanisms provide a flexible means of gene regulation that can result in temporary gene expression changes and contribute to neural plasticity. This hypothesized link between epigenetic regulation and neural plasticity is further strengthened by the highly significant overlap we found between DNA methylation differences and binding sites for neurobiologically relevant transcription factors. This evidence suggests a mechanism whereby DNA methylation can affect brain gene regulation over a very short time period by blocking or enhancing the action of transcription factors.

We also provide evidence for DNA methylation marking genes prior to expression changes, thus potentially poising genes in a ready state that can more immediately respond to an encounter with an intruder. Additionally, we previously reported differences in histone H3 lysine 27 acetylation (H3K27ac), a marker of open chromatin, in mushroom bodies at the same time point after intruder exposure [4], further supporting the conclusion that aggression leads to epigenetic reconfiguration in the brain. Interestingly, neither our DMRs nor suggestive CpGs overlapped regions of differential H3K27ac from Shpigler et al. (2017a), possibly suggesting a different role for these two epigenetic modifications. Lastly, intruder exposure also rapidly alters H3K27ac levels in the diencephalon of freshwater stickleback fish, suggesting an evolutionarily conserved role for epigenetic response to and intruder [51].

There is growing interest in exploring the link between epigenetics, early life stressors and altered behavior. In humans, multiple genes including IL-6 interleukin 6 and the serotonin transporter SLC6A4 (solute carrier family 6 member 4) have altered DNA methylation levels in their promoter regions (reviewed in Waltes et al. [52]) thought to be caused by early life trauma. Model organisms such as rats and chickens have provided direct evidence of artificially elevated levels of cortisol or diminished maternal care early in life leading to altered epigenetic levels of key genes involved in aggressive behavior. Our study adds to this growing literature especially with the finding of dynamic changes in brain methylation associated with aggression-related genes after exposure to an intruder. These changes suggest epigenetic involvement in altering future aggression-related behavior.

Conclusions

Behavioral plasticity allows organisms to adapt to new adverse circumstances, including hostile territorial intrusions. We showed that differences exist between the brain DNA methylation profiles of aggressive and control honey bees, indicating a possible epigenetic basis for plasticity in aggression. Most of the differences were observed two hours after exposure to a foreign intruder but not at an earlier time point, suggesting that they arose as a direct result of experience. Additionally, we detected a few methylation differences between aggressive and control bees five minutes after intruder exposure, possibly reflecting pre-existing individual differences in the epigenetic control of aggressive tendencies. This study provides the first evidence that changes in brain DNA methylation occur in response to territorial intrusion, supporting the conclusion that there is an epigenetic basis to behavioral plasticity in aggression.

Methods

Laboratory intruder assay and bee collection

Honey bee workers used in this study were derived from a queen instrumentally inseminated with semen from a single drone and reared up to adulthood in a colony maintained under typical conditions at the University of Illinois at Urbana-Champaign Bee Research Facility, Urbana, Illinois. Bees in this area are a mixture of European subspecies of Apis mellifera, primarily ligustica. Honeycomb frames containing older pupae were collected from a single colony and maintained in an incubator room at 34 °C, 50 ± 5% humidity. Emerging adults were collected every 24 h and marked on the thorax with Testor’s paint. The bees were placed in groups of ten in vertically oriented petri dishes lined with a layer of beeswax foundation [4]. They were fed 30% sucrose solution and pollen balls made of fresh frozen pollen mixed with 30% sucrose solution. Food was checked daily and replenished as necessary. The groups were kept in the incubator room for 8 days to establish group identity prior to the intruder assay.

Aggression was assessed using a laboratory intruder assay modified from [4], which includes a supplemental video of the assay. A foreign worker bee from another colony was introduced into each petri dish as an intruder, and biting and stinging responses to the intruder were observed and recorded for the next 5 min. The colored paint dots on the bees allowed for assignment of each aggressive event to an individual bee. The two most aggressive bees, defined as those displaying the highest number of biting and stinging events, were collected from each group and flash frozen in liquid nitrogen. Collections were done either 5 min or 120 min following initial intruder exposure. Control groups were exposed to a 0.2 ml plastic tube similar in size to a worker bee, which was removed after 5 min. Two bees were collected from each control group at random either 5 min or 120 min following initial exposure to the tube. After flash freezing, all collected bees were transferred individually to microfuge tubes and stored at − 80 °C.

Sample preparation and whole-genome bisulfite sequencing

The heads were removed from frozen bees and placed into a dry ice/100% ethanol bath. The frons were removed with a scalpel and the heads were completely submerged in RNALater ICE overnight at − 20 °C, allowing the reagent to permeate the brain tissue for 15 to 16 h. Whole brains were dissected from head capsules at room temperature using a scalpel and fine forceps, with care taken to remove the adjacent hypopharyngeal gland. Individual brains were placed in microfuge tubes and stored at − 80 °C until DNA extraction.

Genomic DNA was extracted from whole brains using the Gentra Puregene Tissue Kit (Qiagen) and stored at − 20 °C. Unmethylated lambda phage DNA (Promega) was spiked in as a control for assessment of bisulfite conversion efficiency. The samples were given to the Roy J. Carver Biotechnology Center at the University of Illinois for quality control, library construction, and sequencing. The integrity of the samples was checked using the AATI Fragment Analyzer, which confirmed that the DNA was of high molecular weight. Shotgun DNA libraries were made using the KAPA Library Preparation Kit, then bisulfite treated using the Zymo EZ DNA Methylation-Lightning Kit. The libraries were then amplified for 10 PCR cycles using HiFi Uracil+ DNA Polymerase (KAPA Biosystems), pooled in equimolar quantities, and sequenced on the HiSeq 2500 (Illumina). Paired-end 100-nucleotide sequence reads were generated from six control and six experimental samples for each collection time point (5 min and 120 min). FASTQ files were generated and demultiplexed with the bcl2fastq v2.17.1.14 Conversion Software (Illumina).

Alignment

Paired-end reads were trimmed with the trimmomatic program (v 0.33) using the following parameters: (ILLUMINACLIP:TruSeq3-PE-2.fa:2:15:10 CROP:98 HEADCROP:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:30). Trimmed reads were aligned to the Amel_4.5 honey bee genome assembly with bowtie2 within the bismark program (v 0.16.1) using bowtie2 default settings except for max insert size set to 1000. The bismark program was also used to mark duplicates and extract DNA methylation levels for each CpG. The R package bsseq was used to merge DNA methylation information across all samples and calculate coverage. Bisulfite conversion efficiency was estimated both empirically by assessing number of unconverted cytosines that aligned to the lambda genome from the lambda spike-in control and assessment of typically unmethylated cyotsines across the genome using the MethPipe software [53] (Details: Additional file 6: Table S5).

Differential methylation analysis at individual cytosines

Cytosines that had an average coverage across all samples of 5 or more reads and displayed at least 10% methylation levels in three or more individual bees were defined as methylated (119,557 met this criteria). We chose a minimum of three individuals to avoid spurious methylation in one or two samples while at the same time allow for CpG sites unique to sub-phenotypes that might only be represented by a small subset of bees. To find differences in DNA methylation between aggressive and control bees over time, we first arcsine transformed percentage measurements derived from ratios of reads containing methylated and unmethylated cytosines in the CpG dinucleotide context. Student’s t-test was performed on individual CpG sites and these raw p-values were used to identify suggestive CpGs (p-value ≤0.05) used in subsequent methods. Multiple test correction was performed using the qvalue package in R [54, 55] and a FDR was calculated for each methylated CpG.

Identifying differentially methylated regions based on individual cytosines

To determine whether aggression is associated with changes in methylation at the regional level, we implemented a pipeline for the identification of differentially methylated regions (DMRs). Raw p-values for individual CpG sites within each region designated by the bumphunter package were combined using the comb-p software [56] and regional p-values were corrected for multiple testing. The comb-p package first calculates the correlation between proximal P-values of varying distances. The Stouffer–Liptak–Kechris correction (slk) is then applied which corrects each P-value respective of weighting determined by correlation calculation. This analysis identified regions of coordinated change where multiple CpG sites change in the same direction and pairwise comparisons were made between each combination of time point and stimulus (n = 6 for each group).

Overlap of differential methylation and genomic features / published data

To assess the genome-wide impact of brain DNA methylation on brain gene expression, we used gene expression data from the Shpigler et al. (2017a) [4] study that utilized the same experimental set-up to identify aggressive bees. This study identified brain gene expression differences between aggressive and control bees at 30, 60, and 120 min after the introduction of an intruder bee to a petri dish containing test bees. To identify trends, we included all CpGs with a p-value ≤0.05 in our analysis (suggestive CpGs), which includes 3689 differences at minute-5 and 4323 at minute-120. The suggestive CpGs identified at each time point were largely distinct, with only 135 common CpGs across time points. We quantified the number of DEGs that overlap suggestive CpGs for each time point. To test the significance of this overlap, 1000 permutations of random sets of genes matching the number of significant DEGs were created for each time point and checked for overlap with suggestive CpGs. Empirical p-values for overlap were reported based on permutation results. Significance of overlap between 5-min and 120-min suggestive CpGs within 120-min DEGs was calculated by hypergeometric test.

A similar analysis was performed on DNA methylation data from Cingolani et al. [9], however, WGBS data was originally aligned to the Amel2 version of the genome, which necessitated alignment to the Amel4.5 version of the genome to compare results. Permutation analysis was used to calculate the degree of overlap between DMRs and suggestive CpGs found in our study to CpGs that differed between Africanized and European honey bees by 30%.

Co-localization of DNA methylation and TFBS in promoters

Locations of putative transcription factor binding sites (TFBS) were found by mapping TFBS motifs identified in Drosophila melanogaster to the Amel 4.5 version of the honey bee genome using the program CisGenome [40]. Suggestive CpGs identified at minute-120 and TFBS were quantified within each promoter (2 kb upstream of TSS) for genes of OGS 3.2. Hypergeometric tests were performed to calculate p-values for overlap between DNA methylation and TFBS within promoters, and FDRs were calculated from p-values using the qvalue package in R [54, 55]. Additional analysis of genome-wide enrichment of TFBS around suggestive CpGs irrespective of gene location was performed using the AME program [49], which is part of the MEME suite. For each suggestive CpG identified at 120-min post intruder exposure, 20 base pairs of sequence up and down stream of cytosine was recorded (total 41 bp for each suggestive CpG). Background sets were based on the sequence surrounding 4000 randomly selected methylated CpGs (also 41 bp per cytosine). Target sequences were analyzed for enrichment of TFBS (Combined Drosophila Databases) against background sets.

Abbreviations

5mC:

5-methylcytosine

AHB:

Africanized honey bee

Cf2:

Chorion factor 2

CGRP:

calcitonin gene-related peptide

CpG:

Cytosine-guanine dinucleotide

CRF:

corticotropin-releasing factor

DEG:

Differentially expressed gene

Dfr:

Drifter gene

DH31-R:

Diuretic hormone 31 receptor

DH44:

Diuretic hormone 44

DMR:

Differentially methylated region

DNA:

Deoxyribonucleic acid

DNMT:

DNA methyltransferase

EHB:

European honey bee

FDR:

False discovery rate

Fkh:

Forkhead gene

foxo:

forkhead box, sub-group O

H3K27ac:

histone H3 lysine 27 acetylation

Hr51:

Hormone receptor 51 gene

IL-6:

Interleukin 6

Jigr1:

Jing interacting regulatory gene 1

Kb:

Kilobase

L-LTP:

Late phase of long-term potentiation

Lola:

longitudinals lacking gene

MB:

Mushroom body

NMDA:

N-methyl-D-aspartate

NMDR3:

N-methyl-D-aspartate receptor 3

Sbf:

SET domain binding factor

SET:

Su(var)3–9, Enhancer-of-zeste and Trithorax

SLC6A4:

solute carrier family 6 member 4

Tet:

Ten-Eleven Translocation

TFBS:

Transcription factor binding site

TRP channels:

Transient receptor potential channels

Vvl:

Ventral veins lacking gene

WGBS:

Whole-genome bisulfite sequencing

Zif:

Zinc finger protein

References

  1. Sweatt JD. Neural plasticity and behavior – sixty years of conceptual advances. J Neurochem. 2016;139:179–99.

    Article  CAS  PubMed  Google Scholar 

  2. Sandi C, Haller J. Stress and the social brain: behavioural effects and neurobiological mechanisms. Nat Rev Neurosci. 2015;16:290–304.

    Article  CAS  PubMed  Google Scholar 

  3. Alaux C, Robinson GE. Alarm pheromone induces immediate-early gene expression and slow behavioral response in honey bees. J Chem Ecol. 2007;33:1346–50.

    Article  CAS  PubMed  Google Scholar 

  4. Shpigler HY, Saul MC, Murdoch EE, Cash-Ahmed AC, Seward CH, Sloofman L, et al. Behavioral, transcriptomic and epigenetic responses to social challenge in honey bees. Genes Brain Behav. 2017;16(6):579-91.

  5. Alaux C, Sinha S, Hasadsri L, Hunt GJ, Guzmán-Novoa E, DeGrandi-Hoffman G, et al. Honey bee aggression supports a link between gene regulation and behavioral evolution. Proc Natl Acad Sci U S A. 2009;106:15400–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Li-Byarlay H, Rittschof CC, Massey JH, Pittendrigh BR, Robinson GE. Socially responsive effects of brain oxidative metabolism on aggression. Proc Natl Acad Sci. 2014;111:12533–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Rittschof CC, Bukhari SA, Sloofman LG, Troy JM, Caetano-Anollés D, Cash-Ahmed A, et al. Neuromolecular responses to social challenge: common mechanisms across mouse, stickleback fish, and honey bee. Proc Natl Acad Sci. 2014;111:17929–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Rittschof CC, Robinson GE. Manipulation of colony environment modulates honey bee aggression and brain gene expression. Genes Brain Behav. 2013;12:802–11.

    Article  CAS  PubMed  Google Scholar 

  9. Cingolani P, Cao X, Khetani RS, Chen C-C, Coon M, Sammak A, et al. Intronic non-CG DNA hydroxymethylation and alternative mRNA splicing in honey bees. BMC Genomics. 2013;14:666.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Chandrasekaran S, Ament SA, Eddy JA, Rodriguez-Zas SL, Schatz BR, Price ND, et al. Behavior-specific changes in transcriptional modules lead to distinct and predictable neurogenomic states. Proc Natl Acad Sci. 2011;108:18020–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Wu Y, Patchev AV, Daniel G, Almeida OFX, Spengler D. Early-life stress reduces DNA methylation of the Pomc gene in male mice. Endocrinology. 2014;155:1751–62.

    Article  PubMed  Google Scholar 

  12. Murgatroyd C, Patchev AV, Wu Y, Micale V, Bockmühl Y, Fischer D, et al. Dynamic DNA methylation programs persistent adverse effects of early-life stress. Nat Neurosci. 2009;12:1559–66.

    Article  CAS  PubMed  Google Scholar 

  13. Elliott E, Ezra-Nevo G, Regev L, Neufeld-Cohen A, Chen A. Resilience to social stress coincides with functional DNA methylation of the Crf gene in adult mice. Nat Neurosci. 2010;13:1351–3.

    Article  CAS  PubMed  Google Scholar 

  14. Tsankova NM, Berton O, Renthal W, Kumar A, Neve RL, Nestler EJ. Sustained hippocampal chromatin regulation in a mouse model of depression and antidepressant action. Nat Neurosci. 2006;9:519–25.

    Article  CAS  PubMed  Google Scholar 

  15. Day JJ, Sweatt JD. DNA methylation and memory formation. Nat Neurosci. 2010;13:1319–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Heyward FD, Sweatt JD. DNA methylation in memory formation: emerging insights. Neurosci Rev J Bringing Neurobiol Neurol Psychiatry. 2015;21:475–89.

    CAS  Google Scholar 

  17. Miller CA, Gavin CF, White JA, Parrish RR, Honasoge A, Yancey CR, et al. Cortical DNA methylation maintains remote memory. Nat Neurosci. 2010;13:664–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Halder R, Hennion M, Vidal RO, Shomroni O, Rahman R-U, Rajput A, et al. DNA methylation changes in plasticity genes accompany the formation and maintenance of memory. Nat Neurosci. 2016;19:102–10.

    Article  CAS  PubMed  Google Scholar 

  19. Saunderson EA, Spiers H, Mifsud KR, Gutierrez-Mecinas M, Trollope AF, Shaikh A, et al. Stress-induced gene expression and behavior are controlled by DNA methylation and methyl donor availability in the dentate gyrus. Proc Natl Acad Sci U S A. 2016;113:4830–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Meadows JP, Guzman-Karlsson MC, Phillips S, Holleman C, Posey JL, Day JJ, et al. DNA methylation regulates neuronal glutamatergic synaptic scaling. Sci Signal. 2015;8:ra61-ra61.

  21. Gavin DP, Kusumo H, Sharma RP, Guizzetti M, Guidotti A, Pandey SC. Gadd45b and N-methyl-D-aspartate induced DNA demethylation in postmitotic neurons. Epigenomics. 2015;7:567–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Day JJ, Childs D, Guzman-Karlsson MC, Kibe M, Moulden J, Song E, et al. DNA methylation regulates associative reward learning. Nat Neurosci. 2013;16:1445–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Wang Y, Jorda M, Jones PL, Maleszka R, Ling X, Robertson HM, et al. Functional CpG methylation system in a social insect. Science. 2006;314:645–7.

    Article  CAS  PubMed  Google Scholar 

  24. Herb BR, Wolschin F, Hansen KD, Aryee MJ, Langmead B, Irizarry R, et al. Reversible switching between epigenetic states in honeybee behavioral subcastes. Nat Neurosci. 2012;15:1371–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Lyko F, Foret S, Kucharski R, Wolf S, Falckenhayn C, Maleszka R. The honey bee epigenomes: differential methylation of brain DNA in queens and workers. PLoS Biol. 2010;8:e1000506.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Foret S, Kucharski R, Pellegrini M, Feng S, Jacobsen SE, Robinson GE, et al. DNA methylation dynamics, metabolic fluxes, gene splicing, and alternative phenotypes in honey bees. Proc Natl Acad Sci. 2012;109:4968–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Kucharski R, Maleszka J, Foret S, Maleszka R. Nutritional control of reproductive status in honeybees via DNA methylation. Science. 2008;319:1827–30.

    Article  CAS  PubMed  Google Scholar 

  28. Cardoso-Júnior CAM, Guidugli-Lazzarini KR, Hartfelder K. DNA methylation affects the lifespan of honey bee (Apis mellifera L.) workers - evidence for a regulatory module that involves vitellogenin expression but is independent of juvenile hormone function. Insect Biochem Mol Biol. 2017;92:21–9.

    Article  PubMed  Google Scholar 

  29. Biergans SD, Giovanni Galizia C, Reinhard J, Claudianos C. Dnmts and Tet target memory-associated genes after appetitive olfactory training in honey bees. Sci Rep. 2015;5:16223.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Li-Byarlay H, Li Y, Stroud H, Feng S, Newman TC, Kaneda M, et al. RNA interference knockdown of DNA methyl-transferase 3 affects gene alternative splicing in the honey bee. Proc Natl Acad Sci U S A. 2013;110:12750–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Kennedy AJ, Sweatt JD. Drugging the methylome: DNA methylation and memory. Crit Rev Biochem Mol Biol. 2016;51:185–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Zheng J. Molecular mechanism of TRP channels. Compr Physiol. 2013;3:221–42.

    PubMed  PubMed Central  Google Scholar 

  33. Li X, Liu Y, Salz T, Hansen KD, Feinberg A. Whole-genome analysis of the methylome and hydroxymethylome in normal and malignant lung and liver. Genome Res. 2016;26(12):1730-41.

  34. Shpigler HY, Saul MC, Corona F, Block L, Cash Ahmed A, Zhao SD, et al. Deep evolutionary conservation of autism-related genes. Proc Natl Acad Sci U S A. 2017;114:9653–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Duncan GE, Inada K, Farrington JS, Koller BH, Moy SS. Neural activation deficits in a mouse genetic model of NMDA receptor hypofunction in tests of social aggression and swim stress. Brain Res. 2009;1265:186–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Lin S, Huang Y, Lee T. Nuclear receptor unfulfilled regulates axonal guidance and cell identity of Drosophila mushroom body neurons. PLoS One. 2009;4:e8392.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Zayed A, Robinson GE. Understanding the relationship between brain gene expression and social behavior: lessons from the honey bee. Annu Rev Genet. 2012;46:591–615.

    Article  CAS  PubMed  Google Scholar 

  38. Howse PE. Design and function in the insect brain. In: Barton Browne L. (eds). Exp Anal Insect Behav. Berlin: Springer; 1974.

  39. Minoux M, Holwerda S, Vitobello A, Kitazawa T, Kohler H, Stadler MB, et al. Gene bivalency at Polycomb domains regulates cranial neural crest positional identity. Science. 2017;355:eaal2913.

    Article  PubMed  Google Scholar 

  40. Ji H, Vokes SA, Wong WH. A comparative analysis of genome-wide chromatin immunoprecipitation data for mammalian transcription factors. Nucleic Acids Res. 2006;34:e146.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Sears JC, Broihier HT. FoxO regulates microtubule dynamics and polarity to promote dendrite branching in Drosophila sensory neurons. Dev Biol. 2016;418:40–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Giniger E, Tietje K, Jan LY, Jan YN. Lola encodes a putative transcription factor required for axon growth and guidance in Drosophila. Dev Camb Engl. 1994;120:1385–98.

    CAS  Google Scholar 

  43. Sun X, Morozova T, Sonnenfeld M. Glial and neuronal functions of the Drosophila homolog of the human SWI/SNF gene ATR-X (DATR-X) and the jing zinc-finger gene specify the lateral positioning of longitudinal glia and axons. Genetics. 2006;173:1397–415.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Meier S, Sprecher SG, Reichert H, Hirth F. Ventral veins lacking is required for specification of the tritocerebrum in embryonic brain development of Drosophila. Mech Dev. 2006;123:76–83.

    Article  CAS  PubMed  Google Scholar 

  45. Chang KC, Garcia-Alvarez G, Somers G, Sousa-Nunes R, Rossi F, Lee YY, et al. Interplay between the transcription factor Zif and aPKC regulates neuroblast polarity and self-renewal. Dev Cell. 2010;19:778–85.

    Article  CAS  PubMed  Google Scholar 

  46. Gold KS, Brand AH. Optix defines a neuroepithelial compartment in the optic lobe of the Drosophila brain. Neural Develop. 2014;9:18.

    Article  Google Scholar 

  47. Miguel-Aliaga I, Thor S, Gould AP. Postmitotic specification of Drosophila Insulinergic Neurons from Pioneer neurons. PLoS Biol. 2008;6:e58.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Glastad KM, Hunt BG, Yi SV, Goodisman MA. DNA methylation in insects: on the brink of the epigenomic era. Insect Mol Biol. 2011;20:553–65.

    Article  CAS  PubMed  Google Scholar 

  49. McLeay RC, Bailey TL. Motif enrichment analysis: a unified framework and an evaluation on ChIP data. BMC Bioinformatics. 2010;11:165.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Yin Y, Morgunova E, Jolma A, Kaasinen E, Sahu B, Khund-Sayeed S, et al. Impact of cytosine methylation on DNA binding specificities of human transcription factors. Science. 2017;356:eaaj2239.

    Article  PubMed  Google Scholar 

  51. Bukhari SA, Saul MC, Seward CH, Zhang H, Bensky M, James N, et al. Temporal dynamics of neurogenomic plasticity in response to social interactions in male threespined sticklebacks. PLoS Genet. 2017;13:e1006840.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Waltes R, Chiocchetti AG, Freitag CM. The neurobiological basis of human aggression: a review on genetic and epigenetic mechanisms. Am J Med Genet Part B Neuropsychiatr Genet Off Publ Int Soc Psychiatr Genet. 2016;171:650–75.

    Article  Google Scholar 

  53. Song Q, Decato B, Hong EE, Zhou M, Fang F, Qu J, et al. A reference methylome database and analysis pipeline to facilitate integrative and comparative epigenomics. PLoS One. 2013;8:e81148.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci. 2003;100:9440–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57:289–300.

    Google Scholar 

  56. Pedersen BS, Schwartz DA, Yang IV, Kechris KJ. Comb-p: software for combining, analyzing, grouping and correcting spatially correlated P-values. Bioinforma Oxf Engl. 2012;28:2986–8.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We thank Terry Harrison for maintenance of honey bee colonies; Aditi Warhekar for assistance with behavioral assays and sample collection; Amy Cash Ahmed for preparation of some of the DNA samples; Alvaro Hernandez and the staff of the High-Throughput Sequencing Unit at the Roy J. Carver Biotechnology Center for oversight of library preparation and sequencing; and members of the Robinson lab for helpful feedback on the manuscript.

Funding

This project was supported by National Institutes of Health grant R21-MH107962 to G.E.R. Additonally, M.S.S. was supported by grant SFLife 291812 (principal investigators Lisa Stubbs and G.E.R.) from the Simons Foundation.

Availability of data and materials

The datasets generated and analyzed during the current study are publicly available in NCBI GEO repository under GSE104994.

Author information

Authors and Affiliations

Authors

Contributions

M.S.S. and G.E.R. jointly conceived the study and designed the experiments. M.S.S. performed behavioral assays and prepared biological samples. B.R.H. did most of the data analysis, with contributions from M.S.S and C.J.F., who also developed analysis pipelines for WGBS data. M.S.S., B.R.H. and G.E.R. prepared the manuscript. All authors read and approved the final paper.

Corresponding author

Correspondence to Gene E. Robinson.

Ethics declarations

Ethics approval and consent to participate

Honey bees are not a regulated invertebrate. Ethical use approval is not required for this species.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. DMRs between control and aggressive bees at 5 min post interaction. (XLSX 36 kb)

Additional file 2:

Table S2. DMRs between control and aggressive bees at 120 min post interaction. (XLSX 61 kb)

Additional file 3:

Figure S1. Results of permutation testing for DEG/CpG overlap and comparison with Cingolani et al. 2013 [9]. (a-f) Minute 5 and Minute 120 suggestive CpGs were overlapped with gene expression differences of Aggressive and Control bees from Shpigler et al. (2017a) [4]. To test the significance of the overlap, we chose at random the same number of significantly different DEGs, checked the overlap and repeated this process 1000 times. For each plot, the blue histogram represents the distribution of random trials and the red vertical line represents the true number of overlapping DEGs with suggestive CpGs. Note that both minute 5 and minute 120 suggestive CpGs had a highly significant overlap with only minute 120 DEGs. (g-h) A similar approach was taken to test the significance of the overlap between AHB / EHB differences and minute 120 DMRs and suggestive CpGs. In both cases the true overlap was greater than the 1000 permutations. (PDF 478 kb)

Additional file 4:

Table S3. Details of GO enrichments for Differentially expressed genes from Shpigler et al. 2017a [34] that overlap suggestive CpGs. (XLSX 47 kb)

Additional file 5:

Table S4. Enrichments of Minute-120 suggestive CpGs and transcription factor binding sites in gene regions. (XLSX 190 kb)

Additional file 6:

Table S5. Summary of samples and WGBS sequencing. (XLSX 43 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Herb, B.R., Shook, M.S., Fields, C.J. et al. Defense against territorial intrusion is associated with DNA methylation changes in the honey bee brain. BMC Genomics 19, 216 (2018). https://doi.org/10.1186/s12864-018-4594-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-4594-0

Keywords