Single-base resolution analysis of active DNA demethylation using methylase-assisted bisulfite sequencing

Active DNA demethylation in mammals involves TET-mediated iterative oxidation of 5-methylcytosine (5mC)/5-hydroxymethylcytosine (5hmC) and subsequent excision repair of highly oxidized cytosine bases 5-formylcytosine (5fC)/5-carboxylcytosine (5caC) by Thymine DNA glycosylase (TDG). However, quantitative and high-resolution analysis of active DNA demethylation activity remains challenging. Here we describe M.SssI methylase-assisted bisulfite sequencing (MAB-seq), a method that directly maps 5fC/5caC at single-base resolution. Genome-wide MAB-seq allows systematic identification of 5fC/5caC in Tdg-depleted embryonic stem cells, thereby generating a base-resolution map of active DNA demethylome. A comparison of 5fC/5caC and 5hmC distribution maps indicates that catalytic processivity of TET enzymes correlates with local chromatin accessibility. MAB-seq also reveals strong strand asymmetry of active demethylation within palindromic CpGs. Integrating MAB-seq with other base-resolution mapping methods enables quantitative measurement of cytosine modification states at key transitioning steps of active demethylation pathway, and reveals a regulatory role of 5fC/5caC excision repair in active DNA demethylation cascade.

development and cell-type specific gene expression. In mammals, DNA demethylation (conversion of 5mC to unmodified C) can be achieved either passively through successive rounds of DNA replication in the absence of functional DNA methylation maintenance machinery or actively by Ten-eleven translocation (TET) family of 5mC-modifying enzymes 5,6 . TET proteins are Fe 2+ and 2-oxoglutarate-dependent dioxygenases capable of successively oxidizing 5mC to 5hmC, 5fC and 5caC (Fig. 1a) [7][8][9][10] . Oxidative modification of 5mC by TET enzymes promotes DNA demethylation by either replication-dependent dilution of oxidized cytosines 11 or TDG-mediated 5fC/5caC excision followed by baseexcision repair (BER) 10,12,13 . The active demethylation pathway involving generation and excision repair of 5fC/5caC is of particular interest as it may take place in a wide range of somatic cell-types including post-mitotic cells. Interestingly, only TDG, but not other members of the uracil DNA glycosylase (UDG) superfamily, possesses robust 5fC/5caC excision activity and is indispensable for embryonic development 14,15 , implicating the TET/TDG-dependent DNA demethylation pathway in regulating tissue-specific gene expression and development.
The observation that 5hmC can accumulate to a relatively high level in diverse cell-types, particularly in adult neurons 16,17 , raises the possibility that TET proteins tend to stall at 5hmC and that further oxidation of 5hmC to 5fC/5caC is a rate-limiting step. Therefore, identifying cytosines that are committed to active DNA demethylation requires methods that permit quantitative measurement of TDG-mediated excision of 5fC/5caC at high resolution. In addition, because both 5fC and 5caC are selectively recognized and excised by TDG, determining the strand-specific preference of active DNA demethylation activity requires the ability to simultaneously map both 5fC and 5caC. Recent studies using either 5fC/5caCspecific antibodies or chemical tagging of 5fC have shown that genomic regions containing cytosines undergoing TET/TDG-dependent active demethylation can be identified by analyzing ectopic 5fC/5caC accumulation in Tdg-depleted cells 18,19 . However, 5fC and 5caC maps generated by affinity-enrichment methods are of limited resolution (several hundred base-pairs), only represent relative enrichment, and lack strand distribution information. To address these limitations, we have developed a method named MAB-seq, which allows simultaneous and quantitative mapping of both 5fC and 5caC at baseresolution. Furthermore, our genome-wide MAB-seq analysis of mouse ESCs has provided new insights into catalytic processivity, strand asymmetry and feedback regulation of the TET/TDG-mediated active DNA demethylation pathway. and sequenced as C (Fig. 1a). In MAB-seq, genomic DNA is first treated with the bacterial DNA CpG methyltransferase M.SssI, an enzyme that is originally isolated from Spiroplasma sp. strain MQ1 and is known to efficiently methylate cytosines within CpG dinucleotides 22 . Bisulfite conversion of M. SssI-treated DNA may therefore only deaminate 5fC and 5caC; originally unmodified C within CpGs is protected as 5mC. Subsequent sequencing would reveal 5fC and 5caC as T, whereas C/5mC/5hmC as C ( Supplementary Fig. 2a). Notably, MAB-seq is unable to distinguish 5fC/5caC from unmodified C within non-CpG context due to the poor activity of M. SssI towards C within non-CpG context. This limitation does not affect the application of this technique for two reasons: first, 5hmC is almost exclusively found in the CpG context (>99% in CpGs), even in mouse ESCs and neurons where non-CpG methylation is prevalent 17,23 ; second, recent structural and biochemical analyses indicate that TET proteins have a strong preference for oxidizing 5mC in CpG sites than in non-CpG context 24,25 . Thus, oxidative modification of 5mC by TET proteins occurs predominantly in the CpG context, and MAB-seq may provide a quantitative measurement of the abundance of 5fC/5caC within CpG dyads.
Successful detection of 5fC/5caC using MAB-seq requires complete conversion of C to 5mC by M. SssI as well as efficient bisulfite conversion of 5fC/5caC. We first optimized the reaction conditions and achieved nearly complete (99.2%) conversion of unmodified CpGs to 5mCpGs by M. SssI methylase measured by Sanger sequencing (Supplementary Fig. 1b). Next, we performed high-throughput BS-seq analysis of a synthetic double-stranded DNA (dsDNA) containing CpGs with specific cytosine modifications (5hmC/5fC/5caC). Consistent with previous reports 21 , our analysis showed that 5fC (84.7%) and 5caC (99.5%), but not 5hmC (3.3%), are efficiently deaminated by bisulfite treatment and read as T ( Supplementary Fig. 2b). In addition to unmodified CpGs (C:C), asymmetrically-modified CpGs (5mC/5hmC/5fC/5caC:C) may be present at low levels in the genome 23 . We thus tested MAB-seq in analyzing asymmetrically modified dsDNA (5hmC/5fC/5caC:C). This analysis demonstrated that M. SssI methylase is capable of efficiently methylating unmodified C in hemi-modified CpG dyads ( Supplementary Fig. 2c), validating the capability of MAB-seq in mapping asymmetrically modified 5fC/5caC in a strand-specific manner.
We next performed BS-seq and MAB-seq analysis of unmethylated lambda phage genome (6,224 CpGs within 48,502 bp) using Illumina high-throughput sequencing and sequenced to an average depth of 239x and 305x per cytosine, respectively. In BS-seq, a nearly complete C-to-T conversion within CpG sites was observed (99.9 +/− 0.06%, n=3 experiments) contrasted to a low conversion rate in MAB-seq (2.04 +/− 0.14%, n=9) ( Fig. 1b and Supplementary Fig. 3a). To test whether unprotected CpGs in MAB-seq experiments exhibit random distribution, we analyzed the sequences immediately flanking 67 CpGs (mean methylation: 94.1%) that are not efficiently methylated by M. SssI ( Supplementary Fig. 3bc). We found that these 67 CpGs are not associated with any specific sequences ( Supplementary Fig. 3d), suggesting that M. SssI has minimal sequence preference for catalyzing CpG methylation reactions. Consistent with previous findings 22 , we found that M. SssI only methylates a small fraction of cytosines (1.3%) within non-CpG context ( Fig.  1b and Supplementary Fig. 3a).
To test MAB-seq in analyzing mammalian genomic DNA, we applied this method to examine four 5fC/5caC-enriched loci (Tbx5, Vps26a, Ace, and Slc2a12) that were previously identified in Tdg-depleted mouse ESCs by affinity-enrichment methods 18,19 (Fig. 1c, Supplementary Fig. 4 and 5a). Using mouse ESCs deficient for both Tet1 and Tet2 (largely absent of 5fC/5caC) as a negative control to correct background signals, locus-specific MAB-seq analysis not only located individual CpGs associated with significant levels of 5fC/5caC ( Fig. 1d and Supplementary Fig. 6; P<0.05, Fisher's exact test), but also revealed the absolute level of 5fC/5caC at each CpG ( Fig. 1c and Supplementary Fig. 4). In contrast, active promoter of Rest, where 5fC/5caC is undetectable by affinity methods, is largely absent of MAB-seq signals (left panel in Supplementary Fig. 5a). This analysis also confirmed the positive effect of vitamin C (VC) on the catalytic activity of TET enzymes 26,27 , as VC-treated Tdg-depleted mouse ESCs (shTdg+VC) display even higher level of 5fC/5caC compared to shTdg cells ( Fig. 1c and Supplementary Fig. 4, 6). Notably, a small number of 5fC/5caC-modified CpGs were also identified in the shCtrl and shCtrl+VC samples ( Fig. 1d and Supplementary Fig. 6). Although the levels of 5fC/5caC in shCtrl are much lower than those in shTdg, the results suggest that MAB-seq is highly sensitive and may identify rare 5fC/5caC modifications in wild-type cells. Similar 5fC/5caC distribution patterns were detected for biologically independent replicates, demonstrating the reproducibility and robustness of MAB-seq ( Supplementary Fig. 5b).

Genome-wide MAB-seq analysis of mouse ESCs
Having validated MAB-seq using locus-specific analysis, we next applied MAB-seq to identify cytosines undergoing active DNA demethylation at a genome-scale. Affinityenrichment-based studies of 5fC/5caC distributions in mouse ESCs have shown that a large cohort of 5fC and 5caC peaks overlap with genomic regions enriched for H3K4me1 18, 19 , a histone mark broadly associated with active/poised enhancers, flanking regions of active/ poised gene promoters, and intragenic regions of actively transcribed genes 28 . Because H3K4me1-marked regions enrich for 5fC/5caC signals, we focused our initial genome-scale MAB-seq analysis on H3K4me1-enriched genomic regions captured through chromatin immunoprecipitation (termed H3K4me1-MAB-seq) ( Supplementary Fig. 7a). To establish the baseline of false positive 5fC/5caC signals, we performed H3K4me1-MAB-seq analysis of mouse ESCs deficient for all three Dnmt enzymes (Dnmt1/3a/3b −/− ) or Tet1/2 proteins (Tet1/2 −/− ). We found that median false positive signals in these mutant ESCs, in which both 5fC and 5caC are virtually absent, are very close to the error rate of M. SssI observed in methylating lambda DNA (dashed line in Fig. 2b). In comparison to control knockdown (shCtrl), Dnmt1/3a/3b −/− , Tet1/2 −/− ESCs, 5fC/5caC is present at significantly higher levels in Tdg knockdown (shTdg) mouse ESCs (P<2.2×10 −16 , Wilcoxon rank sum test), and is further increased in VC-treated cells (shTdg+VC in Fig. 2b).
Given that 5fC/5caC is absent in Dnmt/Tet mutant (denoted as 'Neg Ctrl' thereafter) mouse ESCs, MAB-seq signals in these mutant cells may provide an empirical estimate of false discovery rate (FDR). Because the probability that a CpG can be confidently identified as 5fC/5caC-modified is governed by the sequencing depth and abundance of the modification at the cytosine, we modeled the largely stochastic event of M.SssI failure in CpG methylation with a binomial distribution [N as the depth of sequencing at the cytosine and p (2.04%) as the error rate of M. SssI]. Application of this statistical strategy to H3K4me1-MAB-seq datasets identified a total of 127,576 (7.6% out of 1,670,036 CpG dyads with N≥10) 5fC/5caC-modified CpG dyads in VC-treated, Tdg-depleted mouse ESCs (shTdg +VC) with an empirical FDR of 5%. Using a less stringent cutoff (N≥5, FDR<10%), we have identified 267,325 (8.0% out of 3,326,034 CpG dyads) 5fC/5caC-modified CpGs. There are 17.6 times as many 5fC/5caC-modified CpGs overlapping with affinity-identified regions as expected by chance (Fig. 2c), demonstrating the effectiveness of the empirical FDR-based statistical filter.
We next analyzed genomic DNA from mouse ESCs by whole genome (WG)-MAB-seq. Since VC is present at a relatively high level during early development in vivo 29 , baseresolution 5fC/5caC map in VC-treated cells may represent a more physiologically relevant profile of TET/TDG-dependent active demethylation activity. We therefore focused our WG-MAB-seq analysis on shTdg+VC mouse ESCs and sequenced the sample to an average depth of 28.4x per CpG dyad (covering 95.2% of all CpG dyads). Locus-specific analysis of non-enriched and H3K4me1-enriched DNA suggests that 5fC/5caC levels at specific CpGs within selected loci are largely comparable ( Supplementary Fig. 7b, c). Further comparative analysis of both H3K4me1-MAB-seq and WG-MAB-seq experiments indicates that MABseq signals (T/C+T) are highly similar in both experiments at a wide range of promoterproximal and distal genomic elements ( Supplementary Fig. 8a), supporting the validity of using the statistical filter established for H3K4me1-MAB-seq to analyze WG-MAB-seq dataset. Using the empirical P value cutoff established for H3K4me1-MAB-seq datasets with comparable sequencing depth, we identified a total of 675,325 5fC/5caC-modified CpGs (out of 24,872,637 CpGs [N≥10]) in shTdg+VC mouse ESCs through whole genome MAB-seq analysis (with an empirical FDR of 5%). Identified 5fC/5caC-modified CpGs in Tdg-depleted cells correlate well with peaks of 5fC and 5caC enrichment identified by the antibody-based DNA immunoprecipitation (DIP) approach (region 1 in Fig. 2a). Compared to random controls, 5fC/5caC-CpGs significantly overlap with 5fC/5caC-enriched peaks (7.9 times as many as expected by chance). Furthermore, 83.6% of affinity-enrichment identified regions (n=50,923 out of 60,912 covered by WG-MAB-seq) overlap with at least one 5fC/5caC-modified CpGs. In contrast, 80.7% of 5fC/5caC-modified CpGs are not recovered by 5fC/5caC DIP-seq (region 2 in Fig. 2a), suggesting that MAB-seq has a markedly increased sensitivity.

Genomic distribution of active DNA demethylation activity
Previous studies using affinity-enrichment based methods have shown that 5fC/5caC are enriched at poised/active enhancers, Polycomb group protein (PcG) repressed promoters, and gene bodies 18,19 . However, the true abundance of 5fC/5caC cannot be determined from affinity-based approaches, thus precluding quantitative analysis of 5fC/5caC-excisiondependent DNA demethylation events at these gene regulatory elements. In shTdg+VC mouse ESCs, we found that a significant fraction of 5fC/5caC (20.1%) reside in distal regulatory regions (Fig. 2d). Analysis of both raw MAB-seq signals and relative enrichment of called 5fC/5caC at each class of genomic features indicates that 5fC/5caC are more enriched at DNase I hypersensitive sites (observed/random [obs/rand]=5.2), predicted enhancers (obs/rand=4.6), and CTCF-biding sites (obs/rand=2.4) relative to promoter-proximal regions and other genic regions ( Fig. 2e and Supplementary Fig. 8c-d). In support of this observation, regions enriched for pluripotency-related transcription factors (Oct4, Nanog, and Sox2) as well as active enhancers (marked by H3K27ac and/or H3K4me1) are also more enriched with 5fC/5caC than with other genomic elements ( Fig. 2f and Supplementary Fig. 8b). Consistent with previous reports 18, 19 , we found that 5fC/5caC is more enriched at mouse ESC-specific distal enhancers than at other tissue-specific enhancers ( Supplementary Fig. 8e-f), supporting the notion that CpGs within or surrounding active cell-type-specific enhancers tend to undergo 5fC/5caC-excision-dependent active DNA demethylation.
Given that 5fC/5caC are preferentially enriched at active/poised gene regulatory regions where chromatin is generally more accessible, we reasoned that local chromatin structure may influence the occupancy level and/or processivity (the ability to further oxidize 5hmC to 5fC/5caC) of TET enzymes. We first analyzed DNase I hypersensitivity (measured by DNase-seq 30 ) at genomic regions (±50 bp) immediately flanking 5fC/5caC-alone, 5hmCalone and dual modified CpGs. This analysis reveals that 5fC/5caC-alone and dual modified CpGs are associated with significantly higher level of DNase I hypersensitivity signals than 5hmC-alone CpGs ( Fig. 3b and Supplementary Fig. 9c-d). Further analysis indicates that as compared to 5hmC-alone CpGs, 5fC/5caC-alone and dual modified CpGs are associated with higher levels of Benzonase (BNase) sensitivity signals 31 , histone variants (H2A.Z and H3.3) known to destabilize nucleosome structure 31, 32 , TET1 occupancy 33 and pluripotency-related transcription factors (TFs; Oct4, Nanog, and Sox2) 34 ( Fig. 3c and Supplementary Fig. 9d). In contrast, comparable levels of general histone H3 were detected at distinct groups of CpGs (Fig. 3c). These results suggest that TET proteins by default tend to stall at the 5hmC step at most CpGs, but exhibit higher processivity at CpGs associated with a more accessible chromatin state.

Strand asymmetry of 5fC/5caC in palindromic CpGs
While cytosine methylation in palindromic CpG dyads is generally symmetric and exhibits very high heritability upon DNA replication 2 , previous studies suggest that >80% of steadystate 5hmC in human and mouse ESCs are asymmetrically modified 23 . This prompted us to examine whether 5fC/5caC-modified CpGs also show strand biases (Fig. 4a). Focusing on CpG dyads with both strands sufficiently covered by WG-MAB-seq (n=9,261,306 CpG dyads and depth≥10 on both Watson and Crick strands), we found that only 4.97% of called CpG dyads (22,590 out of 454,400 called CpG dyads, FDR=5%) are symmetrically modified with 5fC/5caC. Further analysis focusing on CpG dyads with both strands covered by higher sequencing depth (≥20) reached a similar conclusion (7.92% symmetrically modified with 5fC/5caC). However, because the abundance of 5fC/5caC is low at any given CpG dyad, it is still possible that sequencing depth in WG-MAB-seq might not be sufficient to identify all 5fC/5caC, leading to an underestimation of symmetrically-modified CpGs. To address this issue, we pooled MAB-seq signals of all called 5fC/5caC-modified CpGs and compared the called strand to the opposite strand. The average abundance of 5fC/5caC at the called CpGs is 26.9%, while that of the opposite CpGs is only 8.5% (Fig. 4a). To further confirm the observed strand asymmetry of 5fC/5caC-modified CpG dyads, we analyzed the difference in 5fC/5caC levels of called and opposite strands for each called CpG dyads (FDR=5%, depth≥10). For asymmetrically-modified CpGs, the mean difference in 5fC/5caC levels between called and opposite cytosines is on average 18.6% at each CpG dyad (or 14.4% for CpG dyads with depth≥20) (Fig. 4b). Visualizing the absolute levels of 5fC/5caC on both called and opposite strands in a two-dimensional histogram plot showed a clear shift toward the called cytosines (Fig. 4c). Notably, similar analysis of H3K4me1-MAB-seq datasets also supports the observed strand asymmetry of 5fC/5caC-modified CpG dyads ( Supplementary Fig. 6d-f). Further analysis showed no strand bias for M.SssI in methylating lambda phage genome (97.97% for Watson strand versus 97.90% for Crick strand) or genomic DNA of Dnmt1/3a/3b −/− cells (97.82% for Watson strand and 97.80% for Crick strand). Consistent with genome-scale MAB-seq analysis, strand-specific analysis of Tbx5 and Vps26a loci also revealed that the majority of 5fC/5caC-modified CpG dyads exhibit strand asymmetry [asymmetric (blue)+partially asymmetric (yellow): 86.7% in shTdg and 60% in shTdg+VC] (Fig. 4d-e). Locus-specific MAB-seq analysis of biological replicates suggests that the strand-specific distribution of 5fC/5caC detected at each CpG sites is not entirely stochastic ( Supplementary Fig. 5b), suggesting that strand-specific 5fC/5caC generation and excision at individual CpGs is a regulated process.

Base-resolution mapping of 5fC and 5caC separately
Base-resolution mapping of 5fC has recently been demonstrated using subtraction-based, chemical modification-assisted BS-seq methods such as fCAB-seq 19 and redBS-seq 21 ( Supplementary Fig. 1a). However, relatively low protection rate of 5caC deamination (50-60%) reported in a similar subtraction-based 5caC mapping approach (caCAB-seq) suggests that further optimization is required 20 . To achieve 5caC mapping at base-resolution, we combined MAB-seq with the 5fC mapping method, redBS-seq 21 . After validating redBS-seq using synthetic oligonucleotides ( Supplementary Fig. 10a), we performed locus-specific redBS-seq and BS-seq to quantify 5fC at 73 CpG sites from four loci enriched for 5fC and/or 5caC (Tbx5, Vps26a, Slc2a12 and Ace). The abundance and position of 5caC at these sites can then be determined by subtracting 5fC signals (derived from the difference between redBS-seq and BS-seq) from the levels of 5fC+5caC measured by MAB-seq ( Fig. 5b and Supplementary Fig. 10c-d). Using this integrative approach to map 5fC and 5caC separately, we observed that 5fC and 5caC display largely distinct distribution patterns at individual CpGs within these loci. For instance, within exon2 of Tbx5, a region enriched for both 5fC and 5caC, some CpG sites are only modified with 5fC (CpG #1 at Tbx5 in Fig. 5b), while some others are 5caC-only (CpG #2 at Tbx5 in Fig. 5b). More strikingly, as exemplified by CpG #3 at Tbx5 in Fig. 5b, CpG sites associated with both 5fC and 5caC may exhibit non-overlapping, strand-specific 5fC/5caC distribution. This conclusion is further supported by analysis of an active enhancer within the Vps26a gene ( Supplementary  Fig. 10d). Such diverse distribution patterns of 5fC and 5caC indicate distinct processivity of TETs and/or substrate preference of TDG (5fC versus 5caC) at individual CpGs. Notably, most 5fC-modified CpGs (20 out of 24) identified by redBS-seq are also detected by MABseq as 5fC/5caC-modified ( Supplementary Fig. 10c). The few inconsistent sites between MAB-seq and redBS-seq probably arise from the different principles that the two methods are based upon.
Base-resolution mapping 5caC through integrating MAB-seq and redBS-seq requires two rounds of subtractions, representing a technical challenge for genome-scale analysis. Thus, we explored a subtraction-independent 5caC mapping strategy by taking advantage of the fact that 5fC can be selectively reduced by sodium borohydride (NaBH 4 ) to 5hmC 35 . In this modified version of MAB-seq (termed caMAB-seq), M.SssI-treated DNA was further incubated with NaBH 4 so that only 5caC is read as T in bisulfite sequencing (Fig. 5a). Analyses of defined sequences (lambda DNA or synthetic oligonucleotides) or 5caCenriched genomic loci (e.g. Vps26a) demonstrated the validity of caMAB-seq in direct baseresolution mapping of 5caC ( Fig. 5c and Supplementary Fig. 10e).

Integrative analysis of unmodified C, 5mC/5hmC and 5fC/5caC
Standard BS-seq cannot distinguish unmodified C from 5fC/5caC. When combined with BS-seq, MAB-seq is able to quantify the true abundance of unmodified CpG at baseresolution (Fig. 1a), providing a unique opportunity for exploring the potential role of 5fC/ 5caC excision in regulating upstream (5mC+5hmC) and downstream (unmodified C) steps in the TET/TDG-dependent active DNA demethylation pathway. By performing locusspecific BS-seq and MAB-seq through Illumina deep sequencing, we quantified the levels of unmodified C, 5mC/5hmC and 5fC/5caC at 73 CpGs within four 5fC/5caC-enriched regions. BS-seq analysis showed that Tdg-depletion induces significant decrease in 5mC/ 5hmC within these loci (shTdg versus shCtrl or shTdg+VC versus shCtrl+VC in Fig. 6a). Further analysis shows that significant decrease in 5mC/5hmC was detected at multiple CpGs within each locus (Fig. 6a-b and Supplementary Fig. 11a). In contrast to 5mC/5hmC, integrative analysis of BS-seq and MAB-seq results revealed that unmodified C is significantly increased in response to TDG depletion (Fig. 6c-d and Supplementary Fig.  11b).
Tdg-depletion induced changes in 5mC/5hmC and unmodified C appear to be more pronounced at Tbx5, Ace and Slc2a12 loci when compared to Vps26a locus (an active enhancer), suggesting that active DNA demethylation dynamics at these transcriptionally repressed/poised regions (previously determined by ChIP-seq) is modulated by the generation and/or excision of 5fC/5caC. To test this possibility at a genome-scale, we captured genomic regions marked by H3K27me3, a repressive histone modification associated with transcriptionally poised gene promoters of lineage-specific TFs (e.g. Tbx5), and subjected H3K27me3-enriched DNA to genome-scale MAB-seq and BS-seq analyses. H3K27me3-MAB-seq identified 1,231 genomic intervals enriched for 5fC/5caC signals in shTdg+VC cells, and H3K27me3-BS-seq showed that there was a significant decrease in 5mC/5hmC levels within these 5fC/5caC-enriched regions (P= 3.18×10 −22 , Wilcoxon) ("MAB-seq enriched" in Fig. 6e). For 1,231 randomly sampled regions where 5fC/5caC signals are largely absent, the change in 5mC/5hmC is much less pronounced ("random" in Fig. 6e). Consistent with results of locus-specific analysis, genome-scale analysis supports the notion that inhibition of TDG-mediated 5fC/5caC or 5fC/5caC accumulation itself leads to dysregulation of upstream steps of the cytosine-modifying cascade.
Notably, VC treatment alone can substantially reduce the level of 5mC/5hmC (shCtrl+VC versus shCtrl or shTdg+VC versus shTdg in Fig. 6a and upper panels in Fig. 6f), possibly due to its positive effect on stimulating TET catalytic activity (oxidizing 5mC/5hmC). In addition, VC treatment can induce a significant increase in the level of unmodified C (shCtrl +VC versus shCtrl or shTdg+VC versus shTdg in Fig. 6c). However, VC-induced increase in unmodified C is more pronounced in shCtrl (purple in lower panel of Fig. 6f) than in shTdg (yellow in lower panel of Fig. 6f), indicating that accumulation of unmodified C by VC treatment requires TDG/BER-mediated 5fC/5caC excision step. Together, these results suggest that interfering 5fC/5caC excision step by depleting TDG proteins may result in dysregulation of both upstream (5mC/5hmC generation or oxidization) and downstream (generation and/or methylation of unmodified C) steps of the DNMT/TET/TDG-dependent cytosine-modifying cascade.

DISCUSSION
In this study, we have established an approach, MAB-seq, for quantitative measurement of 5fC/5caC excision repair-dependent active DNA demethylation activity at base-resolution, providing a tool that has broad application in the study of active DNA demethylation. A major advantage of MAB-seq over other subtraction-based 5fC or 5caC mapping methods is the ability to directly determine the location and abundance of 5fC and 5caC in a single experiment. MAB-seq can therefore achieve similar sensitivity for detecting 5fC/5caC with lower sequencing depth and simplifies the computational procedure to confidently identify CpGs marked 5fC/5caC. When combined with Tdg depletion, simultaneous mapping 5fC and 5caC enables quantitative measurement of the generation and excision of 5fC/5caC, providing a direct readout of the TET/TDG-mediated active DNA demethylation activity.
Application of whole genome MAB-seq analysis to mouse ESCs allowed us to investigate the genomic architecture and dynamics of active DNA demethylation activity at single-base resolution across the genome of this stem cell population. While affinity-based 5fC and 5caC mapping methods suggest that 5hmC and 5fC/5caC enriched regions are largely overlapped 18,19 , comparative analysis of base-resolution 5fC/5caC map and 5hmC map reveals that a large fraction of 5fC/5caC and 5hmC are associated with distinct CpGs within 5hmC/5fC/5caC-enriched regions. These results support a model in which TET proteins tend to stall at 5hmC step at most CpGs, but exhibit higher processivity to further oxidize 5hmC to 5fC/5caC at CpGs with higher chromatin accessibility. Moreover, integrating MAB-seq with 5fC-mapping (e.g. redBS-seq) or 5caC-mapping (e.g. caMAB-seq) method not only provides a base-resolution approach to map 5fC and 5caC separately, but also reveals that 5fC and 5caC are frequently non-overlapped at individual CpGs. These findings suggest the possibility that TET exhibits distinct processivity depending on local chromatin accessibility or other yet to be identified regulatory processes. It is also possible that TDG may exhibit distinct activity or substrate preference at different CpGs. In addition, simultaneous mapping of 5fC/5caC by MAB-seq reveals that TET/TDG-dependent active DNA demethylation activity preferentially targets palindromic CpG dyads asymmetrically, supporting the recently proposed asymmetric base-flipping model 24,36 . Lastly, we demonstrate that integrative analysis of MAB-seq and BS-seq datasets not only allows quantitative measurement of cytosine derivatives at all major steps of TET/TDG-mediated active demethylation pathway [iterative oxidation (5mC/5hmC), excision repair (5fC/5caC) and regeneration (unmodified C)], but also reveals that 5fC/5caC excision by TDG may act as a regulatory checkpoint of the active DNA demethylation cascade.
In addition to a robust base-resolution mapping method for 5fC/5caC, we provide here a whole-genome base-resolution map of TET/TDG-dependent active DNA demethylation activity in the mammalian genome. Given the simplicity and cost-effectiveness of MAB-seq, we anticipate that genome-scale MAB-seq analysis can be applied to analyze diverse celltypes in future studies to provide new insights into the mechanism and function of the DNMT-TET-TDG/BER cytosine-modifying cascade. Thus, the 5fC/5caC mapping technology described in this study and other base-resolution mapping methods set the stage for systematic investigation of the functional significance of active DNA demethylation in mammalian development and human diseases.

Data processing of genome-scale MAB-seq/BS-seq
Raw sequencing reads were trimmed for low-quality bases and adaptor sequences using Trimmomatic 37 , and the data quality was examined with FastQC (http:// www.bioinformatics.babraham.ac.uk/projects/fastqc/). The trimmed reads were mapped against the mouse genome (mm9 build) with Bismark 38 . PCR duplicates were removed using the Picard program (http://picard.sourceforge.net/). Methylation call for each CpG was extracted using bismark methylation extractor script from Bismark. Cell line-specific SNPs overlapping with CpGs in the mouse genome (mm9) were filtered out with BisSNP 39 . All programs were performed with default setting. For MAB-seq analysis, raw signals were calculated as % of T/(C+T) at each called CpG. For BS-seq analysis, raw signals were calculated as % of C/(C+T) at each called CpG. Statistics of all genome-scale sequencing libraries is summarized in Supplementary Table 1.

Statistical calling of 5fC/5caC and assessing false discovery rate (FDR) of MAB-seq
For each cytosine within CpG dinucleotides, we counted the number of "T" bases from MAB-seq reads as 5fC/5caC (denoted N T ) and the number of "C" bases as other forms of cytosines (C/5mC/5hmC; denoted N C ). Next, we used the binomial distribution (N as the sequencing coverage (N T + N C ) and p as the error rate (2.04%) of M.SssI methylase) to assess the probability of observing N T or greater by chance. To estimate empirical FDR of calling 5fC/5caC-modified CpGs, we repeated the steps above on MAB-seq signals of merged negative control sample (Dnmt1/3a/3b −/− and Tet1/2 −/− ) in which 5fC/5caC is largely absent. The FDR for a given P-value cutoff of the binomial distribution is the number of 5fC/5caC-called CpGs divided by the number detected in the sample of Tdgdeficient VC-treated mouse ESCs. For calling 5fC/5caC-modified CpG dyads (Fig 2c), reads covering Watson strand and Crick strand of the same CpG dyads were combined. We restricted our analysis to CpG dyads covered by at least five reads. For Fig 3 and 4, strandspecific analysis was performed and we restricted our analysis to CpGs covered by at least ten reads. For Fig 3a, called 5fC/5caC-modified CpGs were visualized in Integrative Genomics Viewer (IGV).

Validation of MAB-seq by synthetic 38-bp oligonucleotides
To monitor the behavior of C and modified C in MAB-seq, 38-bp single-stranded DNA oligos were synthesized (forward strand: 5'-AGCCXGXGCXGXGCXGGTXGAGXGGCXGCTCCXGCAGC-3', reverse (complementary) strand: 5'-GCTGXGGGAGXGGCXGCTXGACXGGXGXGGXGXGGGCT-3', in which X is either unmodified C, 5hmC, 5fC or 5caC). To test whether M.sssI treatment alters the behavior of 5hmC, 5fC and 5caC during bisulfite sequencing, forward strands containing 5hmC, 5fC and 5caC were annealed to reverse strands containing the same modified Cs and ligated to methylated adaptors. The resulting oligos were treated by M.sssI using the same protocol used for locus-specific MAB-seq (described below), followed by bisulfite conversion and deep sequencing to determine their behavior. To test whether M.sssI functions at hemi-5hmC, 5fC and 5caC CpGs, top strands containing 5hmC, 5fC and 5caC were annealed to bottom strand containing unmodified C. The same experimental procedures were undertaken to assess the efficiency of M.sssI in these contexts.

Locus-specific MAB-seq of genomic DNA
1 µg genomic DNA was treated by M.sssI in a 50 µL reaction for four rounds. During each round of treatment, DNA was first treated by M.sssI for 2 hours (1.5 µL M.sssI and 1 µL SAM), and additional M.sssI and SAM were supplemented to treat DNA for another 4 hours (0.5 µL M.sssI, 1 µL SAM), increasing the total concentration of the enzyme to 0.8 unit/µL. DNA was purified by PCI after each round of treatment.
After M.sssI treatment, bisulfite conversion was performed as described above. Selected loci were amplified by PCR using KAPA HiFi Hotstart Uracil+ DNA polymerase (Kapa Biosystems, KK2801) followed by sonication by Bioruptor (Diagenode), library preparation using NEBNext DNA Library Prep Master Mix Set (New England Biolabs) and deep sequencing by illumina HiSeq 2500 sequencer. Alternatively, PCR amplified DNA was cloned into TOPO vectors (Zero Blunt TOPO Cloning Kit, Invitrogen) for standard Sanger sequencing. Primer sequences for all locus-specific MAB-seq experiments were summarized in Supplementary Table 2.

5fC/5caC calling for locus-specific MAB-seq
In a MAB-seq experiment, 5fC and 5caC are read as T while unmodified C, 5mC and 5hmC are read as C. For a CpG site, if we name the number of T reads as N T and the number of C reads as N C , 5fC+5caC level (MAB-seq raw signal before background subtraction) can be calculated as N T /(N T +N C ). Because conversion of unmodified C to 5mC by M.sssI cannot reach 100%, some T reads may come from incomplete methylation rather than real 5fC/ 5caC, leading to an overestimation of 5fC/5caC level. To estimate the level of these background signals resulted from incomplete conversion, Tet1/2 −/− mouse ESCs was examined by MAB-seq, and any 5fC/5caC signal detected in this negative control is treated as background signal (false positive).
To call 5fC/5caC-positive CpG sites in a locus-specific MAB-seq experiment, MAB-seq signals detected in shCtrl, shCtrl+VC, shTdg and shTdg+VC samples were compared to the background signals detected in Tet1/2 −/− sample. For a CpG site in a tested sample, Fisher's exact test was performed using NC and NT of the tested sample and Tet1/2 −/− sample to determine whether MAB-seq signal at this CpG site is significantly different from the corresponding background signal. False discovery rate (FDR) control based on Benjamini-Hochberg procedure was then performed to correct the p values for multiple comparisons, and p<0.05 was used as a cut-off value to generate a list of CpG sites of which MAB-seq signals are significantly different from the background signals detected in Tet1/2 −/− sample.
After that, we further applied the numeric filter that real 5fC/5caC signals should be numerically higher than background signals, generating the list of CpG sites that are 5fC/ 5caC positive.

Calculating 5fC/5caC level in locus-specific MAB-seq
In Tet1/2 −/− sample, certain background (false-positive) signals were indeed detected by MAB-seq. To counteract these noises introduced by the method, these background signals were subtracted from raw MAB-seq signals when calculating real 5fC/5caC levels. For example, if raw MAB-seq signal detected at a CpG site in our sample of interest is a% while the corresponding background signal detected in Tet1/2 −/− sample is b%, and if a≥b, then (a-b)% will be the 5fC/5caC level. If a<b is observed (in some rare cases), then 5fC/5caC level at that CpG site in our sample of interest will be set as 0%. To determine whether the 5fC/5caC level at a group of CpG sites in one sample is different from that in another sample, Wilcoxon signed-rank test for matched pairs was performed.

Analysis of strand asymmetry of 5fC/5caC in locus-specific MAB-seq
For Fig 4d-e, 11 palindromic CpG dyads from Tbx5 locus and 4 palindromic CpG dyads from Vps26a locus were examined by locus-specific MAB-seq to determine whether strand asymmetry of 5fC/5caC exists. 5fC/5caC calling was achieved through performing Fisher's exact test, and the 15 CpG dyads were first categorized into three groups: no significant 5fC/ 5caC signal on either the top and bottom strands; significant 5fC/5caC signal on one strand but not the other; significant 5fC/5caC signals on both strands. The existence of the second group is a direct support for strand asymmetry of 5fC/5caC. As for the third group of CpG dyads, Fisher's exact test was performed to determine whether 5fC/5caC levels differ significantly between the top and bottom strands, and this group was further separated into two groups. For each CpG dyad, if 5fC/5caC levels differ significantly between strands, strand asymmetry of 5fC/5caC exists at this dyad. If not, other methods such as hairpin MAB-seq may help to conclusively determine the state of 5fC/5caC.

Data analysis of locus-specific BS-seq
In BS-seq experiment, unmodified C, 5fC and 5caC are read as T while 5mC and 5hmC are read as C. For a CpG site, if we name the number of T reads as N T and the number of C reads as N C , 5mC+5hmC level at this CpG site (before background subtraction) will be calculated as N T /(N T + N C ). To determine whether the 5mC+5hmC level at a CpG site is altered by Tdg knockdown or vitamin C treatment, Fisher's exact test was performed. To determine whether the 5mC+5hmC level at a group of CpG sites in one sample is different from that in another sample, Wilcoxon signed-rank test for matched pairs was performed.

Calculating the level of unmodified cytosines by combining locus-specific MAB-seq with BS-seq
For a CpG site in our sample of interest, the true level of 5fC+5caC was calculated by subtracting the background signal from the raw MAB-seq signal as described above. The level of 5hmC+5mC was measured by BS-seq. The level of unmodified C equals to 100% -(abundance of 5fC+5caC) -(abundance of 5hmC+5mC). To determine whether the unmodified C level at a group of CpG sites in one sample is different from that in another sample, Wilcoxon signed-rank test for matched pairs was performed.

Locus-specific redBS-seq of genomic DNA
redBS-seq was performed as previously described 21 . In short, 5 µL freshly made sodium borohydride aquerous solution (1M) was added to 250 ng DNA diluted in 15 µL water. The reaction was placed in darkness for an hour, quenched by 10 µL sodium acetate (0.75M, pH=5) and purified by PCI extraction. Bisulfite conversion was then performed as described above, followed by PCR amplification of specific loci and deep sequencing analysis of PCR amplicons.

Analysis of locus-specific redBS-seq data
In a redBS-seq experiment, 5mC, 5hmC and 5fC are read as C while unmodified C and 5caC are read as T. For a CpG site, the difference between the levels of C signal measured by redBS-seq and BS-seq was calculated to represent the level of 5fC. In rare cases when negative values were obtained after this subtraction, 5fC levels were designated as 0%. To call CpG sites with significant levels of 5fC, Fisher's exact test was performed to determine which CpG sites have significantly higher C signals in redBS-seq compared with BS-seq (FDR corrected p value less than 0.05). To calculate the level of 5caC, the level of 5fC measured by redBS-seq was subtracted from the level of 5fC+5caC measured by MAB-seq, and negative values resulted from the subtraction were designated as 0%. To reduce the noise resulted from the two subtractions (one for 5fC calculation and one for 5caC calculation), uncalled sites in redBS-seq or MAB-seq were regarded as having 0% 5fC or 5fC+5caC, respectively.

Locus-specific caMAB-seq and 5caC calling
To perform locus-specific caMAB-seq, genomic DNA was firstly treated by M.sssI as described above and purified through PCI extraction. M.sssI-treated DNA was then reduced by NaBH 4 as described above in the redBS-seq protocol, followed by bisulfite conversion using Qiagen Epitect Bisulfite Kit. Selected loci were then amplified by PCR, and PCR amplicons were sequenced using Illumina deep sequencing. To examine whether unmodified C and 5fC are read as C during caMAB-seq, unmodified lambda DNA and synthesized 5fC oligo were also treated using the same protocol and analyzed by Illumina deep sequencing. In caMAB-seq, 5caC is read as T while C, 5mC, 5hmC and 5fC are read as C. To call 5caC-positive sites in a locus-specific caMAB-seq experiment, Fisher's exact test was performed to determine whether a CpG site in a sample has significantly higher percentage of T compared with background T signals detected in Tet1/2 −/− sample. To calculate the absolute level of 5caC in a locus-specific experiment, the background signal detected in Tet1/2 −/− sample was subtracted from the raw signal detected in a tested sample.

Genome-wide 5mC/5hmC/5fC/5caC DIP-Seq
The antisera for 5fC and 5caC were previously described 40 . For each DIP experiment, 10 µg of sonicated, adaptor ligated genomic DNA from control or Tdg knockdown mouse ESCs (V6.5) was used as input, and 5 µl of 5mC antibody (Eurogentec, BI-MECY-0500), 5 µL of  (a) Schematic diagram of MAB-seq. DNMT methylates unmodified C to generate 5mC, which can be successively oxidized by TET to generate 5hmC/5fC/5caC. Highly oxidized cytosine derivatives, 5fC and 5caC, are repaired by TDG/BER to regenerate unmodified C. (b) M.SssI exhibits robust methylase activity towards unmodified cytosines within CpGs, but has significantly lower activity for cytosines in non-CpG contexts (CHG or CHH, H=A, T, C). M.SssI methylase activity was measured by MAB-seq analysis (Illumina deep sequencing) of unmethylated lambda DNA. Standard bisulfite sequencing (BS-seq) confirmed nearly complete conversion of unmethylated C to T at CpG and non-CpG sites. (c) Locus-specific MAB-seq analysis of 5fC/5caC at Tbx5 by Illumina sequencing in control (shCtrl) and Tdg knockdown (shTdg) mouse ESCs. For comparison, also shown are 5fC/ 5caC antibody DNA immunoprecipitation (DIP) based maps of 5fC and 5caC in control and Tdg-depleted mouse ESCs. DIP-seq tracks are represented in normalized read density (reads per 10 million reads, rp10m) and the vertical axis range of all DIP-seq tracks is from 1 to 25. Black horizontal bars denote 5fC/5caC-enriched regions identified by 5fC/5caC DIP-seq. The level of 5fC/5caC (only Watson strand shown) is displayed as the percentage of total C modified as 5fC/5caC, and background signals detected in Tet1/2 −/− mouse ESCs was subtracted. (d) Statistical calling of 5fC/5caC-modified CpGs in locus-specific MAB-seq analysis. Shown are 5fC/5caC levels (background corrected using raw MAB-seq signals in Tet1/2 −/− ) of 11 CpG sites from the Tbx5 locus in shCtrl, shCtrl+VC, shTdg and shTdg+VC mouse ESCs. shTdg was compared with shCtrl while shTdg+VC was compared with shCtrl+VC. An asterisk indicates that a CpG is statistically enriched for 5fC/5caC (multiple comparison corrected P<0.05, Fisher's exact test).    (a) The average 5fC/5caC level on the called cytosines (n= 454,400, identified in WG-MAB-seq; FDR=5%; sequencing depth on both strands≥10) is significantly higher than that of cytosines on the opposite strand. Called, cytosines enriched for 5fC/5caC in shTdg+VC; opposite, cytosines on the opposite strand. (b) Distribution of differences of WG-MAB-seq signals (in shTdg+Vc) between called and opposite cytosines. Symmetric, both strands are called for enriching 5fC/5caC in shTdg+VC; Asymmetric, only one strand is called. (c) Two-dimensional density-plot of WG-MAB-seq signals at each called/opposite CpG dyad reveals the strand asymmetry of TET/TDG-dependent active DNA demethylation. (d) Locus-specific MAB-seq analysis of two representative loci illustrating the asymmetric distribution of 5fC/5caC. 5fC/5caC levels for the Watson strand are in blue, whereas those for the Crick strand are in red. (e) Bar-graph of relative percentage of CpG dyads (n=15, analyzed by locus-specific MABseq) that are associated with different degrees of strand asymmetry of active DNA demethylation activity in shCtrl, shCtrl+VC, shTdg, shTdg+VC mouse ESCs. (a) Schematic diagram of combining BS-seq, redBS-seq, MAB-seq and caMAB-seq to map 5fC and 5caC individually at single-base resolution. (b) Locus-specific analysis of 5fC and 5caC at the Tbx5 locus indicates that 5fC and 5caC are largely non-overlapping and both modifications exhibit strong strand asymmetry. For comparison, affinity-based 5fC and 5caC maps are shown on the top. In enlarged views, base-resolution maps of 5fC+5caC (measured by MAB-seq), 5fC (measured by redBS-seq) Figure 6. Base-resolution mapping of unmodified C and 5mC/5hmC reveals that TDG-mediated 5fC/5caC excision affects active demethylation dynamics (a) Tdg knockdown and/or VC treatment results in a significant decrease in BS-seq signals (5mC+5hmC). Shown are boxplots summarizing the levels of 5hmC+5mC at 73 CpG sites from Ace, Slc2a12, Tbx5 and Vps26a loci. An asterisk indicates a statistically significant difference between the two indicated groups, as determined by Wilcoxon signed-rank test for matched pairs (multiple comparison corrected P-value<0.05). (b) Shown are relative changes in 5mC+5hmC levels (log 2 ratio) in response to Tdg knockdown at individual CpG sites from Ace and Tbx5 loci. shTdg was compared with shCtrl while shTdg+VC was compared with shCtrl+VC. An asterisk indicates a statistically significant change determined by Fisher's exact test (multiple comparison corrected P-value<0.05). (c) Tdg knockdown or VC treatment leads to a significant increase in the level of unmodified C (measured by subtracting C signals in BS-seq from those of MAB-seq). Shown are boxplots summarizing the levels of unmodified C at 73 CpG sites from Ace, Slc2a12, Tbx5 and Vps26a loci. (d) Shown are relative changes in unmodified C levels (log 2 ratio) in response to Tdg knockdown at individual CpG sites from Ace and Tbx5 loci. shTdg was compared with shCtrl while shTdg+VC was compared with shCtrl+VC. (e) Genome-scale H3K27me3-MAB-seq and H3K27me3-BS-seq analyses in control or Tdgdepleted mouse ESCs show that 5m/5hmC levels are significantly decreased within genomic regions enriched for 5fC/5caC. Boxplots depicting MAB-seq and BS-seq signals (binned to 500bp intervals) are shown for 5fC/5caC-enriched regions. Equal number of randomly sampled regions that are covered by both BS-seq and MAB-seq are analyzed. (f) Relative changes (log 2 ratio) in 5mC+5hmC levels (top panels) and unmodified C levels (bottom panels) in response to VC treatment. shCtrl+VC was compared with shCtrl while shTdg+VC was compared with shTdg. For relative changes in 5mC+5hmC levels, an asterisk indicates a statistically significant change determined by Fisher's exact test (multiple comparison corrected P-value<0.05).