Intronic Non-CG DNA hydroxymethylation and alternative mRNA splicing in honey bees

Whole-genome shotgun bisulfite sequencing validates that CG methylation is primarily located in the exons

We performed whole-genome shotgun sequencing of bisulfite modified DNA (BS-Seq) from AHB and EHB heads and obtained over 20× coverage of both genomes (Table  1; Additional file 1: Figure S1; see Methods). All four of the honey bee BS-Seq studies conducted to date — three published studies [6,8,9], and our experiments reported here — identified ~5-10× more CG methylation in exons than in introns. In the present study, there was ~6% CG methylation in exons and ~1% in introns and intergenic regions in EHB, compared with ~3% and ~0.3%, respectively in AHB (Table  2). Similarly, our re-analysis of the data from ref. [6] identified ~8% CG methylation in exons and ~0.3% methylation in introns in EHB (AHB was not studied; Additional file 2: Table S1). Only 15% (AHB) to 21% (EHB) of the CG methylation is symmetrically methylated in honey bees (Table  1), which is lower than what is observed in mammals (over 90%). Our analysis methods, which we call BS-Miner (See Methods), are sensitive to non-CG sites and identify hemi-methylated DNA using algorithms analogous to those by which heterozygous DNA sequences are identified in whole-genome sequences [20]. Recently, a BS-seq analysis tool called Bismark was developed that does not filter out non-CG methylation [21]. The amount of CHH and CHG methylation identified in the ref. [6] dataset by Bismark was approximately the same as the amount of CG methylation (516,148 versus 540,208, Additional file 2: Table S2), which is consistent with our analyses of our AHB and EHB datasets that show much more CHH and CHG methylation than previously reported.

Whole-genome shotgun bisulfite sequencing identifies non-CG modifications that are enriched in introns

As in the previous studies [6,8,9], we determined that CG methylation is primarily in the exons. However, a second finding from our analyses of both our data and the data from ref. [6] is relatively high levels of non-CG modifications (i.e., either ^mC or ^hmC) in bee heads. Surprisingly, our methods detected about 5-fold more CHH modifications (H = C, A, T) than CG methylation in both AHB and EHB DNA (Table  1). As with the CG methylation, we also saw more than twice as many non-CG modifications in EHB than AHB heads. About 2.5% of the total number of CHH sequences was modified in EHB and about 1.1% in AHB (Table  1).

In order to validate our finding of high levels of non-CG modifications, we reanalyzed previously published honey bee data from ref. [6] with our methods and again found that the amount of non-CG modifications exceeded the amount of CG methylation (Additional file 2: Table S1). The differences in the amounts of non-CG modifications in our data compared with ref. [6] might be caused by bisulfite treatment conditions (we did a single treatment and they did multiple treatments) or differences in strains of bees used in the two studies. We constructed the Illumina libraries using the identical protocol as ref. [6] which used unmethylated oligonucleotides, followed by amplification of only bisulfite converted oligonucleotides, to ensure that only fully bisulfite converted DNA is incorporated into the libraries (see Methods).

We also validated our data analysis pipelines by using two recent independently developed bisulfite methylation analysis program (BS-Map and Bismark) [21-23] on the data from ref. [6] and were again able to identify non-CG methylation sites (Additional file 2: Table S2 and Additional file 1: Figure S2). These independent algorithms also found more non-CG than CG methylation, consistent with our findings.

In contrast to CG methylation, CHH modifications were highest in introns (~4% and ~2% in EHB and AHB) and lowest in exons (~2% and ~0.8% in EHB and AHB) (Table  3). Consistent with this, our re-analysis of the EHB data from ref. [6] identified ~4 times more CHH modifications in introns than exons (2% versus 0.5%; Additional file 2: Table S1).

We also detected CHG modifications at levels lower than both CG and CHH (~1.1% and ~0.3% in EHB and AHB, Table  4). There was only about one-seventh as many CHG modifications as CHH modifications, and very few of the CHG modifications were symmetrical (~3.4% in EHB and ~2.3% in AHB, Table  2). CHG modifications were more uniform across the genome than CG methylation or CHH modifications (Table  2). Coverage for individual chromosomes (Additional file 1: Figure S3) demonstrates that there were no significant biases toward any portion of the genome in the sequencing procedure.

Overall, EHB had almost 2.5× more modified Cs than did AHB (1,591,851 versus 638,400). We also observed ~3-4-fold more CG methylation than in the previous three studies: we detected 253,041 methylated CGs in EHB, compared with 80,000-90,000 in the previous 3 studies (Table  2) [6,8,9]. This appears to be due to higher sensitivity of the analysis program used; as stated above, using our methods on data from ref. [6] identified 334,949 methylated CGs (Additional file 2: Table S1), which is even more than we identified in our data. The significance of EHB having 2.5× more modified Cs than AHB is not known.

Bees have 5-hydroxymethylcytosine

BS-Seq cannot distinguish ^mC from ^hmC because both base pair as C after BS treatment [10]. We used highly sensitive anti-CMS antibodies [13] to determine the levels of ^hmC in the heads and bodies of AHB and EHB. Consistent with the BS-Seq results, we found comparable and statistically indistinguishable levels of ^hmC in EHB and AHB heads (15.2 pmol/μg and 13.5 pmol/μg of genomic DNA) (Additional file 1: Figure S4). For bodies, there was significantly more ^hmC in EHB than AHB (25.7 pmol/μg and 19.3 pmol/μg) (p < 0.05, 2-tailed t-test, Additional file 1: Figure S4).

The 5-hydroxymethylcytosine in bees is enriched in introns

To distinguish ^mC from ^hmC, we immunoprecipitated honey bee head DNA with antibodies against 5-hydroxymethylcytosine (HMeDIP). The immunoprecipitated DNA was then sequenced by next-generation DNA sequencing (HMeDIP-seq). Compared to previous findings that ^mC is found primarily in exons, we found that most of the ^hmC DNA is present in introns (Table  5), where most of the non-CG modifications are also present. This leads to the speculation that ^hmC is enriched in non-CG sites, and that many of the non-CG modifications that are detected by whole genome bisulfite sequencing are actually ^hmC, since bisulfite cannot distinguish ^mC from ^hmC.

Pvu-seq validates the location of ^hmC in introns

To validate the ^hmC findings, we developed a new technique that we call Pvu-Seq. This technique involves digesting the DNA isolated from AHB and EHB heads with the Type 2 restriction endonuclease, PvuRts1I, which cuts near single hydroxymethylcytosine sites [24-26]. The distances between the cleavage sites and the modified cytosine are fixed within a narrow range, with the majority being 11-13 nucleotides away in the top strand and 9-10 nucleotides away in the bottom strand [24,25]. There was an excellent correlation between hmDIP-Seq and Pvu-Seq data; over 89% of the HMeDIP-Seq peaks were also represented by Pvu-Seq peaks. An example of an HMeDIP-Seq peak correlating with a Pvu-Seq reads in AHB is shown in Figure  1b and c. These results confirm the findings made with other techniques and indicate that Pvu-Seq is a valid technique for mapping ^hmC sites in the AHB and EHB genomes.

There are more non-CG modifications in genes with a low CG content

Previous analyses in honey bees have shown that there are two classes of genes with respect to CG methylation: one has a low observed over expected (o/e) CG ratio (i.e., low CG content), is highly methylated, and is enriched in housekeeping genes, and a second has a high o/e ratio (i.e., high CG content), is unmethylated, and is enriched in caste-specific and developmental genes (Figure  2, dashed lines; Additional file 1: Figure S5) [6,8,27-29]. Although non-CG sequences have a unimodal distribution (Additional file 1: Figure S6), rather than a bimodal distribution in the genome, non-CG modifications in introns surprisingly were found primarily in genes with a low o/e CG ratio (Figure  2). Genes with greater than 10% non-CG modifications, such as ^mC and ^hmC, in introns are primarily in low o/e CG genes (i.e., the left peak in the o/e CG ratio plot, dashed lines; Figure  2), whereas genes with zero percent non-CG modifications in introns are in high o/e CG genes (i.e., the right peak in the o/e CG ratio plot, dashed lines; Figure  2). Therefore, the pattern of methylation is dependent on the CG dinucleotide content and not the CHH or CHG trinucleotide content. We speculate that this might indicate that CG methylation is somehow linked to non-CG modifications, possibly via interactions between the maintenance DNA methyltransferase, Dnmt1, and the de novo enzyme, Dnmt3, both of which are present in bees [5], and the bee TET protein. In contrast to Dnmt1, which has almost exclusive specificity for CG sites (however, see ref. [30]), Dnmt3 is responsible for most non-CG methylation in hESC [31]; this has not been examined in bees.

CHH modifications are enriched in the introns of genes that regulate transcription

We found that CHH modifications are highest in the introns of genes in the GO category “regulation of transcription” for both AHB and EHB (Figure  3a). This includes several Homeobox transcription factors, such as Distalless (Figure  3b) and aristaless (not shown). This is in contrast to genes with the highest CG methylation in exons, which were enriched in “housekeeping gene” GO categories, such as mitochondrial, ribosomal, and nucleotide-binding genes ( [6,8,9] and data not shown). As shown for Distalless (Figure  3b), introns in the “regulation of transcription” GO category often had a large amount of CG methylation in addition to the non-CG modifications. This again suggests that CG methylation and non-CG modifications are coordinately regulated.

Non-CG modifications might regulate alternative mRNA splicing

Consistent with the idea that DNA methylation might be involved in regulating mRNA splicing [8,9], we found that splice donors and acceptors were often encoded by DNA with non-CG modifications, such as either ^mC or ^hmC. In bees and other invertebrates, over 90% of splice donors have a G at the first position and a U at the second position (i.e., 5’-AC-3’ on the template DNA strand, where the C on the template strand can be methylated). We identified several hundred modified cytosines at splice donor and acceptor sites on the template strands (346 ^mCs in 321 genes in AHB, and 1677 ^mCs in 1312 genes in EHB) (Figure  4a).

Based on the above numbers, only ~0.66% of the ~56,000 splice junctions were methylated in AHB (346) and ~3.3% in EHB (1627) (Figure  4a). However, the distribution was clearly not random because pathway analyses show that genes with methylated splice sites were most enriched in the GO pathway “phosphoprotein” in both AHB and EHB (FDR < 10E-9 for both; Additional file 2: Table S5). Intriguingly, the GO pathway “alternative splicing” was also significantly enriched in both AHB and EHB (FDR < 0.05 for both; Additional file 2: Table S5). For example, the honey bee gene that was most heavily methylated at splice junctions is GB13778, the ortholog of Drosophila dumpy, whose protein products are involved in cell adhesion in D. melanogaster[33]; it has four methylated splice junctions in AHB and twelve methylated splice junctions in EHB at CHH and CHG sites (Figure  4b). Since dumpy has complex splicing programs in Drosophila (16 distinct mature spliced mRNAs are listed in FlyBase), and there are dozens of dumpy exons in honey bees, it is attractive to speculate that non-CG modifications at splice junctions in honey bees regulate alternative splicing at this locus and others as well.

In hESC lines, methylated non-CG sites are strongly conserved especially within the motif 5’-TA^mCAG-3’ on the non-coding DNA strand at 3’ splice junctions [2]. While both hESCs and bees have non-CG modifications at splice junctions, bees differ from hESCs in several respects. In hESCs, methylation is symmetrical at CAG sites on both the template and non-template strand only at the 3’ splice junction, which is most frequently CAG. However in bees, methylation was primarily asymmetrical at CHH sites at the +1 position on the template strand encoding splice acceptors and the -1 position on the template strand at splice donors, and very few of the CHG sites in bees were symmetrically modified (Table  2). Another difference between hESC and bees is that the 3’ splice sites are predominantly methylated in humans but both 5’ and 3’ splice sites were modified in bees (Figure  4).

Genes with more CHH modifications in EHB than AHB are enriched in behavioral response genes

Genes with significantly more CHH modifications in EHB than AHB were enriched in GO categories that are involved in neurological functions, such as “response to external stimulus”, “substrate specific channel activity”, “exocytosis”, and “neurotransmitter receptor activity” (Additional file 2: Table S6). These categories were highly significant even after correcting for the higher overall CHH modifications in EHB as well as multiple testing using the false discovery rate method (Additional file 2: Table S6) [20]. It is attractive to speculate that differential CHH modifications in introns might partially explain the striking behavioral differences between AHB and EHB, especially in aggression, but this is beyond the scope of the present study. Genetic studies suggest epigenetic regulation of aggression in vertebrates [34-38], and there are extensive aggression-related differences in brain gene expression between the aggressive AHB and the less-aggressive EHB [39]. Since there also are developmental differences between AHB and EHB, in addition to differences in aggressive behavior, these two types of differences would need to be teased apart in future studies.

Gene expression positively correlates with both ^mC and ^hmC levels in the exons

We found a weak, but significant correlation between exon methylation and exon expression. This result was obtained for methylation detected by either bisulfite sequencing or Pvu-Seq (Additional file 1: Figure S9). This correlation was stronger for CG methylation than non-CG.

Exons and Splice junctions with ^mC or ^hmC appear to affect alternative mRNA splicing

As stated above, several studies have suggested that DNA methylation might be involved in regulating mRNA splicing. To determine how our new ^hmC finding might influence our understanding of the relationship between DNA methylation and mRNA splicing, we analyzed exons and splice junctions that differ in either ^mC or ^hmC between AHB and EHB and compared this with RNA-Seq data that we generated from AHB and EHB brains. We found several examples of differential ^mC and ^hmC associated with alternative mRNA splicing, as described below.

Consistent with previous reports, we found several examples of differential CG methylation associated with alternative mRNA splicing (Figure  5a). Consistent with previously published examples [6], CG methylation in DNA encoding exons often correlated with that exon being skipped. For example, for gene GB15706, we observed that a heavily methylated exon in AHB was skipped, whereas an adjacent heavily methylated exon in EHB was skipped (see Additional file 1: Figure S8). We also report for the first time that non-CG modifications, such as ^hmC, also showed a correlation with alternative splicing (Figure  5b). For example, the gene GB18247 had ^hmC on an exon in AHB and that exon was retained in AHB, but that same exon did not have ^hmC on it in EHB and that exon was skipped. In other words, at least for these examples, ^mC on exons correlated with exon skipping, whereas ^hmC on exons correlated with exon retention.

Sequencing

All sequencing was performed using Illumina Genome Analyzer GAIIx with a Paired End Cluster Generation Kit. Image analysis, base calling and sequence extraction was performed using standard Illumina Pipeline v1.6 software. We performed whole-genome shotgun sequencing of bisulfite modified DNA (BS-Seq). DNA sequencing (>20× coverage) also was performed to ensure that the C to U conversions were BS induced, rather than natural single nucleotide polymorphisms (SNPs; not presented here). BS-Seq was performed on Africanized honey bees (AHB) (Apis mellifera scutellata) and European honey bees (EHB) (a mixture of subspecies, primarily A.m. ligustica). The number of lanes sequenced was 11 (pair-end reads) for AHB and 8 for EHB, resulting in 240 million reads for AHB and 317 million reads for EHB. Reads were 76 base-pair long yielding a total of 18.2 giga bases and 24.1 giga bases for AHB and EHB respectively. We note that we used the Illumina Whole Genome Bisulfite Sequencing (WGBS) kit which first ligates non-methylated primers to the genomic DNA prior to bisulfite conversion. After bisulfite conversion, a second set of primers are used that only amplify the fully converted primers (e.g., [46]). We confirmed that only fully-converted primers were amplified in the BS-seq libraries (not shown).

Bioinformatics analysis

The bioinformatics analysis was conducted by two different groups (D.M.R. and S.Z.), using different approaches and without sharing any processed data or results.

Method 1 (BS-Miner)

The reference genome was Amel2, which is ~228 million bases long. At the time we performed our first bisulfite sequencing data analysis, there were only a few BS-Seq mappers available and some of them had a tendency to filter out non-CG methylation. For that reason, we decided to create our own unbiased analysis pipeline, which we called BS-Miner. It should be noted that nowadays many more options exist and those early mappers have been greatly improved, so there is no longer need to develop ad hoc methods. Read mapping and downstream analysis was done using BS-Miner. One of our main biological questions was whether non-CG methylation was present or not so we designed our pipeline using statistical methods well known in standard base calling algorithms. As a result, we obtain better sensitivity but only on 100%, 50% or 0% methylation levels (that is: methylated, hemi-methylated or no methylation). This sensitivity comes at a cost. As expected, the method does not detect methylation as a continuous range (from 0% to 100%) like other algorithms do. This design trade-off was aligned with our research hypothesis.

BS-Miner uses either BWA [47,48] or Bowtie [49] for read alignment. Both alignment programs are based on the Burrows-Wheeler [50] transformation and create SAM output format [51]. There are other tools and methods based on similar approaches [21,22,52]. In this case BWA was selected as the main mapping method in order to have better alignment near insertions and deletions. BS-Miner performs methylation calls by invoking Samtools [51], which uses a probabilistic model [53]. It must be noted that the BAQ model is explicitly disabled by BS-Miner, since some of its assumptions do not apply to methylation calls.

The BcfTools package was invoked to produce methylation calls in VCF format. BS-Miner was set to filter out low quality (Q < 20) methylation calls. After all mapping and filtering steps, the mean coverage was 20.8 and 27.4 for AHB and EHB respectively, which is > 2-fold more coverage than the previously most comprehensive bee study [5]. As a final step, BS-Miner performs several statistical analyses of the methylation results, including ranking of hypo-methylated and hyper-methylated genes by means of the Wilcoxon rank test and Fisher’s exact test. Multiple testing was corrected using False Discovery Rate methodology [20]. Some additional statistics were carried out using custom programs in R programming language [54] Further statistics from BS-Miner as well as additional analysis are available at http://www.mcb.mcgill.ca/~pcingola/bees/.

We also performed reanalysis of our data using Bowtie and a pipeline similar to the one shown in Krueger et al. [55]. Read trimming was performed using Trimmomatic [56]. Quality control was performed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and in-house software. Performing reanalysis using different pipelines (BS-Miner and Bismark [21]) and different filtering strategies, we obtained consistent results. Removing duplicates using Bismark’s remove duplicates module [21], did not seem to change our results significantly.

Method 2 (BSMap)

Using standard software, BSMap [23], we obtained results consistent with previous studies: methylation was primarily at CG dinucleotides in exons and very little non-CG methylation was present (Additional file 1: Figure S1a). There were 61,149,121 uniquely mapped cytosines (Cs) in EHB and 53,443,185 in AHB. More than 88% of EHB Cs and 83% of AHB Cs were covered by at least two sequencing reads (Additional file 1: Figure S2). However, when we included the reads that have DNA methylation in a non-CpG context, we found considerable amounts of CHH and CHG methylation (Additional file 1: Figure S1b). Bisulfite reads were mapped with BSMap to distinguish the methylated cytosine from unmethylated cytosines. The reads of bisulfite converted DNA were mapped to Apis mellifera genome assembly 4 after converting the Cs to Ts. Two mismatches were allowed for the alignment to be made. To reduce the potential erroneous cytosine methylation reads, the read for a cytosine was discarded if there was a mismatch event within the 2-bp surrounding context. There was also another “CHH-filter” which filtered out the entire read if the read contains three consecutive methylated CHH. The analyses are presented in two ways: with CHH-filter and without. For every gene, its 3k upstream region, 3k downstream region, and every exon or intron was divided into 30 bins and the ratio of the number of methylated cytosines over the number of all cytosines was plotted against the bin number. Bins were graphed from upstream (relative to the gene) to downstream.

Comparison of methylation in AHB and EHB

The cytosine methylation in 3k upstream of all genes was averaged to calculate "upstream" methylation levels for every gene in both AHB and EHB; the cytosine methylation levels in gene body (or exons) were averaged to calculate the methylation level of gene body (or exons). To reduce bias from low mapping coverage, genes with less than 100 cytosines covered by BS-Seq reads, in gene body and 3k upstream, were excluded.

Validating BS-Seq with MeDIP-Seq and HMeDIP-Seq analyses

Methylated DNA immunoprecipitation followed by sequencing (MeDIP-Seq or HMeDip-Seq) was performed for a total of 41.3 million 76 bp reads. Reads were aligned using BWA and SamTools. Peak-calling was performed using MACS 1.4 [57] beta version. Genomic DNA from 3 AHB and 3 EHB heads was sheared to 300-600 bp fragments, gel purified, immunoprecipitated with antibody, ligated to the library primers, amplified, and then sequenced. The mC antibodies were mouse monoclonal (Active Motif, Inc.) and the hmC antibodies were rabbit polyclonal (Active Motif, Inc.). Immunoprecipitation was with protein G beads (Active Motif, Inc.) following the manufacturer’s protocol. Additional file 1: Figure S4a shows the peak model based on the “forward reads” [those that align with the “Watson” (plus) strand] and the “reverse reads” [those that align with the “Crick” (minus) strand].

Differential methylation analyses

Counts of methylation sites per gene, transcript, intron, splice sites and other regions of interest were calculated. Gene orthologs were calculated by InParanoid [58,59]. Enrichment was calculated using a greedy Wilcoxon rank sum method, RssGsc (http://www.rssgsc.sourceforge.net), which also performs multiple testing correction using false discovery rate. Fisher’s exact test was also applied as a secondary method, by setting a suitable threshold in the ranked list.

Differencial SNP analysis

Counts of SNPs sites per gene, transcript, intron, splice sites and other regions of interest were calculated. Quantile normalization was applied and ranked genes analyzed for enrichment using the same methods as described in the previous section.

Observed over expected (o/e) analysis

The observed over expected ratio is defined as the number of CG dinucleotides in the reference sequence, divided by the number expected under a uniform random distribution. For genes having multiple transcripts (which are a minority in amel2 reference genome), a weighted average was calculated according to each transcript’s length. This definition is consistent with [60].

Accession numbers for data

The next-generation DNA sequencing experiments were done by the Applied Genomics Technology Core at Wayne State University. The next-generation DNA sequencing data (RNA-seq, Pvu-seq, MeDIP-seq, BS-seq, and HMeDIP-seq) was deposited into the GEO database according to the MINSEQE standards (Minimum Information about a high-throughput SeQuencing Experiment). The GEO database accession number is GSE50990. Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/projects/geo/).

Animal use ethical use issues

Honey bees are not a regulated invertebrate. Therefore, no ethical use approval is necessary.

GEO accession number

(GSE50990; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50990).