5-Hydroxymethylcytosine is associated with enhancers and gene bodies in human embryonic stem cells

Abstract

Background: 5-Hydroxymethylcytosine (5hmC) was recently found to be abundantly present in certain cell types, including embryonic stem cells. There is growing evidence that TET proteins, which convert 5-methylcytosine (5mC) to 5hmC, play important biological roles. To further understand the function of 5hmC, an analysis of the genome-wide localization of this mark is required.

Results: Here, we have generated a genome-wide map of 5hmC in human embryonic stem cells by hmeDIP-seq, in which hydroxymethyl-DNA immunoprecipitation is followed by massively parallel sequencing. We found that 5hmC is enriched in enhancers as well as in gene bodies, suggesting a potential role for 5hmC in gene regulation. Consistent with localization of 5hmC at enhancers, 5hmC was significantly enriched in histone modifications associated with enhancers, such as H3K4me1 and H3K27ac. 5hmC was also enriched in other protein-DNA interaction sites, such as OCT4 and NANOG binding sites. Furthermore, we found that 5hmC regions tend to have an excess of G over C on one strand of DNA.

Conclusions: Our findings suggest that 5hmC may be targeted to certain genomic regions based both on gene expression and sequence composition.

Figure 1

High resolution map of hydroxymethylcytosine in human embryonic stem cells. (a) Genome-browser view of hmeDIP-seq data. Two hmeDIP-seq datasets along with input DNA and 'no antibody' controls are shown. Each track is represented as normalized density of reads (reads/bp/million uniquely mapping reads). (b) Gene density (genes/bp) and 5hmC peak density (peaks/bp) in chromosome 3. (c) 5hmC over genes with different expression levels. 5hmC peak density was plotted over RefSeq genes in 300-bp bins. Plots were smoothed by taking the moving average over ± 2 bins. Published RNA sequencing data [9] were used to rank the genes. TSS, transcription start site; TTS, transcription termination site.

Figure 2

Co-localization of 5hmC with enhancers. (a) 5hmC over enhancers in hESCs. 5hmC peak density was plotted over predicted hESC enhancers [9,21] in 100-bp windows. (b) Verfication of 5hmC at enhancers by measuring 5hmC levels at CCGG sites. Four 5hmC peaks at enhancers and two control regions were tested. MspI cannot cut glucosylated-5hmC but is able to cut 5hmC; therefore, copy numbers of MspI + beta-glucosyltransferase (BGT) represent 5hmC levels. On the other hand, HpaII can only cut unmodified DNA; therefore, copy numbers represent 5hmC + 5mC levels. Background signal was subtracted from the copy number of each sample and then normalized to the undigested glucosylated sample. Genomic locations of tested cytosines are indicated. Error bars represent the standard deviation. (c) Histone modifications over 5hmC peaks. The overlap of 5hmC peaks with previously defined H3K4me1- and H3K27ac-enriched regions [12] was calculated. Random regions with the same number and size distribution as 5hmC peaks were generated and overlap with histone modifications was calculated 100 times. Error bars represent standard deviation. (d) Histone modifications over 5hmC-marked enhancers. Predicted enhancers that overlapped with 5hmC peaks were selected. The overlap of these 5hmC-marked enhancers, as well as all predicted enhancers, with previously defined H3K4me1 and H3K27ac enriched regions [12] was calculated. Random regions with the same number and size distribution as the enhancers were generated and overlap with histone modifications was calculated 100 times. Error bars represent standard deviation. (e) hESC-specific expressed genes significantly overlap with 5hmC peaks. hESC-specific genes were defined as genes that were expressed in hESCs (reads per kilobase of exon model per million mapped reads (RPKM) ≥ 0.5) and silent in IMR90 cells (RPKM = 0) using published RNA-seq data [9]. *P < 0.0001.

Figure 3

Co-localization of 5hmC with transcription factor binding sites. (a) 5hmC regions over CTCF binding sites in hESCs. 5hmC peak density was plotted over the ± 2 kb relative to the center of defined binding sites in 100-bp windows. (b) Over NANOG binding sites in hESCs. (c) Over OCT4 binding sites in hESCs.

Figure 4

5hmC regions are GC-skewed. (a) Base composition of the Watson strand over the 5' and 3' boundaries of 5hmC regions. GC-skew = absolute value (G - C)/(G + C). (b) Base composition of the Watson strand over the centers of 5hmC regions.