Skip to main content

The effect of DNA methylation on bumblebee colony development



Although around 1% of cytosines in bees’ genomes are known to be methylated, less is known about methylation’s effect on bee behavior and fitness. Chemically altered DNA methylation levels have shown clear changes in the dominance and reproductive behavior of workers in queen-less colonies, but the global effect of DNA methylation on caste determination and colony development remains unclear, mainly because of difficulties in controlling for genetic differences among experimental subjects in the parental line. Here, we investigated the effect of the methylation altering agent decitabine on the developmental rate of full bumblebee colonies. Whole genome bisulfite sequencing was used to assess differences in methylation status.


Our results showed fewer methylated loci in the control group. A total of 22 CpG loci were identified as significantly differentially methylated between treated and control workers with a change in methylation levels of 10% or more. Loci that were methylated differentially between groups participated in pathways including neuron function, oocyte regulation and metabolic processes. Treated colonies tended to develop faster, and therefore more workers were found at a given developmental stage. However, male production followed the opposite trend and it tended to be higher in control colonies.


Overall, our results indicate that altered methylation patterns resulted in an improved cooperation between workers, while there were no signs of abnormal worker dominance or caste determination.


The phenotype of an individual is determined, ultimately, by the context-driven interpretation of a given DNA sequence [1]. DNA methylation, of which the most common example is the addition of a methyl group to a cytosine, is a reversible biological process that can change the activity of a DNA segment {Citation} [2] and lead to phenotypic plasticity in modular organisms such as plants [3], or induce reversible and quick adaptation to stress conditions in unicellular fungi [4]. In social insects, where genetic relatedness affects cooperation and reproductive behavior within colonies [5], DNA methylation has been shown to affect important phenotypic features such as caste differentiation and worker reproductive behavior [6, 7].

Genome-wide methylation levels are highly variable and differ substantially between major taxonomical groups. In mammals for example, about 70% of CpG dinucleotides are methylated in somatic cells [8], while the genome of most plants, invertebrates, fungi, or protists shows “mosaic” methylation patterns, where only specific genomic elements are targeted and distinct patterns of methylated and unmethylated domains can be discerned [9, 10]. In insects, DNA methylation levels are often low (less than 1%), but they generally concentrate in gene-coding regions [11]. While DNA methylation at gene promoter regions suggests gene silencing as the main function, as proposed for mammals [12], methylation within gene-coding regions suggests a role in alternative splicing [13]. However, in the particular case of insects such as bumblebees, differentially expressed genes contain lower levels of methylation compared to non-differentially expressed genes, which further indicates that DNA methylation in bumblebees is more related to gene expression than to an alteration of gene function [14, 15].

One of the most defining characteristic of insect societies is the reproductive division of labor, where workers usually do not produce offspring in the presence of a queen [16]. The prevailing theory that explains such ‘altruistic’ behavior is inclusive fitness [17, 18], which states that workers will be selected to nurse their mother’s offspring rather than investing time and energy in their own progeny. In bumblebees, altruistic worker behavior can be explained by the higher relatedness of the female sisters, who would share ¾ of their genome, compared to the relatedness of the workers to their own male progeny. The annual colony life cycle of a bumblebee is therefore divided into a cooperative phase, when female workers are produced while the queen has absolute reproductive dominance, and a highly aggressive competition phase later in the season when the workers and queen compete over male production [19]. Bumblebee queens mate only once, implying that all workers are full sisters and that their genomes therefore have all the necessary information to become “dominant” and show reproductive behavior. However, if the queen dies or is removed, unmated workers can differentiate into reproductive and non-reproductive sub-castes by both their ovary development and aggressive behavior [20].

The genomes of queen-less reproductive workers and queen-less non-reproductive workers have been shown to differ in methylation levels [7], suggesting that worker reproductive behavior may be determined, and inherited, by epigenetic factors. Moreover, queen-less workers whose genomes had been experimentally altered by using an inhibitor of DNA methyl-transferase were more aggressive and more likely to develop ovaries compared with control queen-less workers [7], indicating that DNA methylation is important in this highly plastic reproductive division of labor. However, these results also indicate that variation in DNA methylation levels could affect overall colony development: if workers with low DNA methylation levels gain dominance and are more keen to reproduce, the resulting colonies would show an earlier disruption of the cooperation phase, which in turn would lead to a higher production of males and a decreased production of queens [19], but see [21]. However, the reproductive behavior of bumblebee workers, when separated from the queen under laboratory conditions, can be expected to differ from those in colonies where the founder queen remains present. Most experiments studying epigenetic effects on worker behavior used micro-colonies of full sisters that were obtained and kept separately from single mated Bombus queens [22]. While this ensures that all individuals have the same genotype, the absence of the queen is a rather unnatural situation that may have a profound impact on the results.

Caste determination is another feature of social insects that can be affected by epigenetic factors. Previous work on honeybees has shown that changes in methylation levels are involved in the switch between workers and queens [23]. The comparison of larval heads between queens and workers of honeybees show a total of 2399 genes with significant differences of methylation [24]. In addition, substantial differentially methylated genes were found among different castes in the termite Zootermopsis nevadensis [25] and the ant Camponotus floridanus [26]. However, no association between caste and methylation has been found in some primitive wasps, such as Polistes spp. ([27]). Recently, [15] found differences in methylation levels between reproductive castes of bumblebee workers, with some differentially methylated genes involved in behavior and reproductive processes. Their results also showed high inter-colony variance in methylation levels, suggesting that different couples of queens and males transmit different methylomes to their progeny [28,29,30], which in turn will lead to developmental differences at the colony level [31]. DNA methylation could also be involved in worker vs. (daughter) queen development by fertilized eggs in bumblebees.

In this study, we tested the hypothesis that DNA methylation had a significant effect on colony development of the bumblebee Bombus terrestris. By experimentally exposing a pure genetic line of B. terrestris founders to the methylation disruptor decitabine through sugar water provisions, we first investigated the effect of DNA methylation on the developmental fate of larvae, and how this affected colony development. Second, to identify the specific loci that were affected by the addition of decitabine and the biological pathways these genes were involved in, brain tissue samples were collected from adult workers to be examined for DNA methylation at single base resolution using whole genome bisulfite sequencing (WGBS).


Temporal succession of main colony events

We obtained 6/6 and 5/6 developed colonies in control and treated queens, respectively. The time needed to lay the first eggs or to produce the first pupae did not differ significantly between founder queens that were assigned to the different treatment levels (W = 15.0, P = 1; W = 12.5, P = 0.7112 for eggs and pupae, respectively). Correspondingly, the number of days needed to produce the first workers in a colony also did not differ significantly between treatments (W = 12.5, P = 0.7138), nor did date of first male appearance (W = 14, P = 0.926).

Colony size and production of males

Overall, colonies supplemented with decitabine had a significantly higher brood size than control colonies (χ2 = 23.36, P < 0.001, Fig. 1a). In agreement with the end of the cooperative phase of the colony, differences were more pronounced when colonies were counted 8 weeks after colony start-up (Z = − 4.33, P < 0.0001), although they remained significant when colonies were examined two weeks later (Z = − 2.60, P = 0.046, Fig. 1a). Larval mortality did not differ among control and treated colonies (χ2 = 0.14, P = 0.705). However, queens treated with decitabine were more active at laying eggs before the competition point, as indicated by the significantly higher number of egg cups in treated colonies at week 8 (χ2 = 24.12, P < 0.001). Correspondingly, treatment also positively affected worker production (χ2 = 16.25, P < 0.001, Fig. 1b). This effect was consistent at both assessment weeks, although for this parameter developmental differences accumulated over time, yielding significant differences for worker production for the last counting week only (Fig. 1b, z = − 1.871, P = 0.240, and Z = − 3.497, P = 0.003 for weeks 8 and 10 after colony start-up, respectively).

Fig. 1

(a) Average brood size (sum of egg cups and larvae), (b) average number of workers and (c) average number of males per colony for control and decitabine-treated colonies, counted 8 and 10 weeks after start-up. Symbol depict ls-means, plotted at the original scale, and vertical lines show 95% confidence intervals

The administration of decitabine tended to decrease the number of males produced after the competition point (week 10, Fig. 1c). However, the overall effect of the treatment on male production did not reach statistical significance (χ2 = 1.614, P = 0.204).

Worker reproductive behaviour

Random dissections at week 8 showed developed ovaries as the most observed stage of worker reproductive status at that point of colony development, and therefore it was observed in 7/10 (control) and 8/10 of the dissected workers (treated colonies). Remaining specimens showed either incipient (1/10 of treated colonies, 0/10 in control) or no ovary development (3/10 and 2/10 for control and treated workers, respectively). The frequency of occurrence of each ovary status category did not differ among treatments (odds ratio = 0.599, P = 1; odds ratio = 0.000, P = 1; and odds ratio = 3.611, P = 0.582 for developed, incipient and no developed ovaries, respectively).

Egg dumping behaviour was observed in 2/6 (control) and 1/5 (treated) of the colonies, which led to similar occurrence of this behaviour among experimental groups (odds ratio = 1.879, P = 1).

Methylation differences

Overall, the mean mapping efficiency was (mean ± SD) 39.5 ± 2.9%. The percentage of methylated CpG’s in control conditions was a mean of 0.5 ± 0.2% (mean ± SD). Similar levels were also found in non-CpG methylation contexts, a mean of 0.4 ± 0.1% of CHG’s were methylated and 0.4 ± 0.1% CHH’s (mean ± SD, referring ‘H’ to any base other than guanine). For the set of workers treated with decitabine, methylation rates were almost two-fold, with an average of 0.9 ± 0.1, 0.7 ± 0.1, and 0.8 ± 0.1 (mean ± SD), for CpG, CHG and CHH methylation types, respectively.

The current B. terrestris annotation file (Refseq accession no. GCF_000214255.1) was used as a reference to calculate the proportion of methylated reads for the C-context that were located in annotated loci. Variation in CpG methylation between treatment levels was much higher than variation caused by colonies (Fig. S1 A,B). Workers from treated colonies showed higher percentages of loci with significant levels of methylation (0.08% vs 0.04%, as group average, Table 1). Of the significantly methylated loci, methylation fraction (number of C reads / total number of reads) appears to be similar between treated and control groups. (Table 1).

Table 1 Overview of mapping efficiency, methylation rates in CpG, CHG and CHH contexts, number of significantly methylated loci and methylation fraction per loci for the 14 workers subjected to WGBS

A total of 22 loci (methylation differences > 10%, q = 0.05) showed significant methylation differences in the CpG context between treated and control groups. Of these, 20 differentially methylated loci between treated and control bees that could be successfully annotated. These loci were associated with GO enrichment for terms including regulation of oogenesis, oocyte mRNA localisation, neuronal function and various metabolic processes (Table 2). No differentially methylated loci were identified in a non-CpG context.

Table 2 GO enrichment analysis- biological processes (BP) terms of the set of differentially methylated loci for treated and untreated workers at 10% methylation difference threshold


Social insects are among the most successful on Earth, due to their ability to divide tasks among castes and to cooperate (see [32] and references therein). In many cases, individuals with different functions within the colony are genetically identical, which immediately raises the question whether and how epigenetic regulation contributes to this phenotypic differentiation. Previous studies on honeybees have shown that DNA methylation affects larval differentiation into worker-queen castes [23, 24]. In this study, we have experimentally manipulated the methylomes of founder queens of B. terrestris belonging to the same single parental line by using decitabine in order to infer the role of DNA methylation on phenotypic features such as caste determination and reproductive behaviour. The resulting colonies were monitored during ten weeks and their offspring was compared with colonies receiving a control treatment. Furthermore, the efficacy of the treatment was tested for a random subset of workers by analysing their methylomes using WGBS (see graphical summary in Supplementary Material).

Effects of decitabine on colony development

Chemically-induced DNA demethylation has been already proven to contribute to phenotypic variation in insects [33, 34]. In this study, we used full colonies (queen-right) instead of queen-less micro-colonies of bumblebee workers to assess the effect of altered methylation on colony development and caste determination, which is impossible to assess when using unfertilized workers. Our results showed that colonies that received a continuous supply of decitabine grew better than control colonies receiving no decitabine. In absence of genetic differences among queen founders, and considering that there were no temporal differences for the succession of main colony events, such an effect was mainly due to a higher egg-laying activity by the queen itself, combined with the higher maturation success of the queen brood. Gene ontology for the set of differentially methylated genes in the queen female offspring included oogenesis regulation and oocyte mRNA localization, consistent with the suggestion that founder queens got an enlarged progeny from their set of eggs. Moreover, decitabine is supposed to have no effect on post-mitotic adult specimens, such as the founder queen in this particular setup (7).

In our experiment, we did not analyse direct effects of decitabine on the queen methylomes, as this would have precluded the completion of the experiment itself, and therefore we cannot assess how founder queens responded to decitabine at a molecular level. However, workers were certainly exposed to the effect of the chemical, and therefore a more collaborative behaviour by the treated workers (less dominance, more nursing activity, etc.) might lead to the development of a bigger brood by the queen. This finding is in contrast with the higher aggressiveness, dominance and reproductive behaviour exhibited by workers when exposed to decitabine in absence of the founder queen (7). In queen-right colonies, this behavior would translate into early disruption of the cooperation phase and more abundant male production at the colony level [35]. In our experiment we found no signs of an earlier competition point in treated colonies, and male production tended to be higher in control colonies. Consistent with these findings, we did not find any indication of worker aggressiveness or altered behavior, and our dissections did also not indicate differences in worker ovary development among groups. These two papers show that altering methylation in different social contexts (Queenright vs queenless) has different effects. This intriguingly mirrors the predictions for genomic imprinting in different social contexts as predicted by Queller and Strassmann [36]. DNA Methylation is a major component of the imprinting systems in mammals and flowering plants [37].

The better colony development observed here also appears to be in contrast to the findings of Ellers et al. [38], who suggested that decitabine has an overall anti-metabolic activity that leads to a decline of the physical condition of the insects, and therefore results in phenotypic differences among treatments that are not dependent on DNA methylation. The concentration of decitabine that was used in our experiments was the same as the one used by Amarasinghe et al. (7), who showed that this concentration had no deleterious effects on bumblebee workers. Administration of decitabine, or the solvent itself (acetic acid), had a large effect on queen survival and egg laying when it was administered immediately after queen start-up. Due to its antimicrobial activity, low concentrations of acetic acid are commonly used in the food industry as preserving agent, suggesting that addition of acetic acid could have affected the symbiotic gut flora of the queens [39, 40]. Such effects were not observed when the sugar water treatments were offered two weeks later. We hypothesize that post-hibernating queens represent a sensitive stage in the colony cycle in terms of survival [41] and therefore the applied treatments, or the resulting changes in gut flora, may have affected colony initiation or even caused death. Moreover, our data have shown that the use of decitabine translates into differences in methylation patterns between the two experimental groups, including processes related to oogenesis, which also contradicts a cytotoxic role of the compound.

Use of decitabine did not induce early queen production, indicating that altered methylation patterns did not affect caste determination. Queen production in Bombus hypnorum can be induced by exposing the last larvae instar to juvenile hormone [42]. However, exposition to juvenile hormone had no proven effect in B. terrestris, as caste is known to be determined very early in larval development [43, 44]. Caste determination in bumblebees has been classically assumed to be controlled by the action of queen pheromones and the endocrine profile of the larvae itself [43, 45]. Further research is needed to elucidate if endocrine differences might be impacted by differential methylation in larvae, and the efficacy of such mechanism compared to the direct exposition of larvae to queen pheromones [46].

DNA methylation patterns in bumblebees

In invertebrates, levels of DNA methylation are much lower than in vertebrates [47, 48]. In insects, DNA methylation is concentrated in gene bodies and associated with more stable patterns of gene expression [48]. The average cytosine methylation (CpG) in the genome of honeybees is 0.7%, and much lower figures for the other methylation types (less than 0.2%) [49]. In Bombus terrestris, we found an average CpG methylation of 0.5% in the control group, and similar methylation rates for the other methylation types. These results largely correspond with recent findings in Bombus terrestris audax, where non-CpG methylation rates were 0.4 (CHG) and 0.5% (CHH) [15]. Compared to CpG methylation, however, CHG and CHH methylation were not particularly enriched in coding regions, which indicates that they have a minor regulatory role. Even though our experimental bees were treated to induce methylation differences, we found a very limited amount of loci with significant methylation levels for non-CpG types, which seems to corroborate that these types of methylation have a lower impact on gene expression.

Because we used founder queens from a pure genetic line, they do not differ (at least significantly) in their genetic information, making them excellent models to assess phenotypic effects of the methylation differences we experimentally induced. The wasp Nasonia has emerged as a model for DNA methylation studies in insects, due to its naturally inbred nature [30]. Despite the controlled variation at the genetic level, methylomes might still be subject to other sources of variation [50] and the stability and inheritance of methylomes has been subject to debate [51]. Despite this controversy, the latest evidence from honeybees suggests that even though gene body CpG methylation can oscillate during development, it is kept mostly invariable in the germline, likely to preserve function and methylation patterns over generations [49]. Experiments with Nasonia have also indicated stable inheritance of methylation status through generations [30]. This also agrees with recent findings that have shown high-inter-colony variation in methylation [15].

Efficacy of decitabine to alter methylation patterns

Although decitabine clearly affected the methylation patters on callow workers, it was not obvious that there was a reduction in overall cytosine methylation. While Amarasinghe et al. (7) attributed this finding to an artefact of the detection technique (methylation sensitive AFLP, known as MSAP, [52], our WGBS data also indicated higher methylation levels for the treated bees. Decitabine has been recently used in wasps to manipulate methylation levels, and the resulting methylomes were analysed by WGBS [34]. This study also showed that the chemical does not work as a pure demethylation agent (thus visible in the global % of cytosine methylation), as it has been previously described (reviewed in [53, 54]). Considering that the effect of the drug seems to be context-dependent [34], future studies that aim to manipulate phenotypes by using decitabine should carefully check the effect of the compound at the molecular level using the best resolution available.

The topical administration of another DNA methyltransferase inhibitor such as RG108 on honeybee workers induced reduction in global levels of DNA methylation, which was confirmed by an ELISA-based methylation assay. Such overall reduction resulted in a phenotypic consequence for these workers, that showed increased lifespan ([55]). A similar methodology was previously used by Biergans et al. [56], who administered RG108 or Zebularine to the honeybee thorax, proving a significant decrease in global methylation levels -here measured by capillary electrophoresis- for both agents.

We detected higher rates of cytosines that were methylated along the genome in treated bees, and correspondingly a larger list of loci that were significantly methylated within the treated set of workers. Although we have not checked the expression of these loci, recent reports have associated WGBS with expression patterns, and found a negative correlation between methylation status of a given loci and levels of expression [15, 49]. Therefore, it is likely that treated bees showed higher expression levels of the differentially methylated genes. Our GO enrichment analyses were suggestive of a non-random effect of decitabine in the identification of DML related to oogenesis regulation and oocyte mRNA localisation. Interestingly, increased expression of these pathways would be in agreement with the positive response that the compound produced in terms of brood development and production of workers.


There is increasing evidence that DNA methylation has a pronounced impact on the phenotype expressed by insects [57]. In the particular case of social insects, epigenetically-driven phenotypic differences may be particularly important, as they may affect the different roles that colony members with identical genetic information play within the colony organization. In our experiments, we have used an inbred line of bumblebee queen founders, and epigenetic differences were experimentally induced by adding the DNA methyl-transferase inhibitor decitabine to sugar provisions. Addition of the chemical resulted in altered DNA methylation patterns that led to a set of differentially methylated loci (with smaller methylation levels at the treated group) including some oogenesis. In contrast to previous research, queen production over worker production or early worker reproduction were not induced, and colonies showed a better cooperative behaviour, which led to higher worker production and less males. Considering that the genome of Bombus is rather simple and well annotated, more targeted work (targeted knock-downs) is needed to establish direct links between phenotypes (caste, reproduction) to causal CpG methylation on specific genes.


Experimental setup

A total of 24 lab-reared B. terrestris queens (Biobest Group, Westerlo, Belgium) from one single genetically pure, inbred line were used for all experiments. Males were obtained from separate colonies of the same original population, and they were discarded after one mating event. Fecundated queens that survived hibernation were set in separate cages and kept at 26 °C and 60% humidity on a diet of 50% v/v Biogluc® and Gamma irradiated honeybee collected, multifloral pollen ad libitum until they developed into full colonies.

A stock solution of 5-aza-20-deoxycytidine (decitabine) was made by dissolving 5 mg of decitabine (Sigma Aldrich, Belgium) in 2 ml of 1: 1 v/v acetic acid: distilled water solution. 10 mM decitabine was added to sugar water (0.0925% v/v) (7) and fed to two test groups of six colonies each at either one or three weeks after queen start-up (12 colonies in total). The control group (12 colonies) was fed with standard sugar water, plus the solvent (acetic acid), also at 0.0925% v/v. Solutions were provided to each colony on a weekly basis throughout the entire experiment. Neither the treatment nor the control colonies that received the sugar treatment from the first week after the start of the experiment did develop properly. Just 1 out of 6 queens receiving the treatment managed to lay eggs and start a colony, while this number was 0 out of 6 in the control group. As a result, these colonies were removed from subsequent analyses. Because similar effects were found for the control and treatment colonies, we conclude that the addition of the solvent to the sugar solution had a deleterious effect on the survival/fitness of queens shortly after hibernation. No such effects were observed when the feeding solutions were provided to the colonies two weeks later and all colonies followed the subsequent development stages. Exposure 2 weeks after queen start-up, as done in our experiment, still guaranteed that all resulting larvae were exposed to decitabine.

Colony development and behavioral records

Colony development was monitored over a 10-week period. For each colony, we assessed the developmental time by determining the timing of first egg-laying, first pupation and first emergence of adults (workers, males). For each colony, we also counted the number of adult workers, pupae, larvae, dead larvae, and egg cups at week 8 and week 10. Counts at week 8 and 10 represent the colony development before and after the competition point, respectively. Counts of dead larvae at week 10 was not possible because it is nearly impossible to obtain accurate estimates in fully occupied nest boxes. Colony development was categorized into no development (no brood being produced), and successfully developed colonies, that follow the normal timing of development for the subsequent phases [19]. In addition, the behavior of the colony was recorded by observing colonies twice a week for 15 min throughout the trial. During these observations, we annotated queen dominance, worker aggressiveness as well as indications of worker reproduction (competition point), egg dumping, and nursing behavior. Egg dumping was expressed through bees placing abandoned eggs on top of the sugar water reservoir.

DNA methylation on workers

After eight weeks, two colonies were randomly chosen per treatment. From each selected colony, up to five coetaneous adult workers were collected. Each worker was dissected using a fresh Ringers solution (Sodium chloride 2.25 g/L, Potassium chloride 0.10 g/L, Calcium chloride 6H2O 0.12 g/L, Sodium bicarbonate 0.05 g/L) to extract the brain tissue. At the same time, ovary development was assessed following the scale provided by [35], which we simplified to a three-level score of no development, intermediate and fully developed ovaries. DNA was then extracted from each flash-frozen brain sample using the EZNA Insect kit (Omega Bio-Tek), following the instructions from the manufacturer. DNA quality and quantity were determined by Nanodrop and Qubit® fluorospectrometers. Samples that yielded less than 1 μg DNA/μl were discarded from subsequent analyses.

DNA extraction and library

Bisulphite conversion of genomic DNA combined with next-generation sequencing (WGBS) was used to measure the methylation state of the whole genome, the methylome, at single-base resolution [58]. Fourteen WGBS non-directional libraries (seven workers per treatment, taken from two colonies per treatment, with three to four workers per individual colony) were prepared by BGI Tech Solutions Co (Hong Kong). This involves DNA fragmentation, adapter ligation, bisulphite treatment, size selection and amplification. Resulting libraries (14, one per bee brain sample) were sequenced by using 100 bp paired-end bisulfite sequencing on a HiSeq 2000 machine (Illumina, Inc.) by BGI Tech Solutions Co.

Data analyses

Colony development

To test whether colony development differed between founders that were assigned to different treatments, we used a Wilcoxon test with the number of days before egg laying, pupae appearance, male and worker emergence as dependent variables, while treatment was considered as fixed factor.

A generalized linear model with Poisson distribution was used to investigate whether brood size (the sum of egg cups and larvae) and the number of workers at eight and ten weeks after queen start-up differed between treatments. Treatment and counting week, as well as their interaction, were treated as fixed factors and brood size and the number of workers as dependent variables. Post-hoc analyses (Tukey HSD) were conducted to see whether the dependent variables differed significantly between treatments for each counting week. A similar model was used to investigate whether the number of dead larvae differed between treatments at week 8. Here treatment was the only fixed factor.

We calculated contingency tables to investigate whether the number of individuals assigned to the respective categories of ovary development differed between treatments. Significance was estimated using a Fisher’s exact test. Similar analyses were conducted to see whether dumping behaviour differed between treatments.

DNA methylation analyses

Poor quality reads/bases and adapter sequences were removed using the function of BBMap ( Bismark v0.18.1 [59] was used to align reads to the B. terrestris genome (GCF_000214255.1 [60]) using the non-directional protocol,remove PCR duplicates, and filter non-converted reads. The Bismark BAM file output was processed using methylKit [61] to remove base calls with a quality below Q30 and filtered to remove low (< 10 reads) and high (> 99.9th percentile) coverage loci. This processed file was loaded in methylKit and filtered to exclude loci that were present in less than four samples per group. A binomial test was used to make per-loci methylation status calls, using a 1% error rate. Only loci with significant levels of methylation in at least one sample were subsequently tested for differential methylation between groups using the Chi-squared test in methylKit, controlling for colony as a covariate and correcting for overdispersion. A minimum methylation difference of 10% was necessary for a locus to be considered differentially methylated, using an FDR-adjusted q-value of 0.05. Differentially methylated loci (DML) were further curated manually and loci where the result was driven by a single methylated sample excluded.

Differentially methylated loci (DML) were extracted from the genome and annotated using a custom-made database [62]. GO enrichment was conducted against all RNA features in the bumblebee genome using GOStats [63] to conduct a hypergeometric test, with significant GO terms identified using Benjamini-Hochberg correction (adjusted p-value 0.05).

Availability of data and materials

Sequence Read Archive (SRA) data analysed during the current study are available in the NCBI repository,, in read-only format. Other data and scripts generated or analysed during this study are included in the submission as supplementary information files.



DNA sequence context characterized by cytosines followed by two consecutive H residues, where H correspond to Adenine, Thymine or Cytosine


DNA sequence context characterized by cytosines followed by H residues (which correspond to Adenine, Thymine or Cytosine) and Guanine


DNA sequence context characterized by Cytosines followed by Guanine residues


Differentially Methylated Loci


DeoxyriboNucleic Acid


DNA (cytosine-5)- methyltransferase 1


Gene Ontology




Messenger RiboNucleic Acid


Whole Genome Bisulfite Sequencing


  1. 1.

    Wong AH, Gottesman II, Petronis A. Phenotypic differences in genetically identical organisms: the epigenetic perspective. Hum Mol Genet. 2005;14:R11–8.

    CAS  PubMed  Article  Google Scholar 

  2. 2.

    Hirsch S, Baumberger R, Grossniklaus U. Epigenetic variation, inheritance, and selection in plant populations. Cold Spring Harb Symp Quant Biol. 2012;77:97–104.

    CAS  PubMed  Article  Google Scholar 

  3. 3.

    Herrera CM, Bazaga P. Epigenetic correlates of plant phenotypic plasticity: DNA methylation differs between prickly and nonprickly leaves in heterophyllous Ilex aquifolium (Aquifoliaceae) trees. Bot J Linn Soc. 2013;171:441–52.

    Article  Google Scholar 

  4. 4.

    Herrera CM, Pozo MI, Bazaga P. Jack of all nectars, master of most: DNA methylation and the epigenetic basis of niche width in a flower-living yeast. Mol Ecol. 2012;21:2602–16.

    CAS  PubMed  Article  Google Scholar 

  5. 5.

    Pegoraro M, Marshall H, Lonsdale ZN, Mallon EB. Do social insects support Haig’s kin theory for the evolution of genomic imprinting? Epigenetics. 2017;12:725–42.

    PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    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.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  7. 7.

    Amarasinghe HE, Clayton CI, Mallon EB. Methylation and worker reproduction in the bumble-bee (Bombus terrestris). Proc R Soc B Biol Sci. 2014;281:20132502.

    Article  Google Scholar 

  8. 8.

    Tost J. DNA methylation: an introduction to the biology and the disease-associated changes of a promising biomarker. Mol Biotechnol. 2010;44:71–81.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  9. 9.

    Suzuki MM, Kerr AR, De Sousa D, Bird A. CpG methylation is targeted to transcription units in an invertebrate genome. Genome Res. 2007;17:625–31.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  10. 10.

    Zemach A, McDaniel IE, Silva P, Zilberman D. Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science. 2010;328:916–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  11. 11.

    Wang X, Wheeler D, Avery A, Rago A, Choi J-H, Colbourne JK, et al. Function and evolution of DNA methylation in Nasonia vitripennis. PLoS Genet. 2013;9:e1003872.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  12. 12.

    Feng S, Cokus SJ, Zhang X, Chen P-Y, Bostick M, Goll MG, et al. Conservation and divergence of methylation patterning in plants and animals. Proc Natl Acad Sci. 2010;107:8689–94.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  13. 13.

    Bonasio R, Li Q, Lian J, Mutti NS, Jin L, Zhao H, et al. Genome-wide and caste-specific DNA methylomes of the ants Camponotus floridanus and Harpegnathos saltator. Curr Biol. 2012;22:1755–64.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Glastad KM, Hunt BG, Goodisman MA. Evolutionary insights into DNA methylation in insects. Curr Opin Insect Sci. 2014;1:25–30.

    PubMed  Article  PubMed Central  Google Scholar 

  15. 15.

    Marshall H, Lonsdale ZN, Mallon EB. Methylation and gene expression differences between reproductive castes of bumblebee workers. Evol Lett. 2019;3:485–99.

    PubMed  PubMed Central  Article  Google Scholar 

  16. 16.

    Bourke AF. Dominance orders, worker reproduction, and queen-worker conflict in the slave-making ant Harpagoxenus sublaevis. Behav Ecol Sociobiol. 1988;23:323–33.

    Article  Google Scholar 

  17. 17.

    Hamilton WD. Genetical evolution of social behavior I. J Theor Biol. 1964; 7:1:16.

  18. 18.

    Hamilton WD. The genetical evolution of social behaviour. II J Theor Biol. 1964;7:17–52.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  19. 19.

    Duchateau MJ, Velthuis HHW. Development and reproductive strategies in Bombus terrestris colonies. Behaviour. 1988;107:186–207.

    Article  Google Scholar 

  20. 20.

    Maleszka R. Epigenetic integration of environmental and genomic signals in honey bees: the critical interplay of nutritional, brain and reproductive networks. Epigenetics. 2008;3:188–92.

    PubMed  Article  Google Scholar 

  21. 21.

    Lopez-Vaamonde C, Koning JW, Jordan WC, Bourke AF. No evidence that reproductive bumblebee workers reduce the production of new queens. Anim Behav. 2003;66:577–84.

    Article  Google Scholar 

  22. 22.

    Schmid-Hempel R, Schmid-Hempel P. Female mating frequencies in Bombus spp. from Central Europe. Insect Soc. 2000;47:36–41.

    Article  Google Scholar 

  23. 23.

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

    CAS  PubMed  Article  Google Scholar 

  24. 24.

    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.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  25. 25.

    Glastad KM, Gokhale K, Liebig J, Goodisman MA. The caste-and sex-specific DNA methylome of the termite Zootermopsis nevadensis. Sci Rep. 2016;6:1–14.

    Article  CAS  Google Scholar 

  26. 26.

    Alvarado S, Rajakumar R, Abouheif E, Szyf M. Epigenetic variation in the Egfr gene generates quantitative variation in a complex trait in ants. Nat Commun. 2015;6:1–9.

    Article  CAS  Google Scholar 

  27. 27.

    Standage DS, Berens AJ, Glastad KM, Severin AJ, Brendel VP, Toth AL. Genome, transcriptome and methylome sequencing of a primitively eusocial wasp reveal a greatly reduced DNA methylation system in a social insect. Mol Ecol. 2016;25:1769–84.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  28. 28.

    Burggren W. Epigenetic inheritance and its role in evolutionary biology: re-evaluation and new perspectives. Biology. 2016;5:24.

    PubMed Central  Article  Google Scholar 

  29. 29.

    Remnant EJ, Ashe A, Young PE, Buchmann G, Beekman M, Allsopp MH, et al. Parent-of-origin effects on genome-wide DNA methylation in the cape honey bee (Apis mellifera capensis) may be confounded by allele-specific methylation. BMC Genomics. 2016;17:226.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  30. 30.

    Wang X, Werren JH, Clark AG. Allele-specific transcriptome and methylome analysis reveals stable inheritance and cis-regulation of DNA methylation in Nasonia. PLoS Biol. 2016;14:e1002500.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  31. 31.

    Maleszka R. Epigenetic code and insect behavioural plasticity. Curr Opin Insect Sci. 2016;15:45–52.

    PubMed  Article  Google Scholar 

  32. 32.

    Gordon DM. From division of labor to the collective behavior of social insects. Behav Ecol Sociobiol. 2016;70:1101–8.

    PubMed  Article  Google Scholar 

  33. 33.

    Cook N, Pannebakker BA, Tauber E, Shuker DM. DNA methylation and sex allocation in the parasitoid wasp Nasonia vitripennis. Am Nat. 2015;186:513–8.

    PubMed  Article  Google Scholar 

  34. 34.

    Cook N, Parker DJ, Turner F, Tauber E, Pannebakker BA, Shuker DM. Genome-wide disruption of DNA methylation by 5-aza-2’deoxycytidine in the parasitoid wasp Nasonia vitripennis. BioRxiv. 2019;437202.

  35. 35.

    Duchateau MJ, Velthuis HHW. Ovarian development and egg laying in workers of Bombus terrestris. Entomol Exp Appl. 1989;51:199–213.

    Article  Google Scholar 

  36. 36.

    Queller DC, Strassmann JE. The many selves of social insects. Science. 2002;296:311–3.

    CAS  PubMed  Article  Google Scholar 

  37. 37.

    Zilberman D. The evolving functions of DNA methylation. Curr Opin Plant Biol. 2008;11:554–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  38. 38.

    Ellers J, Visser M, Mariën J, Kraaijeveld K, Lammers M. The importance of validating the demethylating effect of 5-aza-dC in model species. Am Nat. 2019;194:422–31.

    PubMed  Article  Google Scholar 

  39. 39.

    Billiet A, Meeus I, Van Nieuwerburgh F, Deforce D, Wäckers F, Smagghe G. Impact of sugar syrup and pollen diet on the bacterial diversity in the gut of indoor-reared bumblebees (Bombus terrestris). Apidologie. 2016;47:548–60.

    CAS  Article  Google Scholar 

  40. 40.

    Pozo MI, Bartlewicz J, van Oystaeyen A, Benavente A, van Kemenade G, Wäckers F, et al. Surviving in the absence of flowers: do nectar yeasts rely on overwintering bumblebee queens to complete their annual life cycle? FEMS Microbiol Ecol. 2018;94: fiy196.

  41. 41.

    Beekman M, van Stratum P, Lingeman R. Diapause survival and post-diapause performance in bumblebee queens (Bombus terrestris). Entomol Exp Appl. 1998;89:207–14.

    Article  Google Scholar 

  42. 42.

    Röseler PF, Röseler I. Morphological and physiological differentiation of the caste in the bumblebee species Bombus hypnorum (L.) and Bombus terrestris (L.). Zool Jahrb Physiol Bd. 1974;78:175–98.

    Google Scholar 

  43. 43.

    Röseler P-F. Unterschiede in der Kastendetermination zwischen den Hummelarten Bombus hypnorum und Bombus terrestris/differences in the caste determination between the bumblebee species Bombus hypnorum and Bombus terrestris. Z Für Naturforschung B. 1970;25:543–8.

    Article  Google Scholar 

  44. 44.

    Cnaani J, Robinson GE, Hefetz A. The critical period for caste determination in Bombus terrestris and its juvenile hormone correlates. J Comp Physiol A. 2000;186:1089–94.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  45. 45.

    Cnaani J, Schmid-Hempel R, Schmidt JO. Colony development, larval development and worker reproduction in Bombus impatiens Cresson. Insect Soc. 2002;49:164–70.

    Article  Google Scholar 

  46. 46.

    Bortolotti L, Duchateau MJ, Sbrenna G. Effect of juvenile hormone on caste determination and colony processes in the bumblebee Bombus terrestris. Entomol Exp Appl. 2001;101:143–58.

    CAS  Article  Google Scholar 

  47. 47.

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

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  48. 48.

    Bewick AJ, Vogel KJ, Moore AJ, Schmitz RJ. Evolution of DNA methylation across insects. Mol Biol Evol. 2016;34:654–65.

    PubMed Central  Google Scholar 

  49. 49.

    Harris KD, Lloyd JP, Domb K, Zilberman D, Zemach A. DNA methylation is maintained with high fidelity in the honey bee germline and exhibits global non-functional fluctuations during somatic development. Epigenetics Chromatin. 2019;12:62.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  50. 50.

    Seymour DK, Becker C. The causes and consequences of DNA methylome variation in plants. Curr Opin Plant Biol. 2017;36:56–63.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  51. 51.

    Bonduriansky R. Day T. Extended Heredity: A New Understanding of Inheritance and Evolution. Princeton University Press; 2020.

    Google Scholar 

  52. 52.

    Schulz B, Eckstein RL, Durka W. Scoring and analysis of methylation-sensitive amplification polymorphisms for epigenetic population studies. Mol Ecol Resour. 2013;13:642–53.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  53. 53.

    De Vos D, van Overveld W. Decitabine: a historical review of the development of an epigenetic drug. Ann Hematol. 2005;84:3.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  54. 54.

    Linnekamp JF, Butter R, Spijker R, Medema JP, van Laarhoven HWM. Clinical and biological effects of demethylating agents on solid tumours–a systematic review. Cancer Treat Rev. 2017;54:10–23.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  55. 55.

    Cardoso-Júnior CA, 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. 2018;92:21–9.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  56. 56.

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

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  57. 57.

    Corona M, Libbrecht R, Wheeler DE. Molecular mechanisms of phenotypic plasticity in social insects. Curr Opin Insect Sci. 2016;13:55–60.

    PubMed  Article  PubMed Central  Google Scholar 

  58. 58.

    Krueger F, Kreck B, Franke A, Andrews SR. DNA methylome analysis using short bisulfite sequencing data. Nat Methods. 2012;9:145.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  59. 59.

    Krueger F. Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications bioinformatics. 2011;27:1571–2.

    CAS  PubMed  PubMed Central  Google Scholar 

  60. 60.

    Sadd BM, Barribeau SM, Bloch G, De Graaf DC, Dearden P, Elsik CG, et al. The genomes of two key bumblebee species with primitive eusocial organization. Genome Biol. 2015;16:76.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  61. 61.

    Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, et al. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13: R87.

  62. 62.

    Bebane PS, Hunt BJ, Pegoraro M, Jones AC, Marshall H, Rosato E, et al. The effects of the neonicotinoid imidacloprid on gene expression and DNA methylation in the buff-tailed bumblebee Bombus terrestris. Proc R Soc B. 2019;286:20190718.

    CAS  PubMed  Article  Google Scholar 

  63. 63.

    Falcon S, Gentleman R. Using GOstats to test gene lists for GO term association. Bioinformatics. 2006;23:257–8.

    PubMed  Article  CAS  Google Scholar 

Download references


We thank Stef Rutten from Biobest Group for his timely help during the experiment, and his long and careful work on the bumblebee genetic lines.

This research used the SPECTRE High Performance Computing Facility at the University of Leicester.


We are grateful to the Fonds Wetenschappelijk Onderzoek (FWO) for providing funding for this research via application 12A0716N. This funding agency has no role in the design of the study, the interpretation of the data, and the preparation of the manuscript.

Author information




MIP, HJ, EM and JMGZ conceived the study. MIP and HJ wrote the initial manuscript and all coauthors contributed to the last version. MIP and BH analyzed the data. GVK and MIP performed the experiments. FW provided the experimental subjects. All authors have read and have approved the submitted version of this manuscript.

Corresponding author

Correspondence to María I. Pozo.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1: Figure S1. (A)

Principal component Analyses (PCA) plot generated using methylKit function PCASamples, showing CpG methylation for the 7 control bee samples and the seven treated bee samples. (B). Dendogram generated using methylKit function clusterSamples showing sample cluster by treatment. Red labels indicate control samples and blue labels represent decitabine-treated samples. Black labels in dendogram depict the colony of origin of each sample.

Additional file 2.

Table 1 Overview of mapping efficiency, methylation rates in CpG, CHG and CHH contexts, number of significantly methylated loci and methylation fraction per loci for the 14 workers subjected to WGBS.

Additional file 3.

Text files with data and code used for analyses. We provide datasets on colony developmental data (countings: ctres, temporal overview: events3c) and ovary development (ovarydev). Additional text files are available for the R code (R script) and WGBS analyses (supp_code).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Pozo, M.I., Hunt, B.J., Van Kemenade, G. et al. The effect of DNA methylation on bumblebee colony development. BMC Genomics 22, 73 (2021).

Download citation


  • DNA methylation
  • Social insects
  • Epigenetics
  • Colony development