Nondestructive enzymatic deamination enables single-molecule long-read amplicon sequencing for the determination of 5-methylcytosine and 5-hydroxymethylcytosine at single-base resolution

Abstract

The predominant methodology for DNA methylation analysis relies on the chemical deamination by sodium bisulfite of unmodified cytosine to uracil to permit the differential readout of methylated cytosines. Bisulfite treatment damages the DNA, leading to fragmentation and loss of long-range methylation information. To overcome this limitation of bisulfite-treated DNA, we applied a new enzymatic deamination approach, termed enzymatic methyl-seq (EM-seq), to long-range sequencing technologies. Our methodology, named long-read enzymatic modification sequencing (LR-EM-seq), preserves the integrity of DNA, allowing long-range methylation profiling of 5-methylcytosine (5mC) and 5-hydroxymethylcytosine (5hmC) over multikilobase length of genomic DNA. When applied to known differentially methylated regions (DMRs), LR-EM-seq achieves phasing of >5 kb, resulting in broader and better defined DMRs compared with that previously reported. This result showed the importance of phasing methylation for biologically relevant questions and the applicability of LR-EM-seq for long-range epigenetic analysis at single-molecule and single-nucleotide resolution.

© 2021 Sun et al.; Published by Cold Spring Harbor Laboratory Press.

Figures

Figure 1.

5mC and 5hmC detection by enzymatic deamination method. (A) Principle of the methodology: genomic DNA can either be treated with TET2 and BGT (left) to protect both 5mC and 5hmC or with BGT alone (right) to protect 5hmC. Subsequent deamination by APOBEC3A followed by PCR amplification allows the distinction between the unprotected substrate (read as T) from the protected cytosine derivatives (read as C). The TET2 and BGT treatment results in the distinction of 5mC and 5hmC from C, whereas BGT treatment results in the distinction of 5hmC from C and 5mC. (B) Deaminated cytosines from unmethylated lambda genome display no observable sequence preference by APOBEC(5mC) deamination method. (C) False-positive methylation calling rate (nonconversion error rate) of each cytosine dinucleotides sequence context (CpA [CA], CpC [CC], CpG [CG], and CpT [CT]) estimated from the unmethylated lambda genome for the enzymatic deamination method (APOBEC(5mC), two WGBS performed in this study, i.e., BS kit 1 and BS kit 2, and six published WGBS experiments sampled from the ENCODE Project) (Supplemental Table S3). (D) Deamination of 5mC in the fully methylated XP12 genome results in no observable sequence preference by APOBEC(5hmC) enzymatic deamination method. (E) Distribution patterns of 5mCpG (blue) and 5hmCpG (red: 50 ng library; pink: 1 ng library) at various protein/DNA interaction sites. The absolute (smooth lines) and normalized (dotted lines) 5hmC and 5mC levels in the CpG context are depicted around TET1, RNA polymerase II, and CTCF binding sites, as well as at active transcription chromatin mark (H3K4me3), repressive chromatin mark (H3K27me3), active enhancer mark (H3K27ac), and general enhancer (H3K4me1 in the absence of H3K4me3) regions. Unbound sites that are randomly sampled from the reference genome server as a control. (F) Pearson's correlation between 5hmC measured using sequencing of enzymatically deaminated DNA (x-axis) versus LC-MS (y-axis) for various genomic DNA. There are two technical replicates of the APOBEC(5hmC) sequencing method for each sample. 5hmC levels are presented as 1000 percentage, and both axes use the log scale.

Figure 2.

Enzymatic deamination preserves the integrity of the DNA. (A) qPCR results show the quantities of undamaged amplifiable DNA templates of different sizes after the enzymatic deamination (green) and bisulfite treatments (orange and blue). All quantifications are normalized to the values obtained for the enzymatic deamination experiments. (B) Agilent 2100 Bioanalyzer trace on RNA 6000 pico chip comparing equal amounts of mouse E14 genomic DNA sheared to an average of 15 kb and treated with sodium bisulfite (green), APOBEC(5hmC) (red), or APOBEC(5mC) (blue) over the control ssDNA (cyan). Bisulfite treatment fragmented the DNA to an average of 800 bp, whereas enzymatically treated DNA shows no notable size differences compared with control DNA. (C) Agarose gel images of end-point PCR of six amplicons ranging from 388–4229 bp illustrating upper amplicon size limit for sodium-bisulfite-, APOBEC(5mC)-, or APOBEC(5hmC)-treated E14 genomic DNA. (D) The 731-bp amplicons from the E14 genomic DNA shown in C were cloned and sequenced, and the methylation status was determined by bisulfite treatment (left), the enzymatic deamination method for 5mC (center), and the enzymatic deamination method for 5hmC (right) (Supplemental Data S1). Open and closed circles indicate unmethylated and methylated CpG sites, respectively.

Figure 3.

5mC and 5hmC phasing using long-read sequencing. (A) Scatter plots and Pearson's correlations of calculated methylation (top) and hydroxymethylation (bottom) levels of all CpG sites within the 5378-bp region from the mouse E14 genome between the three sequencing platforms: PacBio, Nanopore, and Illumina. (B) Dot plots showing methylation (left) and hydroxymethylation (right) levels of individual CpG sites within the 5378-bp region calculated by the LR-EM-seq method using three major sequencing platforms: Illumina (red), Nanopore (green), and PacBio (blue). The fitted lines are drawn using the LOESS method. (C) Single-base single-molecule cytosine modification maps of the 5378-bp region generated by the LR-EM-seq method coupled with PacBio SMRT sequencing (top) and Nanopore sequencing (bottom). Methylated (left) and hydroxymethylated CpG sites are depicted by red dots, and unmodified CpG sites are depicted by beige dots.

Figure 4.

Phasing of 5mC and 5hmC by LR-EM-seq. (A) Single-base, single-molecule CpG methylation (middle) and hydroxymethylation (bottom) profile of a 4.6-kb region of the imprinted Inpp5f_v2 gene locus (top) in the mouse brain. Red dots represent modified sites, and beige dots represent unmodified sites. This region overlaps with the promoter for the Inpp5f_v2 gene and contains a previously reported DMR (orange box). The shaded area in the dot plots corresponds to the known DMR. (B) Correlation matrix of CpG modification state: (top) 5mC; (bottom) 5hmC. Each location on the matrix represents the correlation of any two CpG sites across the amplicon and the correlation strength is depicted by color: red indicates correlation=1; blue, correlation=−1; white, no correlation. The known DMR is indicated by a black outline.

Figure 5.

Phasing of 5mC with heterozygous variants using LR-EM-seq. (A) Phasing of 5mC with SNP of a 3.1-kb region in the imprinted Inpp5f_v2 gene promoter of the mouse cortex brain from a F1 offspring of a cross between two inbred mouse strains (129X1/SvJ male and Cast/EiJ female). Methylation state of individual CpG sites at the single-molecule level is denoted by either a beige dot (unmodified) or a red dot (methylated). The heterozygous SNP near the 5′ end of the region was either highlighted in red for paternal allele (A) or blue for maternal allele (G). The orange boxes denote previously identified DMRs. Our result not only confirmed the existence of the imprinted DMR but also revealed much extended boundaries of the imprinted DMR. (B) Phasing of 5mC and SNP in the imprinted promoter of the Gnas gene in the mouse cortex from a cross between the inbred mouse strains 129X1/SvJ (male) and Cast/EiJ (female). Methylation state of individual CpG sites at single-molecule level is denoted by either a beige dot (unmodified) or a red dot (methylated). The heterozygous SNP was highlighted in red for paternal allele (A) and blue for maternal allele (G). The orange box denotes a previously identified DMR. Our result confirmed the existence of the imprinted DMR and further extended this DMR in both directions particularly into the CpG island.