BACKGROUND
methylation of dna is increasingly appreciated as a potent way to regulate gene expression, but a comprehensive understanding of the diversity of methylation mechanisms has not yet been achieved.
for example, methylation that does not occur at cytosine-guanosine dinucleotide sequences  is an underappreciated and poorly understood form of epigenetic regulation.
while rare in mammalian somatic cells, non‒cg methylation occurs on 20‒40% of the cytosines in human embryonic stem cells
<cit> , and is thought to be involved in pluripotency
<cit> .
a recent comparative analysis of dna methylation across hesc lines found that heavily methylated non‒cg sites are strongly conserved, especially within the motif tamcag at 3’ splice junctions
<cit> , suggesting a role in splicing or alternative splicing of mrna transcripts.
also, ctcf induced rna polymerase ii pausing was shown to link alternative mrna splicing to dna methylation at a ctcf binding site in genomic dna encoding an intron
<cit> .
recently it was shown that rna interference knockdown of dna methyltransferase  <dig>  affects gene alternative splicing in the honey bee
<cit> .
the honey bee  is emerging as a new model to study effects of methylation on genome function because, unlike drosophila melanogaster, it possesses a fully functioning methylation system
<cit> .
three studies
<cit>  have reported that honey bees have cg methylation primarily at exon coding regions, and we have confirmed these studies here.
these three studies have also reported that honey bees have little, if any, non-cg methylation
<cit> .
however, since all three of the previous studies’ experimental design filtered out much of the non-cg methylation , it is still an open question as to whether there are significant amounts of non-cg methylation in bees.
the are several reasons for filtering out non-cg methylation:  non-cg methylation is much less abundant than cg methylation in mammals;  there are several times more non-cg sequences  than cg sequences, and focusing analyses on cg methylation is simpler;  non-cg methylation is often in less complex regions of the genome, such as introns, and is therefore difficult to map with short-read next-generation sequencing technology; and  since bisulfite works less well on double-stranded dna than single-stranded dna, less complex regions might form snap-back structures that are resistant to bisulfite conversion.
hydroxymethylation of dna is a newly discovered form of epigenetic regulation.
it has been found recently in embryonic stem cells  and in the brains of mammals
<cit> .
in mammals, tet proteins have been shown to be dioxygenases that convert mc to 5-hydroxymethylcytosine .
honey bees have a tet ortholog
<cit> , but hmc has not yet been reported.
genome-wide wide mapping in mouse escs has revealed that hmc is enriched at the start sites of genes whose promoters bear histone  <dig> lysine  <dig> trimethylation  and histone  <dig> lysine  <dig> trimethylation  marks
<cit> .
in human escs, this dual mark is derived from separable subpopulations of self-renewing and lineage-biased escs within the heterogeneous unfractionated esc population
<cit> .
the most common chemical approach to study dna methylation is treating single-stranded dna with bisulfite, but bisulfite cannot distinguish mc from hmc because both base-pair as cytosine after bisulfite treatment.
the modification mc mostly remains in this form after bisulfite treatment, whereas hmc is converted to cytosine methylene sulfonate  after bisulfite treatment, which has the same hydrogen bond donor and acceptor configuration as cytosine for base pairing to guanine
<cit> .
in this paper, in addition to whole-genome shotgun bisulfite sequencing of honey bee head dna, we developed new biochemical and bioinformatics tools to analyze non-cg methylation and hydroxymethylation.
we sharpened our analyses by comparing bees endemic to north america that are derived from a mixture of european subspecies of apis mellifera  with bees derived from the african subspecies apis mellifera scutallata, introduced to south america in the last century .
ehb and ahb are attractive candidates for comparative molecular analysis because they differ in many physiological and behavioral traits, including aggression.
differences in methylation between ehb and ahb that are reported here may be related to these physiological and behavioral differences, and could motivate further studies, beyond the scope of the present paper.
we show that non-cg methylation and hydroxymethylation are both enriched in introns, and cytosine modification at splice junctions might help regulate alternative splicing and gene expression.
RESULTS
we asked whether non-cg methylation and hydroxymethylation occurs in honey bees using unbiased approaches that do not filter out non-cg methylation.
we compared ahb and ehb mainly to demonstrate the robustness of our experimental and bioinformatics approaches.
to standardize our comparisons, we analyzed head methylomes only from ahb and ehb “guard” bees, a specialized group of individuals that patrol the hive entrance and respond first to a threat to their colony.
most of the tissues in the bee head consist of brain.
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  from ahb and ehb heads and obtained over 20× coverage of both genomes .
all four of the honey bee bs-seq studies conducted to date — three published studies
<cit> , 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 ~ <dig> %, respectively in ahb .
similarly, our re-analysis of the data from ref.
<cit>  identified ~8% cg methylation in exons and ~ <dig> % methylation in introns in ehb .
only 15%  to 21%  of the cg methylation is symmetrically methylated in honey bees , which is lower than what is observed in mammals .
our analysis methods, which we call bs-miner , 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
<cit> .
recently, a bs-seq analysis tool called bismark was developed that does not filter out non-cg methylation
<cit> .
the amount of chh and chg methylation identified in the ref.
<cit>  dataset by bismark was approximately the same as the amount of cg methylation , which is consistent with our analyses of our ahb and ehb datasets that show much more chh and chg methylation than previously reported.
*the total number of cytosines represent both the uniquely mapped and the non-uniquely mapped cytosines .
whole-genome shotgun bisulfite sequencing identifies non-cg modifications that are enriched in introns
as in the previous studies
<cit> , 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.
<cit>  is relatively high levels of non-cg modifications  in bee heads.
surprisingly, our methods detected about 5-fold more chh modifications  than cg methylation in both ahb and ehb dna .
as with the cg methylation, we also saw more than twice as many non-cg modifications in ehb than ahb heads.
about  <dig> % of the total number of chh sequences was modified in ehb and about  <dig> % in ahb .
in order to validate our finding of high levels of non-cg modifications, we reanalyzed previously published honey bee data from ref.
<cit>  with our methods and again found that the amount of non-cg modifications exceeded the amount of cg methylation .
the differences in the amounts of non-cg modifications in our data compared with ref.
<cit>  might be caused by bisulfite treatment conditions  or differences in strains of bees used in the two studies.
we constructed the illumina libraries using the identical protocol as ref.
<cit>  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 .
we also validated our data analysis pipelines by using two recent independently developed bisulfite methylation analysis program
<cit>  on the data from ref.
<cit>  and were again able to identify non-cg methylation sites .
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  and lowest in exons  .
consistent with this, our re-analysis of the ehb data from ref.
<cit>  identified ~ <dig> times more chh modifications in introns than exons .
we also detected chg modifications at levels lower than both cg and chh .
there was only about one-seventh as many chg modifications as chh modifications, and very few of the chg modifications were symmetrical .
chg modifications were more uniform across the genome than cg methylation or chh modifications .
coverage for individual chromosomes  demonstrates that there were no significant biases toward any portion of the genome in the sequencing procedure.
overall, ehb had almost  <dig> × more modified cs than did ahb .
we also observed ~3-4-fold more cg methylation than in the previous three studies: we detected  <dig>  methylated cgs in ehb, compared with  <dig> - <dig>  in the previous  <dig> studies
<cit> .
this appears to be due to higher sensitivity of the analysis program used; as stated above, using our methods on data from ref.
<cit>  identified  <dig>  methylated cgs , which is even more than we identified in our data.
the significance of ehb having  <dig> × 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
<cit> .
we used highly sensitive anti-cms antibodies
<cit>  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  .
for bodies, there was significantly more hmc in ehb than ahb  .
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 .
the immunoprecipitated dna was then sequenced by next-generation dna sequencing .
compared to previous findings that mc is found primarily in exons, we found that most of the hmc dna is present in introns , 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  <dig> restriction endonuclease, pvurts1i, which cuts near single hydroxymethylcytosine sites
<cit> .
the distances between the cleavage sites and the modified cytosine are fixed within a narrow range, with the majority being 11- <dig> nucleotides away in the top strand and 9- <dig> nucleotides away in the bottom strand
<cit> .
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  cg ratio , is highly methylated, and is enriched in housekeeping genes, and a second has a high o/e ratio , is unmethylated, and is enriched in caste-specific and developmental genes
<cit> .
although non-cg sequences have a unimodal distribution , 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 .
genes with greater than 10% non-cg modifications, such as mc and hmc, in introns are primarily in low o/e cg genes , whereas genes with zero percent non-cg modifications in introns are in high o/e cg genes .
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, dnmt <dig>  and the de novo enzyme, dnmt <dig>  both of which are present in bees
<cit> , and the bee tet protein.
in contrast to dnmt <dig>  which has almost exclusive specificity for cg sites , dnmt <dig> is responsible for most non-cg methylation in hesc
<cit> ; 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 .
this includes several homeobox transcription factors, such as distalless  and aristaless .
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 .
as shown for distalless , 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.
<cit> .
b, distalless is a hox gene and a transcription factor with high amounts of chh modifications in the introns of both ahb and ehb.
non-cg modifications might regulate alternative mrna splicing
consistent with the idea that dna methylation might be involved in regulating mrna splicing
<cit> , 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 .
we identified several hundred modified cytosines at splice donor and acceptor sites on the template strands  .
based on the above numbers, only ~ <dig> % of the ~ <dig>  splice junctions were methylated in ahb  and ~ <dig> % in ehb  .
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 .
intriguingly, the go pathway “alternative splicing” was also significantly enriched in both ahb and ehb .
for example, the honey bee gene that was most heavily methylated at splice junctions is gb <dig>  the ortholog of drosophila dumpy, whose protein products are involved in cell adhesion in d. melanogaster <cit> ; it has four methylated splice junctions in ahb and twelve methylated splice junctions in ehb at chh and chg sites .
since dumpy has complex splicing programs in drosophila , 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’-tamcag-3’ on the non-coding dna strand at 3’ splice junctions
<cit> .
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 + <dig> position on the template strand encoding splice acceptors and the - <dig> position on the template strand at splice donors, and very few of the chg sites in bees were symmetrically modified .
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 .
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” .
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
<cit> .
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
<cit> , and there are extensive aggression-related differences in brain gene expression between the aggressive ahb and the less-aggressive ehb
<cit> .
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 .
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 .
consistent with previously published examples
<cit> , cg methylation in dna encoding exons often correlated with that exon being skipped.
for example, for gene gb <dig>  we observed that a heavily methylated exon in ahb was skipped, whereas an adjacent heavily methylated exon in ehb was skipped .
we also report for the first time that non-cg modifications, such as hmc, also showed a correlation with alternative splicing .
for example, the gene gb <dig> 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.
1: figure s <dig>  b, top blue lines indicated alternatively spliced exons in ahb and ehb, the red-and-blue stacked lines underneath the exons show pvu-seq analyses of ahb and ehb heads.
note that hmc also correlates with alternative mrna splicing, but in manners that differ from mc .
details of the rna-seq analyses are given in additional file
1: figure s <dig>
CONCLUSIONS
our findings underscore the diversity of dna methylation mechanisms that exist.
non-cg modifications were only recently discovered in hesc and now we report them in the distantly related honey bee.
we speculate that cytosine methylation at exons and splice junctions on the dna might affect the mrna splicing machinery.
it is important to learn how these mechanisms work in relation to known regulators of splicing, such as histone acetylation
<cit> , histone  <dig> lysine  <dig> methylation
<cit> , histone  <dig> lysine  <dig> methylation
<cit> , and histone  <dig> lysine  <dig> methylation
<cit> .
understanding how non-cg hydroxymethylation might affect alternative splicing is an exciting new area of research.
our data are consistent with a model in which dna is methylated at cg sites by the maintenance dna methyltransferase, dnmt <dig>  at exons, and at non-cg sites by the de novo dna methyltransferase, dnmt <dig>  at introns.
in contrast to dnmt <dig>  mammalian dnmt <dig> is known to methylate both cg and non-cg sites in mammalian stem cells
<cit> .
we speculate that the honey bee tet enzyme primarily recognizes the non-cg sites in the introns, thereby enriching dna hydroxymethylation in the introns.
we further speculate that the mrna splicing machinery, as well as the histone modification machinery, distinguishes exons and introns by somehow recognizing the patterns of cg methylation in exons and non-cg hydroxymethylation in introns.
our identification of go categories related to protein phosphorylation that were enriched for genes with methylated splice junctions is consistent with a similar finding in a recent study of species-specific alternative exons
<cit> .
the authors present evidence that argues that alternative splicing is used to alter protein phosphorylation, which can alter protein stability, subcellular localization, activity, and other properties
<cit> .
further research is needed to determine the mechanism by which splice junction methylation and hydroxymethylation affect mrna splicing.

 methylation of dna is increasingly appreciated as a potent way to regulate gene expression, but a comprehensive understanding of the diversity of methylation mechanisms has not yet been achieved. for example, methylation that does not occur at cytosine-guanosine dinucleotide sequences is an underappreciated and poorly understood form of epigenetic regulation. while rare in mammalian somatic cells, non‒cg methylation occurs on 20‒40% of the cytosines in human embryonic stem cells, and is thought to be involved in pluripotency. a recent comparative analysis of dna methylation across hesc lines found that heavily methylated non‒cg sites are strongly conserved, especially within the motif tamcag at 3’ splice junctions, suggesting a role in splicing or alternative splicing of mrna transcripts. also, ctcf induced rna polymerase ii pausing was shown to link alternative mrna splicing to dna methylation at a ctcf binding site in genomic dna encoding an intron. recently it was shown that rna interference knockdown of dna methyltransferase affects gene alternative splicing in the honey bee. the honey bee is emerging as a new model to study effects of methylation on genome function because, unlike drosophila melanogaster, it possesses a fully functioning methylation system. three studies have reported that honey bees have cg methylation primarily at exon coding regions, and we have confirmed these studies here. these three studies have also reported that honey bees have little, if any, non-cg methylation. however, since all three of the previous studies’ experimental design filtered out much of the non-cg methylation, it is still an open question as to whether there are significant amounts of non-cg methylation in bees. the are several reasons for filtering out non-cg methylation: non-cg methylation is much less abundant than cg methylation in mammals; there are several times more non-cg sequences than cg sequences, and focusing analyses on cg methylation is simpler; non-cg methylation is often in less complex regions of the genome, such as introns, and is therefore difficult to map with short-read next-generation sequencing technology; and since bisulfite works less well on double-stranded dna than single-stranded dna, less complex regions might form snap-back structures that are resistant to bisulfite conversion. hydroxymethylation of dna is a newly discovered form of epigenetic regulation. it has been found recently in embryonic stem cells and in the brains of mammals. in mammals, tet proteins have been shown to be dioxygenases that convert mc to 5-hydroxymethylcytosine. honey bees have a tet ortholog, but hmc has not yet been reported. genome-wide wide mapping in mouse escs has revealed that hmc is enriched at the start sites of genes whose promoters bear histone lysine trimethylation and histone lysine trimethylation marks. in human escs, this dual mark is derived from separable subpopulations of self-renewing and lineage-biased escs within the heterogeneous unfractionated esc population. the most common chemical approach to study dna methylation is treating single-stranded dna with bisulfite, but bisulfite cannot distinguish mc from hmc because both base-pair as cytosine after bisulfite treatment. the modification mc mostly remains in this form after bisulfite treatment, whereas hmc is converted to cytosine methylene sulfonate after bisulfite treatment, which has the same hydrogen bond donor and acceptor configuration as cytosine for base pairing to guanine. in this paper, in addition to whole-genome shotgun bisulfite sequencing of honey bee head dna, we developed new biochemical and bioinformatics tools to analyze non-cg methylation and hydroxymethylation. we sharpened our analyses by comparing bees endemic to north america that are derived from a mixture of european subspecies of apis mellifera with bees derived from the african subspecies apis mellifera scutallata, introduced to south america in the last century. ehb and ahb are attractive candidates for comparative molecular analysis because they differ in many physiological and behavioral traits, including aggression. differences in methylation between ehb and ahb that are reported here may be related to these physiological and behavioral differences, and could motivate further studies, beyond the scope of the present paper. we show that non-cg methylation and hydroxymethylation are both enriched in introns, and cytosine modification at splice junctions might help regulate alternative splicing and gene expression.  we asked whether non-cg methylation and hydroxymethylation occurs in honey bees using unbiased approaches that do not filter out non-cg methylation. we compared ahb and ehb mainly to demonstrate the robustness of our experimental and bioinformatics approaches. to standardize our comparisons, we analyzed head methylomes only from ahb and ehb “guard” bees, a specialized group of individuals that patrol the hive entrance and respond first to a threat to their colony. most of the tissues in the bee head consist of brain. 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 from ahb and ehb heads and obtained over 20× coverage of both genomes. all four of the honey bee bs-seq studies conducted to date — three published studies, 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 ~ %, respectively in ahb. similarly, our re-analysis of the data from ref. identified ~8% cg methylation in exons and ~ % methylation in introns in ehb. only 15% to 21% of the cg methylation is symmetrically methylated in honey bees, which is lower than what is observed in mammals. our analysis methods, which we call bs-miner, 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. recently, a bs-seq analysis tool called bismark was developed that does not filter out non-cg methylation. the amount of chh and chg methylation identified in the ref. dataset by bismark was approximately the same as the amount of cg methylation, which is consistent with our analyses of our ahb and ehb datasets that show much more chh and chg methylation than previously reported. *the total number of cytosines represent both the uniquely mapped and the non-uniquely mapped cytosines. whole-genome shotgun bisulfite sequencing identifies non-cg modifications that are enriched in introns as in the previous studies, 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. is relatively high levels of non-cg modifications in bee heads. surprisingly, our methods detected about 5-fold more chh modifications than cg methylation in both ahb and ehb dna. as with the cg methylation, we also saw more than twice as many non-cg modifications in ehb than ahb heads. about % of the total number of chh sequences was modified in ehb and about % in ahb. in order to validate our finding of high levels of non-cg modifications, we reanalyzed previously published honey bee data from ref. with our methods and again found that the amount of non-cg modifications exceeded the amount of cg methylation. the differences in the amounts of non-cg modifications in our data compared with ref. might be caused by bisulfite treatment conditions or differences in strains of bees used in the two studies. we constructed the illumina libraries using the identical protocol as ref. 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. we also validated our data analysis pipelines by using two recent independently developed bisulfite methylation analysis program on the data from ref. and were again able to identify non-cg methylation sites. 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 and lowest in exons. consistent with this, our re-analysis of the ehb data from ref. identified ~ times more chh modifications in introns than exons. we also detected chg modifications at levels lower than both cg and chh. there was only about one-seventh as many chg modifications as chh modifications, and very few of the chg modifications were symmetrical. chg modifications were more uniform across the genome than cg methylation or chh modifications. coverage for individual chromosomes demonstrates that there were no significant biases toward any portion of the genome in the sequencing procedure. overall, ehb had almost × more modified cs than did ahb. we also observed ~3-4-fold more cg methylation than in the previous three studies: we detected methylated cgs in ehb, compared with - in the previous studies. this appears to be due to higher sensitivity of the analysis program used; as stated above, using our methods on data from ref. identified methylated cgs, which is even more than we identified in our data. the significance of ehb having × 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. we used highly sensitive anti-cms antibodies 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. for bodies, there was significantly more hmc in ehb than ahb. 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. the immunoprecipitated dna was then sequenced by next-generation dna sequencing. compared to previous findings that mc is found primarily in exons, we found that most of the hmc dna is present in introns, 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 restriction endonuclease, pvurts1i, which cuts near single hydroxymethylcytosine sites. the distances between the cleavage sites and the modified cytosine are fixed within a narrow range, with the majority being 11- nucleotides away in the top strand and 9- nucleotides away in the bottom strand. 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 cg ratio, is highly methylated, and is enriched in housekeeping genes, and a second has a high o/e ratio, is unmethylated, and is enriched in caste-specific and developmental genes. although non-cg sequences have a unimodal distribution, 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. genes with greater than 10% non-cg modifications, such as mc and hmc, in introns are primarily in low o/e cg genes, whereas genes with zero percent non-cg modifications in introns are in high o/e cg genes. 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, dnmt and the de novo enzyme, dnmt both of which are present in bees, and the bee tet protein. in contrast to dnmt which has almost exclusive specificity for cg sites, dnmt is responsible for most non-cg methylation in hesc; 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. this includes several homeobox transcription factors, such as distalless and aristaless. 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. as shown for distalless, 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.. b, distalless is a hox gene and a transcription factor with high amounts of chh modifications in the introns of both ahb and ehb. non-cg modifications might regulate alternative mrna splicing consistent with the idea that dna methylation might be involved in regulating mrna splicing, 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. we identified several hundred modified cytosines at splice donor and acceptor sites on the template strands. based on the above numbers, only ~ % of the ~ splice junctions were methylated in ahb and ~ % in ehb. 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. intriguingly, the go pathway “alternative splicing” was also significantly enriched in both ahb and ehb. for example, the honey bee gene that was most heavily methylated at splice junctions is gb the ortholog of drosophila dumpy, whose protein products are involved in cell adhesion in d. melanogaster; it has four methylated splice junctions in ahb and twelve methylated splice junctions in ehb at chh and chg sites. since dumpy has complex splicing programs in drosophila, 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’-tamcag-3’ on the non-coding dna strand at 3’ splice junctions. 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 + position on the template strand encoding splice acceptors and the - position on the template strand at splice donors, and very few of the chg sites in bees were symmetrically modified. 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. 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”. 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. 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, and there are extensive aggression-related differences in brain gene expression between the aggressive ahb and the less-aggressive ehb. 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. 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. consistent with previously published examples, cg methylation in dna encoding exons often correlated with that exon being skipped. for example, for gene gb we observed that a heavily methylated exon in ahb was skipped, whereas an adjacent heavily methylated exon in ehb was skipped. we also report for the first time that non-cg modifications, such as hmc, also showed a correlation with alternative splicing. for example, the gene gb 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. 1: figure s b, top blue lines indicated alternatively spliced exons in ahb and ehb, the red-and-blue stacked lines underneath the exons show pvu-seq analyses of ahb and ehb heads. note that hmc also correlates with alternative mrna splicing, but in manners that differ from mc. details of the rna-seq analyses are given in additional file 1: figure s  our findings underscore the diversity of dna methylation mechanisms that exist. non-cg modifications were only recently discovered in hesc and now we report them in the distantly related honey bee. we speculate that cytosine methylation at exons and splice junctions on the dna might affect the mrna splicing machinery. it is important to learn how these mechanisms work in relation to known regulators of splicing, such as histone acetylation, histone lysine methylation, histone lysine methylation, and histone lysine methylation. understanding how non-cg hydroxymethylation might affect alternative splicing is an exciting new area of research. our data are consistent with a model in which dna is methylated at cg sites by the maintenance dna methyltransferase, dnmt at exons, and at non-cg sites by the de novo dna methyltransferase, dnmt at introns. in contrast to dnmt mammalian dnmt is known to methylate both cg and non-cg sites in mammalian stem cells. we speculate that the honey bee tet enzyme primarily recognizes the non-cg sites in the introns, thereby enriching dna hydroxymethylation in the introns. we further speculate that the mrna splicing machinery, as well as the histone modification machinery, distinguishes exons and introns by somehow recognizing the patterns of cg methylation in exons and non-cg hydroxymethylation in introns. our identification of go categories related to protein phosphorylation that were enriched for genes with methylated splice junctions is consistent with a similar finding in a recent study of species-specific alternative exons. the authors present evidence that argues that alternative splicing is used to alter protein phosphorylation, which can alter protein stability, subcellular localization, activity, and other properties. further research is needed to determine the mechanism by which splice junction methylation and hydroxymethylation affect mrna splicing. 
additionally, using a new technique, pvu-seq, we further confirmed that hydroxymethylation is enriched in introns at non-cg sites. therefore, we used an unbiased analysis of bisulfite sequencing data combined with molecular and bioinformatics approaches to distinguish 5-methylcytosine from 5-hydroxymethylcytosine. previous whole-genome shotgun bisulfite sequencing experiments showed that dna cytosine methylation in the honey bee is almost exclusively at cg dinucleotides in exons. by doing this, we have performed the first whole genome analyses of dna modifications at non-cg sites in honey bees and correlated the effects of these dna modifications on gene expression and alternative mrna splicing. we confirmed, using unbiased analyses of whole-genome shotgun bisulfite sequencing data, with both new data and published data, the previous finding that cg dna methylation is enriched in exons in honey bees. cytosine hydroxymethylation at non-cg sites might have more functional significance than previously appreciated, and in honey bees these modifications might be related to the regulation of alternative mrna splicing by defining the locations of the introns. using antibodies against 5-hydroxmethylcytosine, we confirmed that dna hydroxymethylation at non-cg sites is enriched in introns. however, we also found evidence that cytosine methylation and hydroxymethylation at non-cg sites is enriched in introns.
 previous whole-genome shotgun bisulfite sequencing experiments showed that dna cytosine methylation in the honey bee is almost exclusively at cg dinucleotides in exons. however, the most commonly used method, bisulfite sequencing, cannot distinguish 5-methylcytosine from 5-hydroxymethylcytosine, an oxidized form of 5-methylcytosine that is catalyzed by the tet family of dioxygenases. furthermore, some analysis software programs under-represent non-cg dna methylation and hydryoxymethylation for a variety of reasons. therefore, we used an unbiased analysis of bisulfite sequencing data combined with molecular and bioinformatics approaches to distinguish 5-methylcytosine from 5-hydroxymethylcytosine. by doing this, we have performed the first whole genome analyses of dna modifications at non-cg sites in honey bees and correlated the effects of these dna modifications on gene expression and alternative mrna splicing.  we confirmed, using unbiased analyses of whole-genome shotgun bisulfite sequencing data, with both new data and published data, the previous finding that cg dna methylation is enriched in exons in honey bees. however, we also found evidence that cytosine methylation and hydroxymethylation at non-cg sites is enriched in introns. using antibodies against 5-hydroxmethylcytosine, we confirmed that dna hydroxymethylation at non-cg sites is enriched in introns. additionally, using a new technique, pvu-seq, we further confirmed that hydroxymethylation is enriched in introns at non-cg sites.  cytosine hydroxymethylation at non-cg sites might have more functional significance than previously appreciated, and in honey bees these modifications might be related to the regulation of alternative mrna splicing by defining the locations of the introns. honey beesdna methylationdna hydroxymethylationepigenetics 