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 volume 17, Article number: 226 (2016)
Intersexual genomic conflict sometimes leads to unequal expression of paternal and maternal alleles in offspring, resulting in parent-of-origin effects. In honey bees reciprocal crosses can show strong parent-of-origin effects, supporting theoretical predictions that genomic imprinting occurs in this species. Mechanisms behind imprinting in honey bees are unclear but differential DNA methylation in eggs and sperm suggests that DNA methylation could be involved. Nonetheless, because DNA methylation is multifunctional, it is difficult to separate imprinting from other roles of methylation. Here we use a novel approach to investigate parent-of-origin DNA methylation in honey bees. In the subspecies Apis mellifera capensis, reproduction of females occurs either sexually by fertilization of eggs with sperm, or via thelytokous parthenogenesis, producing female embryos derived from two maternal genomes.
We compared genome-wide methylation patterns of sexually-produced, diploid embryos laid by a queen, with parthenogenetically-produced diploid embryos laid by her daughters. Thelytokous embryos inheriting two maternal genomes had fewer hypermethylated genes compared to fertilized embryos, supporting the prediction that fertilized embryos have increased methylation due to inheritance of a paternal genome. However, bisulfite PCR and sequencing of a differentially methylated gene, Stan (GB18207) showed strong allele-specific methylation that was maintained in both fertilized and thelytokous embryos. For this gene, methylation was associated with haplotype, not parent of origin.
The results of our study are consistent with predictions from the kin theory of genomic imprinting. However, our demonstration of allele-specific methylation based on sequence shows that genome-wide differential methylation studies can potentially confound imprinting and allele-specific methylation. It further suggests that methylation patterns are heritable or that specific sequence motifs are targets for methylation in some genes.
Sexual reproduction is generally characterized by the equal transmission of nuclear genomes from mothers and fathers to offspring, and by equal expression of maternal and paternal alleles in offspring [1, 2]. However violations of Mendelian inheritance can arise when the genetic interests of the mother and father diverge [2, 3]. Males benefit if their offspring can secure more resources from the offspring’s mother, whereas females normally benefit if their offspring are provisioned equally regardless of paternity. These sexual conflicts can lead to parent-of-origin effects, arising from genomic imprinting, where offspring phenotype is altered depending on whether an allele is inherited from a mother or a father . Imprints are laid down during gametogenesis, during an epigenetic reprogramming event that results in the non-equivalence of paternal and maternal genomes . For example, mice embryos derived from two maternal genomes or two paternal genomes via pronuclear transfers are non-viable due to parent-specific imprinting [6, 7].
Parent-specific epigenetic modifications that give rise to genomic imprinting often occur via DNA methylation, which involves the addition of a methyl group to cytosine residues of DNA . In addition to its role in imprinting, DNA methylation varies between cell types due to cellular differentiation , and may also vary according to genotype due to allele-specific methylation . DNA methylation can alter the expression and regulation of genes . When methylation occurs in promoter regions, it results in suppression of gene expression. When methylation occurs in the exons and introns of gene bodies, it increases or maintains gene expression and may mediate alternate splicing [12–14]. Gene body methylation is conserved among plants, invertebrates and mammals, and may therefore represent an ancestral function of DNA methylation [15, 16]. Not all invertebrates have DNA methylation, but if it is present it is generally restricted to gene bodies and affects only a small percentage of cytosines genome-wide .
Colonies of Hymenopteran insects (bees, wasps and ants) provide a situation that is ripe with potential for intra-genomic conflict and thus for imprinting [18–20]. Hymenoptera are haplo-diploid. Females arise from fertilized eggs and are diploid, while males derive from unfertilized eggs by arrhenotokous parthenogenesis and are haploid. Workers almost never lay eggs when a queen is present but are physiologically able to produce haploid eggs that give rise to males. This ability allows workers of some species to occasionally ‘cheat’ and lay haploid eggs that produce viable males [21, 22]. In many eusocial hymenopterans the queens are polyandrous . Because females mate with multiple males, any male who can epigenetically modify his female offspring so that they are more likely to become reproductively active will be favored by selection, as his offspring will outcompete those of males that do not do so [18, 20].
Unlike mammals where parthenogenesis generates inviable offspring due to imprinting [24–26], some Hymenopterans can produce diploid female offspring from unfertilized eggs by thelytokous parthenogenesis . Thelytoky, the asexual production of female offspring, is rare in honey bees generally but ubiquitous in workers of the subspecies Apis mellifera capensis (hereafter Capensis) [28, 29]. In Capensis, mated queens reproduce sexually, identically to all other honey bee subspecies, by the union of egg and sperm nuclei within a newly-laid egg. In contrast to other honey bee species, Capensis workers lay diploid eggs that develop into females ([30, 31]; Fig. 1 ), by the fusion of two maternally-derived pronuclei, as if one pronucleus acted as a sperm, again within a newly-laid egg . Embryogenesis otherwise occurs normally as it would in a fertilized embryo, resulting in the production of diploid females.
Thelytokous parthenogenesis as described above generates embryos that are ‘pseudoclones’ of their worker parent. That is, thelytokous offspring receive one copy of the two grandparental genomes and are largely identical to their thelytokous parent, unless recombination during meiosis has resulted in loss of heterozygosity . While thelytoky enables workers to genetically re-incarnate themselves, it is currently unknown how the epigenetic component of the genome is inherited following thelytokous meiosis. Sexually-produced embryos laid by a queen inherit both a maternally- and a paternally-derived set of chromosomes, whereas thelytokous eggs contain pronuclei that have undergone gametogenesis within the ovaries of female workers. Thus thelytokous embryos begin life with two maternally-inherited pronuclei, (Fig. 1) and may show different patterns of methylation to sexually-produced embryos. Therefore Capensis provides the ideal system in which to examine parent-of-origin effects, as it naturally generates diploid bi-parental and diploid uni-parental embryos.
As yet we have only indirect evidence for genomic imprinting and parent-of-origin effect in social insects in general, and honey bees in particular. However three observations suggest the presence of imprinting in the honey bee. First, methylation patterns differ significantly between unfertilized eggs and sperm, suggesting that queens and drones differentially methylate their gametes . Second, reciprocal crosses between honey bee subspecies show substantial parent-of-origin effects for reproductive traits, again suggestive of paternal imprinting [34–36]. Finally, in reciprocal crosses of Africanised and European honey bees, 46 transcripts showed significant parent-of-origin expression effects . A second study using Africanised and European honey reciprocal crosses determined that hybrid workers showed paternally biased phenotypes and expression patterns, particularly relating to reproduction, suggesting that in honey bees worker reproduction is driven by patrigenes .
In addition to its potential role in genomic imprinting, DNA methylation has additional functions in Hymenoptera [39–41]. In honey bees it is involved in larval development , and differentiation into queen and worker castes . Methylation patterns differ between the brains of honey bee queens and workers , as well as between adult workers during behavioural maturation . The majority of genomic methylation is associated with constitutively-expressed, highly-conserved genes [12, 46–49] and, approximately half of all genes are methylated . Methylation may provide a mechanism by which gene dosage is regulated between the haploid and diploid castes . How methylation is regulated, and the mechanisms by which DNA methylation is directed to particular sites and genes is unknown.
Here we examine the methylation patterns of fertilized, sexually produced, bi-parental embryos laid by Capensis queens with asexually produced, parthenogenetic, uni-parental embryos produced by their thelytokous daughter workers (Fig. 1). We determine whether thelytokous embryos with two maternal genomes have different patterns of methylation to sexually-produced embryos with a maternal and paternal genome. To control for genotype, experimental queens were artificially inseminated, each with sperm from a single drone. Therefore, the worker progeny were genetically homogenous, each worker carrying one allele from the inseminating male, and one of the queen’s two alleles. At a cohort level, the workers carried just three alleles, two that were maternally derived, and one that was paternally derived. On average, the worker-laid embryos were genetically identical to those of the queen-derived embryos, and differed only in epigenetic modifications. The presence of three alleles in common through two generations enabled us to follow the methylation patterns of each allele in both fertilized and thelytokous embryos.
If paternal imprinting exists, we would predict significantly more methylated sites in sexually-produced embryos due to methylation present in the paternally-inherited sperm pronucleus. In thelytokous embryos, we expect fewer methylated sites overall. Conversely, for maternally imprinted alleles we would predict significantly higher levels of methylation at specific maternally imprinted sites in thelytokous embryos, due to a double dosage of two maternally imprinted genomes. An absence of differential methylation between thelytokous and fertilized diploid embryos would be strong evidence that imprinting does not occur in the honey bee. Evidence of allele-specific methylation would show that epigenetic modifications are heritable in the honey bee or that they are invariably established according to cis-mediated genotype, but could also confound our ability to detect imprinting.
Sample collection, DNA extraction and sequencing
We produced three colonies (1–3) by insemination of three virgin Capensis queens, each with the semen of a single Capensis drone. The parents were sourced in Stellenbosch, South Africa. The inseminated queens were introduced into separate queenless host colonies. After the introduced queens began laying we collected between 300 and 1200 fertilized embryos from each colony (0–70 h old, hereafter ‘fertilized embryos’) into 99 % ethanol.
Once all host workers had been replaced by worker-offspring of the inseminated queens (about two months), we removed the queens from the colonies. The now queenless Capensis workers commenced laying thelytokous eggs approximately 1 week after the removal of the queen. These Capensis daughter workers produced diploid eggs via thelytokous parthenogenesis. We collected between 300 and 900 worker-laid embryos from each colony (0–70 h old, hereafter ‘thelytokous embryos’), again into ethanol.
We used microsatellite analysis to confirm that individual embryos produced by workers from this colony were all thelytokous diploid females and not arrhenotokous haploid males . We also confirmed that DNA from the fertilized and thelytokous embryos contained the same microsatellite alleles at five unlinked loci, and that there were no more than three alleles at each microsatellite confirming two alleles were derived from the queen and one allele from the sperm from the single drone . This proves that our colonies had not been invaded by reproductive parasites, as often happens in queenless A. m. capensis colonies [51, 52].
Our design resulted in three independent sets of paired embryos (‘fertilised’ and ‘thelytokous’) that were diploid and genetically (but not necessarily epigenetically) equivalent (Fig. 1a). Embryos were collected over similar time frames, so if methylation patterns change during embryonic development, these changes would be equivalent in both cohorts.
Colony 1 generated sufficient fertilized (1280) and thelytokous (870) embryos for whole genome bisulfite sequencing. To extract DNA we removed the ethanol, rinsed the embryos in water, and resuspended them in G2 buffer (Qiagen). Embryos were lysed in a tissue lyser and treated with 2 mg/ml RNAse at 37 °C for 30 min, Proteinase K at 56 °C for two hours, and DNA purified using 20G genomic tips (Qiagen). We quantified the amount of DNA in the extracts using the Qubit Flurometer (Life Technologies). The final yield for fertilized and thelytokous embryo samples was 8.3 and 2.5 μg, respectively. We added lambda bacteriophage DNA (0.1 % (w/w), Promega) to each sample as a spike-in control to evaluate the efficacy of the bisulfite conversion . Bisulfite conversion, library construction and sequencing were performed at the Beijing Genomics Institute as described in Drewell et al. .
Embryos from colonies 2 and 3 were reserved for bisulfite PCR validation of candidate genes.
Whole genome bisulfite sequence mapping
We used BSMAP v2.74 to determine methylation state of cytosine residues from the whole genome sequencing . Bisulfite reads were aligned to the A. mellifera genome assembly 2.0  as described in Drewell et al. .
A. m. capensis is a different subspecies to the reference genome - A. m. ligustica [55, 56]. If there is significant polymorphism between the Capensis and reference genomes then the proportion of bisulfite reads that could be successfully aligned to the reference genome would be reduced. To assess the extent of variation between the Capensis and reference genomes, we obtained A. m. capensis genome sequence data from Wallberg et al,  and generated a ‘Capensis2.0’ reference genome. Briefly, we aligned the A. m capensis reads  to the A. mellifera genome assembly 2.0 using BWA v0.5.9 . The alignment of the 75 bp Capensis reads resulted in 16.8-fold coverage of the reference genome. We utilized the alignment to determine SNPs and attempted to generate an alternate ‘Capensis2.0’ reference. File conversions, SNP calling and creation of an alternate reference was performed with SAMtools v0.1.16  and GATK v3.3 . We repeated BSMAP bisulfite mapping to the generated Capensis2.0 reference. Compared to the original mapping results, fold coverage and methylation calls (CG/CHG/CHH) were not significantly improved (Additional file 1: Table S1). We therefore continued analysis based on the alignment to the A. mellifera Assembly 2.0, so that our data can be directly compared with existing methylomes [33, 42, 44].
We determined the bisulfite conversion rate in our DNA samples by calculating the C-to-T conversion rate in the lambda spike-in control. Using BSMAP we aligned reads to the bacterophage lambda genome (Genbank J02459.1). Based on the number of converted cytosines to total lambda cytosine coverage, we estimated the bisulfite conversion rate as 0.9970 for the fertilized embryo sample, and 0.9965 for the thelytokous embryo sample.
The number of methylated cytosines in each sample was calculated as in Drewell et al. . At each cytosine site, we compared the number of bisulfite-converted reads (unmethylated) and the number of non-converted reads (methylated). We performed a binomial test to determine significant methylation, with a probability of one minus the lambda phage DNA conversion rate, correcting for multiple comparisons .
To quantitate differential methylation between the two samples, the output from BSMAP was analysed using the R package methylKit . To compare the frequency of methylated cytosines at CpG sites, we first excluded sites with less than 5 fold coverage. Sites containing a methylation difference of 25 % or greater, with a P value of 0.05 or less were deemed to be differentially methylated, and annotated against the official gene set from A. mellifera assembly 2.0.
Bisulfite PCR validation
Bisulfite PCR was conducted on samples of fertilized and thelytokous eggs from the two remaining colonies, 2 and 3. DNA was extracted from 300 embryos per sample, using the DNeasy blood and tissue kit (Qiagen). Two independent bisulfite conversion reactions were performed per sample (Zymo EZ DNA Methylation-Direct kit). Bisulfite PCR assays were designed for five candidate genes, Stan (GB18207), Stoned B (GB17165), Sap30 (GB18386), Syd (GB15356) and Pcl (GB16657; Table 3). Nested primers were selected spanning a region containing as many differentially methylated sites as possible within each candidate gene and PCR products were amplified (Additional file 1: Table S2). Previous studies have noted the difficulties of amplifying regions greater than 500–600 bp in length from bisulfite converted DNA [63, 64]. Using the two-step nested PCR and KAPA 2G Robust DNA polymerase (KAPA biosystems), we were able to amplify 793 bp (Stan), 650 bp (Stoned B), 1022 bp (Sap30), 1224 bp (Syd) and 1040 and 1375 bp (Pcl), fragments from the candidate genes (Additional file 1: Table S2). For three genes (Stan, Stoned B and Sap30), PCR products from each independent bisulfite reaction, for each sample were separately cloned into the TOPO vector (Invitrogen) and a minimum of 20 individual amplicons were sequenced for each sample (Macrogen, Seoul). Fragments amplified from the other two genes (Syd and Pcl) could not be cloned or did not generate sufficient clones for analysis after repeated attempts.
CpG methylation frequency in the Capensis genome
On average we obtained 8.6-fold and 10.1-fold genome coverage for fertilized and thelytokous embryo samples, respectively, achieving a similar level of coverage to previous honey bee methylation studies [33, 44, 45]. Genome-wide, the proportion of total CG sites covered by 2 to 30 reads was 79.4 % for fertilized embryos, and 80.3 % for thelytokous embryos (Table 1). There were 114,156 methylated cytosines in the CG context in fertilized embryos and 99,923 in thelytokous embryos. The proportion of CG sites in the genome that were methylated was 0.72 % for fertilized embryos and 0.62 % for thelytokous embryos, consistent with previous studies of the honey bee methylome [33, 44]. Taking into account the total genome-wide number of CG sites covered by 2–30 reads in both samples, the number of methylated CG sites was significantly greater in fertilized embryos than in thelytokous embryos (χ 21 = 1133.0, p <0.001).
We define methylated sites as those where the proportion of non-converted bisulfite reads was significantly higher (as determined by a binomial test) than the background bisulfite non-conversion rate (as measured in the lambda control). Ignoring sites that had less than two reads in either sample, 74,498 CG sites were significantly methylated in both samples (Table 1, Fig. 2a). There were 21,202 CG sites methylated in fertilized embryos but not thelytokous embryos, and 16,004 sites methylated in thelytokous embryos but not fertilized embryos (Table 1, Fig. 2a).
Across methylated CG regions, individual cytosines varied in methylation frequency, ranging from 10 to 100 % methylation (Fig. 3a). The methylated sites that were common to both fertilized and thelytokous embryos had a median methylation frequency of 83 % in thelytokous embryos, and 80 % in fertilized embryos (Fig. 3b). In contrast, the methylated sites that were unique to either fertilized or thelytokous embryos tended to have less methylation. In thelytokous embryos the median methylation frequency of sites showing some methylation was 67 %. In fertilized embryos the median methylation frequency was 60 % (Fig. 3b). This indicates that Capensis embryos share a core set of common CG sites that are highly methylated, regardless of whether the embryos arise sexually or thelytokously. Sites that are methylated only in fertilized or thelytokous embryos tend to have lower methylation frequency (Fig. 3b).
We identified over 50,000 cytosine methylation sites in a non-CG context (‘CH’ methylation: cytosines methylated adjacent to A, T or C; Additional file 1: Table S1). Some of these sites were artifacts arising from vector/mammalian contamination that are present in the Amel2.0 reference sequence (769 cytosines; Additional file 1: Table S1). To examine whether the remaining non-CG sites corresponded to actual Capensis CG sites that are absent from the Amel2.0 reference genome (thus artificially appearing as non-CG methylation), we examined the SNPs present from our Capensis2.0 alignment derived from Wallberg et al. . Only 361 (fertilized embryo) and 113 (thelytokous embryo) non-CG sites corresponded to SNPs in the Capensis2.0 reference sequence (Additional file 1: Table S1). It is likely that at least some non-CG methylated sites correspond to SNPs present in our source population that were not present in the Capensis population used by Wallberg et al. . Finally, a small number of false positives are expected to arise from incomplete bisulfite conversion and errors generated by the next generation sequencing process . Despite these caveats we observed a large number of non-CG methylated cytosines in fertilized and thelytokous embroys (28 % of methylated Cs in fertilized embryos, 12 % in thelytokous embryos, Additional file 1: Table S1). Note, however, that very few non-CG methylated sites (561) were present in both samples.
Differentially methylated cytosines
Comparison of significantly methylated sites between genomes that have variable coverage may yield false positives due to insufficient coverage in one of the genomes being compared. To avoid biased differential methylation calls based on coverage differences, we determined that the median coverage across methylated CG sites in fertilized and thelytokous embryos was five-fold (Additional file 1: Figure S1), so we restricted our comparison to sites with five or more reads in both samples. Using the R package methylKit , we determined the number of significantly differentially methylated CG sites between fertilized and thelytokous embryos. A total of 2,121 sites had significantly different methylation (Fig. 4). 2014 of those sites contained at least a 25 % difference in methylation frequency. If a differentially methylated CG site contained 25 % higher methylation in one sample compared to the other, such sites were considered hypermethylated. Fertilized embryos had significantly more hypermethylated sites (1412) compared to thelytokous embryos (602; χ 21 = 325.8, P <0.001, Fig. 4).
Genes containing differentially methylated sites
We examined the genomic context of differentially methylated sites and found that in both samples over 75 % of sites were located in exons and over 90 % in gene bodies (introns or exons, Table 2). We could therefore assign the majority of differentially methylated sites to annotated genes. As expected from the observation that there are significantly more hypermethylated CG sites in fertilized embryos, there was a significant bias towards hypermethylation of genes in fertilized embryos (χ 21 = 123.6, p <0.001, Table 2, Fig. 2b). The top genes containing at least four hypermethylated sites in either fertilized or thelytokous embryos, are shown in Table 3 (for a full list of all genes with hypermethylated sites see Additional file 2 ). Eight genes with four or more hypermethylated sites in fertilized and thelytokous embryos were present among the top 381 differentially methylated genes (DMGs) from our previous comparison of haploid eggs and sperm . In that study, these eight genes all had significantly higher methylation in haploid eggs compared to sperm. The 14th DMG in haploid eggs relative to sperm , Stan (GB18207) contained 7 hypermethylated sites in fertilized embryos (Table 3, Fig. 5). The top two DMGs in haploid eggs, Stoned B (GB17165) and Sap-30 (GB18386) , were both among the most hypermethylated genes in fertilized Capensis embryos with four and five hypermethylated sites, respectively (Table 3).
Bisulfite PCR validation of candidate DMGs
If genomic imprinting occurs in Capensis, we predict that genes inherited from males and females will be differentially methylated at some sites. If so, paternally methylated genes would be hypermethylated in fertilized (queen laid) embryos relative to unfertilized (worker-laid) embryos. Accordingly, we determined if patterns of methylation differed depending on whether they were inherited from a male or a female. When examining whole genome bisulfite sequencing data, it is difficult to trace methylation back to an original paternal or maternal allele due to the removal of many informative SNPs during conversion of unmethylated cytosines to thymine. The problem is exacerbated by short reads, which often uncouple any remaining informative SNPs from methylation haplotypes. Therefore, we used bisulfite PCR encompassing larger fragments of candidate genes to investigate allele-specific methylation. We selected three candidate genes that contained a number of hypermethylated sites in fertilized eggs, Stan (GB18207, Fig. 6), Stoned B (GB17165, Additional file 1: Figure S2) and Sap30 (GB18386, Additional file 1: Figure S3) to examine the level of methylation in the fertilized and thelytokous embryos collected from our two remaining colonies, 2 and 3. To determine the parent of origin of each allele, we sequenced the same gene in the fathering male. Thus for fertilized eggs we could identify which allele was paternal and which was maternal. In thelytokous eggs we could determine whether an allele was originally derived from the grandmother or the grandfather.
Sequences of Stan showed three distinct alleles in each colony, via informative G-to-A SNPs. Within a colony, fertilized and thelytokous embryos carried the same three alleles (Colony 2, Fig. 6b; Colony 3, Fig. 6c). Based on sequences of Stan obtained from the drone fathers used as sperm donors in colonies 2 and 3, we assigned the paternally-derived allele as allele P1, and inferred the maternally derived queen alleles as allele M2 and M3 for each colony (Fig. 6b & c).
The methylation patterns of each Stan allele varied considerably. In Colony 2, allele P1 showed the highest level of methylation, with an average of 51 % of the 31 possible CG sites methylated in fertilized embryos and 56 % in thelytokous embryos (Fig. 6b, Table 4). Allele M2 showed the lowest levels of methylation, with 10 % (fertilized) and 3 % (thelytokous), and M3 was intermediate with 34 % (fertilized) and 28 % (thelytokous). In Colony 3, P1 had intermediate methylation levels (fertilized 28 %, thelytokous 32 %), M2 the highest (fertilized 61 %, thelytokous 51 %) and M3 the lowest (fertilized 2 %, thelytokous 5 %; Table 4). Thus while the percentage of methylation and the location of methylated sites between alleles differed greatly, within-allele methylation was similar between fertilized and thelytokous samples in the two colonies (Fig. 6b&c).
The degree of methylation in Stoned B in Colony 2 and Colony 3 was much lower than the methylation present in the whole genome bisulfite sequence of Colony 1 fertilized and thelytokous eggs, and not all of the three alleles could be identified based on SNPs. Thus a comparison of allele-specific methylation could not be made for this gene, although there was some evidence for differential methylation in amplicons where two alleles could be distinguished (Additional file 1: Figure S2).
Due to lack of informative SNPs, we could not distinguish between the three alleles in Sap30, although for this gene methylation levels were similarly high to the whole genome bisulfite sequence in colony 1 (Additional file 1: Figure S3). In colony 2, three of the 32 CG sites present were significantly more highly methylated in fertilized embryos and one more highly methylated in thelytokous embryos; and in colony 3, two CG sites were significantly more highly methylated in fertilized embryos (Additional file 1: Figure S3).
If imprinting occurs in honey bees, we would expect to see a difference in the methylation patterns of bi-parental, fertilized embryos produced by union of an egg and sperm, compared to uni-parental, thelytokous embryos produced by parthenogenesis. Here we have shown that embryos derived from the union of a paternal and maternal genome had about 10 % more total methylated CG sites (Table 1), twice as many hypermethylated CG sites (Fig. 4) and twice as many hypermethylated genes (Fig. 2b) than genetically identical embryos that arose from the union of two maternal genomes. Candidate genes that are parentally-imprinted may therefore exist among the hypermethylated genes. This finding enhances earlier theoretical [18–20], phenotypic [34–36], and genomic  evidence that honey bees use their DNA methylation system to imprint certain genes in a parent-of-origin specific manner.
Biologically and technically replicated bisulfite sequencing of one of the top hypermethylated genes in fertilized eggs, Stan, revealed allele-specific methylation. The three alleles in each of colonies 2 and 3 had differing methylation patterns that were strongly correlated with genotype. This allele-specific methylation was the same between fertilized and theytokous embryos, and is therefore inconsistent with paternal or maternal imprinting: methylation was the same regardless of the sex of the parent of origin (Fig. 6).
While the most intensely studied examples of allele-specific methylation relate to parent-of-origin imprinting, it has been reported that the most widespread form of allele-specific DNA methylation in humans is determined by DNA sequence in cis [10, 66]. Sequence-dependent or cis-mediated allele-specific methylation has been mapped to the human genome  revealing many SNPs and genomic regions that are associated with variation in DNA methylation . While allele-specific methylation resulting from genomic imprinting is a non-Mendelian process, cis-mediated allele-specific methylation that is genetically determined is inherited in a Mendelian manner. An epigenetic mechanism that is genetically determined creates another level of genome regulation where the genotype and epigenotype are both involved in generating the phenotype.
The mechanisms behind genetically-determined allele-specific methylation are unclear . The allele-specific methylation observed in Stan included one G-to-A SNP that resulted in the removal of a CG site, thus physically interrupting the process of methylation. However the majority of methylation differences were not as clear and may be a result of long-range affects of polymorphisms that modify chromatin structure and affinity of DNA binding proteins .
Mono-alleleic methylation leading to allele-specific gene expression has previously been reported in ants  and bumblebees , and in both cases the authors hypothesised that differences were due to parent-of-origin effects. An alternate explanation is that the observed mono-alleleic methylation was due to cis-mediated allele-specific methylation, and that sequence-specific methylation is widespread in Hymenoptera. There is one other report of cis-mediated allele-specific methylation in honey bees , where differentially methylated obligatory epialleles of the AmLAM locus were correlated with sequence variation and resulted in different transcription levels. If allele-specific methylation is as widespread as it is in the human genome, it may have a major impact on the epigenetic contribution to phenotype in Hymenoptera.
Whole genome bisulfite sequence analysis indicated that Stan was hypermethylated in fertilized embryos compared to thelytokous embryos at multiple CG sites but bisulfite PCR validation in Colonies 2 and 3 did not support this conclusion. The difference in our genomic comparison may have arisen from unequal sampling of each of the three alleles present in Colony 1, such that more highly methylated alleles were over-represented in the genomic reads in our fertilized egg sample. In the whole genome bisulfite data, the mean coverage for Stan was 9.2 reads in fertilized embryos and 12 reads for thelytokous embryos. At this level of coverage, it is possible that the three alleles were unequally covered.
Whole genome bisulfite sequencing has previously been used to compare methylation patterns between social insect castes with contrasting results [44, 69, 72]. Allele specific methylation may confound genome-wide analysis of differential methylation in any species where there is cis-mediated allele-specific methylation. This kind of problem may be particularly acute in social insect species like honey bees where queens mate with multiple males. For example, if we compare groups of workers of different ages or tasks and conclude that their methylation patterns are different, our conclusion may be confounded by non-equal representation of alleles that have inherently different patterns of methylation. Recent evidence suggests that lack of biological replication could result in sample-specific, rather than caste-specific differences, leading to incorrect conclusions arising from genome-wide methylation patterns . Where possible, an experimental design that enables biological replication or alternative methods of validation is required to support conclusions of differential methylation. A combination of high-coverage analysis , single drone inseminated queens  and reciprocal crosses [34, 37] should allow specific alleles to be followed and provide a means to definitively distinguish between parent-of-origin and cis-mediated allele-specific methylation.
We have shown that global DNA methylation patterns differ between diploid embryos produced sexually and by thelytokous parthenogenesis. Our experiment minimized differences between embryos that were unrelated to uni- or bi-parental origin. These patterns of differential methylation between bi-parental and uni-parental eggs are consistent with imprinting. Nonetheless we caution that at least some of these patterns may be a consequence of cis-mediated allele specific methylation. In particular, our demonstration of allele-specific methylation in Stan means that not all epigenetic differences between between bi-parental and uni-parental embryos can be attributed to parent of origin effects.
Whole genome bisulfite sequencing data generated in this study have been deposited in Genbank BioProject under the following accession numbers: (http://www.ncbi.nlm.nih.gov/bioproject/PRJNA282063).
Studies were conducted using Apis mellifera capensis, an invertebrate animal that is neither endangered or protected and does not require ethics approval under the Animal Research Act 1985, NSW Government, Australia. The apiaries for sample collection were maintained at the ARC-Plant Protection Research Institute, Stellenbosch, South Africa. Samples were preserved in ethanol and imported under quarantine import permit IP14019459 to our quarantine approved facility (QAP #N2083) at the University of Sydney, Australia where they were processed in accordance with quarantine procedures.
Apis mellifera capensis
Differentially methylated gene
Single nucleotide polymorphism
Burt A, Trivers R. Genes in conflict: the biology of selfish genetic elements. Cambridge, Massachusetts: Harvard University Press; 2006.
Normark BB. Perspective: maternal kin groups and the origins of asymmetric genetic systems - genomic imprinting, haplodiploidy, and parthenogenesis. Evolution. 2006;60(4):631–42.
Haig D. The kinship theory of genomic imprinting. Annu Rev Ecol Syst. 2000;31:9–32.
Lawson HA, Cheverud JM, Wolf JB. Genomic imprinting and parent-of-origin effects on complex traits. Nat Rev Genet. 2013;14(9):609–17.
Pires ND, Grossniklaus U. Different yet similar: evolution of imprinting in flowering plants and mammals. F1000Prime Rep. 2014;1(6):63.
McGrath J, Solter D. Completion of mouse embryogenesis requires both the maternal and paternal genomes. Cell. 1984;37(1):179–83.
Surani MA, Barton SC, Norris ML. Development of reconstituted mouse eggs suggests imprinting of the genome during gametogenesis. Nature. 1984;308(5959):548–50.
Paulsen M, Ferguson-Smith AC. DNA methylation in genomic imprinting, development and disease. J Pathol. 2001;195(1):97–110.
Bird A. DNA methylation patterns and epigenetic memory. Genes Dev. 2002;16(1):6–21.
Meaburn EL, Schalkwyk LC, Mill J. Allele-specific methylation in the human genome: Implications for genetic studies of complex disease. Epigenetics. 2010;5(7):578–82.
Li E, Beard C, Jaenisch R. Role for DNA methylation in genomic imprinting. Nature. 1993;366(6453):362–5.
Foret S, Kucharski R, Pittelkow Y, Lockett GA, Maleszka R. Epigenetic regulation of the honey bee transcriptome: unravelling the nature of methylated genes. BMC Genomics. 2009;10:472.
Wojciechowski M, Rafalski D, Kucharski R, Katarzyna M, Maleszka J, Bochtler M, et al. Insights into DNA hydroxymethylation in the honeybee from in-depth analyses of TET dioxygenase. Open Biology. 2014;4(8):140110.
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(1):666.
Feng S, Cokus S, Zhang X, Chen P-Y, Bostick M, Goll MG, et al. Conservation and divergence of methylation patterning in plants and animals. PNAS. 2010;107(19):8689–94.
Zemach A, McDaniel IE, Silva P, Zilberman D. Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science. 2010;328:916–9.
Glastad KM, Hunt BG, Yi SV, Goodisman MAD. DNA methylation in insects: on the brink of the epigenomic era. Insect Mol Biol. 2011;20(5):553–65.
Drewell RA, Lo N, Oxley PR, Oldroyd BP. Kin conflict in insect societies: a new epigenetic perspective. TREE. 2012;27(7):367–73.
Haig D. Intragenomic conflict and the evolution of eusociality. J Theor Biol. 1992;156(3):401–3.
Queller DC. Theory of genomic imprinting conflict in social insects. BMC Evol Biol. 2003;3:15.
Bourke AFG. Worker reproduction in the higher Eusocial Hymenoptera. Q Rev Biol. 1988;63:291–311.
Oldroyd BP, Smolenski AJ, Cornuet J-M, Crozier RH. Anarchy in the beehive. Nature. 1994;371(6500):749.
Keller L, Reeve H. Genetic variability, queen number and polyandry in social hymenoptera. Evolution. 1994;48(3):694–704.
Engelstädter J. Constraints on the evolution of asexual reproduction. BioEssays. 2008;30(11-12):1138–50.
Rougier N, Werb Z. Minireview: parthenogenesis in mammals. Mol Reprod Dev. 2001;59(4):468–74.
Kono T, Obata Y, Wu Q, Niwa K, Ono Y, Yamamoto Y, et al. Birth of parthenogenetic mice that can develop to adulthood. Nature. 2004;428(6985):860–4.
Rabeling C, Kronauer D. Thelytokous parthenogenesis in eusocial Hymenoptera. Ann Rev Entomol. 2013;58:273–92.
Goudie F, Oldroyd BP. Thelytoky in the honey bee. Apidologie. 2014;45(3):306–26.
Goudie F, Allsopp MH, Solignac M, Beekman M, Oldroyd BP. The frequency of arrhenotoky in the normally thelytokous Apis mellifera capensis worker and the Clone reproductive parasite. Insect Soc. 2015;62(3):325–33.
Allsopp MH, Beekman M, Gloag RS, Oldroyd BP. Maternity of replacement queens in the thelytokous Cape honey bee Apis mellifera capensis. Behav Ecol Sociobiol. 2010;64(4):567–74.
Verma S, Ruttner F. Cytological analysis of the thelytokous parthenogenesis in the cape honeybee (Apis mellifera capensis Escholtz). Apidologie. 1983;14(1):41–57.
Oldroyd BP, Allsopp MH, Gloag RS, Lim J, Jordan LA, Beekman M. Thelytokous parthenogeneis in unmated queen honeybees (Apis mellifera capensis): central fusion and high recombination rates. Genetics. 2008;180(1):359–66.
Drewell RA, Bush EC, Remnant EJ, Wong GT, Beeler SM, Stringham JL, et al. The dynamic DNA methylation cycle from egg to sperm in the honey bee Apis melifera. Development. 2014;141(13):2702–11.
Oldroyd BP, Allsopp MH, Roth KM, Remnant EJ, Drewell RA, Beekman M. A parent-of-origin effect on honeybee worker ovary size. Proc R Soc B. 2014;281(1775):20132388.
Linksvayer TA, Rueppell O, Siegel A, Kaftanoglu O, Page Jr RE, Amdam GV. The genetic basis of transgressive ovary size in honeybee workers. Genetics. 2009;183(2):693–707.
Beekman M, Allsopp MH, Holmes MJ, Lim J, Noach-Pienaar L-A, Wossler TC, et al. Racial mixing in South African honeybees: the effects of genotype mixing on reproductive traits of workers. Behav Ecol Sociobiol. 2012;66(6):897–904.
Kocher SD, Tsuruda JM, Gibson JD, Emore CM, Arechavaleta-Velasco ME, Queller DC, et al. A search for parent-of-origin effects on honey bee gene expression. G3. 2015;5(10):1657–62.
Galbraith DA, Kocher SD, Glenn T, Albert I, Hunt GJ, Strassmann JE, et al. Testing the kinship theory of intragenomic conflict in honey bees (Apis mellifera). Proc Natl Acad Sci. 2016;113(4):1020–5.
Welch M, Lister R. Epigenomics and the control of fate, form and function in social insects. Curr Opin Insect Sci. 2014;1:31–8.
Maleszka R. Epigenetic integration of environmental and genomic signals in honey bees: the critical interplay of nutritional, brain and reproductive networks. Epigenetics. 2008;3(4):188–92.
Glastad KM, Chau LM, Goodisman MAD. Epigenetics in social insects. Adv Insect Physiol. 2015;48:227–69.
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. PNAS. 2012;109(13):4968–73.
Kucharski R, Maleszka J, Foret S, Maleszka R. Nutritional control of reproductive status in honeybees via DNA methylation. Science. 2008;319(5871):1827–30.
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(11):e1000506.
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(10):1371–3.
Wang X, Werren JH, Clark AG. Genetic and epigenetic architecture of sex-biased expression in the jewel wasps Nasonia vitripennis and giraulti. Proc Natl Acad Sci. 2015;112(27):E3545–54.
Hunt BG, Glastad KM, Yi SV, Goodisman MAD. Patterning and Regulatory Associations of DNA Methylation Are Mirrored by Histone Modifications in Insects. Genome Biol Evol. 2013;5(3):591–8.
Hunt BG, Brisson JA, Yi SV, Goodisman MAD. Functional conservation of DNA methylation in the pea aphid and the honeybee. Genome Biol Evol. 2010;2:719–28.
Sarda S, Zeng J, Hunt BG, Yi SV. The evolution of invertebrate gene body methylation. Mol Biol Evol. 2012;29(8):1907–16.
Glastad KM, Hunt BG, Yi SV, Goodisman MAD. Epigenetic inheritance and genome regulation: is DNA methylation linked to ploidy in haplodiploid insects? Proc R Soc B. 2014;281(1785):20140411.
Holmes MJ, Oldroyd BP, Allsopp MH, Lim J, Wossler TC, Beekman M. Maternity of emergency queens in the Cape honey bee, Apis mellifera capensis. Mol Ecol. 2010;19(13):2792–9.
Chapman NC, Beekman M, Allsopp M, Rinderer TE, Lim J, Oxley PR, et al. Inheritance of thelytoky in the honey bee Apis mellifera capensis. Heredity. 2015;114(6):584–92.
Beeler SM, Wong GT, Zheng JM, Bush EC, Remnant EJ, Oldroyd BP, et al. Whole genome DNA methylation profile of the jewel wasp (Nasonia vitripennis). G3. 2014;4(3):383–8.
Xi Y, Li W. BSMAP: whole genome bisulfite sequence MAPping program. BMC Bioinformatics. 2009;10:232.
The Honeybee Genome Sequencing Consortium. Insights into social insects from the genome of the honeybee Apis mellifera. Nature. 2006;443(7114):931–49.
Zayed A, Whitfield CW. A genome-wide signature of positive selection in ancient and recent invasive expansions of the honey bee Apis mellifera. PNAS. 2008;105(9):3421–6.
Wallberg A, Han F, Wellhagen G, Dahle B, Kawata M, Haddad N, et al. A worldwide survey of genome sequence variation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat Genet. 2014;46(10):1081–8.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-genertaion DNA sequencing data. Genome Res. 2010;20(9):1297–303.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57(1):289–300.
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(10):R87.
Warnecke PM, Stirzaker C, Song J, Grunau C, Melki JR, Clark SJ. Identification and resolution of artifacts in bisulfite sequencing. Methods. 2002;27(2):101–7.
Kucharski R, Foret S, Maleszka R. EGFR gene methylation is not involved in Royalactin controlled phenotypic polymorphism in honey bees. Sci Rep. 2015;5:14070.
O’Rawe JA, Ferson S, Lyon GJ. Accounting for uncertainty in DNA sequencing data. Trends Genet. 2015;31(2):61–6.
Tycko B. Allele-specific DNA methylation: beyond imprinting. Human Mol Genet. 2010;19(R2):R210–20.
Kerkel K, Spadola A, Yuan E, Kosek J, Jiang L, Hod E, et al. Genomic surveys by methylation-sensitive SNP analysis identify sequence-dependent allele-specific DNA methylation. Nat Genet. 2008;40(7):904–8.
Mendizabal I, Keller TE, Zeng J, Yi SV. Epigenetics and evolution. Integr Comp Biol. 2014;54(1):31–42.
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(19):1755–64.
Lee KD, Lonsdale ZN, Kyriakidou M, Nathanael D, Amarasinghe HE, Mallon EB. Monoallelic methylation and allele specific expression in a social insect. bioRxiv. 2015, (early edition).
Wedd L, Kucharski R, Maleszka R. Differentially methylated obligatory epialleles modulate context-dependent LAM gene expression in the honey bee Apis mellifera. Epigenetics. 2015, (in press).
Libbrecht R, Oxley Peter R, Keller L, Kronauer Daniel Jan C. Robust DNA Methylation in the Clonal Raider Ant Brain. Current Biology. 2016, (in press).
We thank Eliot Bush and Suzannah Beeler from Harvey Mudd College for bioinformatics advice, and members of the Behaviour and Genetics of Social Insects Laboratory, University of Sydney for helpful comments on the manuscript.
This work was funded by Australian Research Council grants DP120101915 to BPO, DE140100199 to AA and DP150100151 to BPO and AA; National Institute of Health (GM090167) and National Science Foundation (IOS-0845103) grants to RAD and a Howard Hughes Medical Institute Undergraduate Science Education Program grant (52007544) to the Biology department at Harvey Mudd College.
The authors declare that they have no competing interests.
EJR conceived and designed the study, performed molecular lab work, analysed the data and wrote the manuscript. AA and PB analysed data and wrote the manuscript. GB performed molecular lab work. MB and MHA collected samples and wrote the manuscript. CS and RAD provided resources and drafted the manuscript. BPO coordinated the study, designed the study, instrumentally inseminated the queens, collected samples and wrote the manuscript. All authors gave final approval for publication.
Figure S1. Comparison of fold coverage at methylated CG sites in thelytokous and fertilized embryos. Both samples exhibit a median coverage of 5 reads, thus our comparison of differential methylation was restricted to CG sites with at least 5 reads in both samples. Figure S2. Direct sequencing of bisulfite PCR products of Stoned B (GB17165) from colony 2 and 3 fertilized and thelytokous embryos. Red ovals indicated methylated cytosines, blue squares indicate SNPs. Figure S3. Direct sequencing of bisulfite PCR products of Sap30 (GB18386) from colony 2 and 3 fertilized and thelytokous embryos. Red ovals indicate methylated cytosines, open circles indicate total number of CG sites in the region, and dashed lines indicate incomplete sequence reads. Also shown is the methylation patterns in fertilized and thelytokous embryos from whole genome bisulfite sequencing of Colony 1. Asterisk indicates CG sites that are hypermethylated in fertilized embryos, dagger indicates CG sites hypermethylated in thelyokous embryos. Table S1. Non-CG methylation in Fertilised and Thelytokous embryos and comparison to Apis mellifera Capensis SNPs. Table S2. Primers used in bisulfite nested PCRs of Stan (Fig. 6), Stoned B (Additional file 1: Figure S2), Sap30 (Additional file 1: Figure S3, Syd and Pcl. (PPTX 607 kb)
Genes with hypermethylated sites (25 % > methylation difference) in Fertilized embryos. Genes with hypermethylated sites (25 % > methylation difference) in Thelytokous embryos. Genes with a combination of hypermethylated sites (25 % > methylation difference) from both Fertilized and Thelytokous embryos. (XLSX 60 kb)
About this article
Cite this article
Remnant, E.J., Ashe, A., Young, P.E. 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 17, 226 (2016). https://doi.org/10.1186/s12864-016-2506-8