Sequencing DNA methylation and hydroxymethylation at co-occurring chromatin features

Domain: neuroscience
Published: 2026-02-10
Editorial outcome: accept

Abstract
Epigenetic modifications govern chromatin dynamics and cell state. However, current methods cannot simultaneously resolve the presence of multiple DNA modifications at co-occurring chromatin-associated features. It is thus not clear how these features are physically coupled and how their combinations regulate genome function. To address this key question, we report 6-base-CUT&Tag, a method for simultaneous 6-base DNA sequencing at target chromatin features. Using 6-base-CUT&Tag to profile 5-methylcytosine (5mC) and 5-hydroxymethylcytosine (5hmC) at co-occurring histone modifications in mouse embryonic stem cells (mESCs), we identify feature-dependent 5mC/5hmC signatures previously unresolvable with untargeted or bisulfite-based workflows. We show that DNA methylation and hydroxymethylation are specifically coupled with the H3K4me1 mark in mESC enhancers and that H3K4me1-derived signatures robustly distinguish different enhancer functional states.

Introduction
Chemical modifications of DNA bases and histone tails regulate physical and functional properties of chromatin, including DNA accessibility, genomic architecture, and binding of transcription regulators1–4. Methylation of cytosine bases on DNA is critical to establish the epigenetic landscape and governs gene activity and inheritance5. Cytosine bases can be methylated by DNA methyltransferases (DNMTs), that transfer a methyl group from the donor S-adenosyl methionine to form 5-methylcytosine (5mC)6, which can be further oxidised to form 5-hydroxymethylcytosine (5hmC) through the action of the Ten-Eleven Translocation family of proteins (TETs1-3)7. TET enzymes can further oxidise 5hmC to other intermediates like 5-formylcytosine (5fC) and 5-carboxylcytosine (5caC)8, each of which can be excised by thymine DNA glycosylase, followed by restoration of C, completing a demethylation cycle (Fig. 1a). 5mC is the most abundant cytosine modification (2–4% of all Cs in mouse embryonic stem cells [mESCs]), and 5hmC is the next most abundant cytosine modification (0.1–0.2% of all Cs in mESCs), with 5fC and 5caC being present at substantially lower levels (~0.002% and ~0.0003%, respectively)8. While 5hmC can be an intermediate inactive DNA demethylation, it has also been shown to be a stable epigenetic mark that is associated with transcriptional activity and cell fate9–14. The functional roles of 5mC and 5hmC depend on their genomic context. Broadly, 5mC enrichment at promoters and enhancers is associated with transcriptional repression15,16, whereas 5hmC is enriched at active regulatory regions—particularly enhancers—and also within the gene bodies of transcriptionally active genes17,18.

The development of sequencing methods to detect 5mC and 5hmC has been an important focus in epigenomics. While early methods such as bisulfite sequencing19 cannot distinguish 5hmC from 5mC, later approaches—including oxidative bisulfite sequencing20, TET-assisted bisulfite sequencing21, APOBEC-coupled epigenetic sequencing22, pyridine borane sequencing methods23,24 and third-generation sequencing platforms (e.g., PacBio, Oxford Nanopore)25– have enabled base-resolution mapping of both 5mC and 5hmC. More recent methods have now leveraged state-of-the-art enzymology and chemistry to enable bisulfite-free and simultaneous detection of 5mC and 5hmC on the same DNA fragment26–28.

Despite such advancements, a key challenge in the field has been to determine how 5mC and 5hmC are coupled in space and time to other chromatin features such as histone marks, transcription factors and chromatin remodelers. Establishing this relationship is critical to appreciate how epigenetic features are coupled, when and where they co-occur, and how they determine causal relationships and regulatory roles. Current strategies that map DNA modifications and chromatin features largely rely on comparing profiling data from independent experiments for different features. This approach does not necessarily determine the co-occurrence of these features at the molecular level (Fig. 1b). Earlier attempts to address this limitation have coupled chromatin immunoprecipitation sequencing (ChIP-seq) and cleavage under targets and tagmentation (CUT&Tag) with bisulfite treatment to detect cytosine methylation on DNA fragments directly associated with chromatin proteins29–32. However, these methods have major limitations that include substantial loss of DNA caused by bisulfite-mediated degradation33, which precludes analysis of smaller sample quantities. Crucially, the inability to discriminate between C, 5mC, and 5hmC bases at read level in the same workflow34 is a major gap in our understanding of how combinations of epigenetic features influence genome regulation.

Here, we introduce 6-base-CUT&Tag sequencing (abbreviated 6B-C&T), a method to simultaneously map G, A, T, C, 5mC, and 5hmC bases at co-occurring chromatin-associated features. 6B-C&T begins with targeted antibody enrichment of specific chromatin features via CUT&Tag that are purified from a heterogeneous pool, followed by an updated 6-base library preparation workflow, sequencing, and analysis (Fig. 1c). This approach substantially reduces sequencing depth requirements and cost relative to whole-genome sequencing analysis (see “Methods”, Supplementary Data 1). We demonstrate that 6B-C&T reveals 5mC/5hmC signatures which are feature-dependent and have not been recoverable in previous untargeted or bisulfite-based workflows.

Results
6-base-CUT&Tag (6B-C&T) development
As proof-of-concept for 6-base-CUT&Tag (6B-C&T), we used E14TG2A mouse embryonic stem cells (mESCs) and initially chose H3K27ac, a highly prevalent histone mark in active chromatin, as the feature to profile (Fig. 1c). Using an antibody selective for H3K27ac, we first performed antibody-directed pA-Tn5 tagmentation with a uracil-containing hairpin mosaic-end adaptor, termed ME2U. Tagmented DNA is then subjected to circularisation to generate dumbbell-like molecules that are resistant to exonuclease treatment. Uracil-specific excision reagent (USER) treatment digests the uracil sites of the adaptor, allowing the strands to separate and the residual hairpin sequences on each strand to refold independently. This produces DNA suitable for our 6-base-seq enzymatic conversion workflow26 (Fig. 1c, Supplementary Fig. 1). After Illumina paired-end sequencing, the four canonical DNA bases and two cytosine epigenetic states can be resolved on individual reads to provide the G, A, T, C, 5mC, and 5hmC sequence context at H3K27ac sites.

An important element in the development of 6B-C&T was to avoid interference from partially tagmented or randomly cleaved DNA, which would be a major source of non-specific background. To address this, we took advantage of the ‘scar’ sequences left on both strands of the DNA fragment after USER digestion of the inserted ME2U hairpin. These residual sequences are copied and deaminated prior to PCR amplification (Fig. 2a) and thus serve as an internal benchmark to identify valid read pairs derived from double-sided hairpin constructs. Furthermore, circularisation of such constructs makes them resistant to exonuclease digestion, as opposed to non-specifically cleaved or partially tagmented DNA fragments. Motivated by related strategies to remove similar DNA fragments35,36, we determined that a mixture of T7 Exonuclease (dsDNA-specific), RecJf, and Exonuclease I (ssDNA-specific) achieved the highest recovery of scarred read pairs relative to a protocol with no enzyme treatment or with other enzyme combinations (Fig. 2a).

In our initial H3K27ac 6B-C&T in mESCs, scar-filtered read pairs mapped to regions that strongly overlapped with H3K27ac sites determined independently using standard CUT&Tag (Fig. 2b, c; Supplementary Figs. 2, 3). To verify the accuracy of DNA modification measurement, we included three different synthetic spike-ins as fiducials (Methods) that were also subjected to 6-base-seq base conversions alongside tagmented mESC DNA. For each DNA control (5mC-modified Lambda DNA, 5hmC-modified oligonucleotide and unmodified pUC19 DNA), we achieved high sensitivity (98-99.5%) in identifying the true cytosine modification state (Fig. 2d). These results show that 6B-C&T can faithfully recapitulate the genomic distribution of a target feature and simultaneously detect 5mC and 5hmC bases with high accuracy.

Having validated 6B-C&T on H3K27ac, we set out to profile 5mC and 5hmC in mESCs at each of several histone modifications marking distinct types of regulatory elements37. For this, we included H3K27ac (active enhancers and promoters), H3K4me3 (promoters), H3K4me1 (enhancers), and H3K27me3 (Polycomb-repressed, bivalent chromatin), along with a non-specific IgG control. As a reference, we also performed whole-genome 6-base-seq, which provides an untargeted picture representative of the total pool of DNA fragments. Using this approach, we confirmed that the global relative 5mC/5hmC abundance on genomic DNA was consistent with reported values in mESCs (Supplementary Fig. 4a). Genomic enrichments obtained by 6B-C&T were highly reproducible across experimental replicates for each histone mark (Pearson’s r≥0.97) and showed good agreement (Pearson’s r≥0.86) when compared to standard CUT&Tag (Supplementary Fig. 5). CpG modification levels (i.e. the fraction of CpGs detected as 5mC or 5hmC) at each histone mark also showed excellent reproducibility (Pearson’s r≥0.90, Supplementary Fig. 6). Finally, we assessed how well 6B-C&T correlates with CUT&Tag-BS, a method that profiles total cytosine methylation (modC) on CUT&Tag DNA fragments via bisulfite conversion32. Both genomic feature enrichment and CpG modification levels for 6B-C&T agreed well with publicly available CUT&Tag-BS data for H3K4me1 and H3K27me3 (Pearson’s r≥0.83, Supplementary Fig. 7).

Measuring 5mC and 5hmC at distinct histone modifications
Methylation is known to be depleted in histone mark-enriched regulatory regions relative to the other regions in the mESC genome38. Consistent with this, we observed that average 5mC and 5hmC levels in 6B-C&T were lower than in the whole-genome (untargeted) 6-base-seq experiment (Supplementary Fig. 4a). This difference is due to the 6B-C&T method specifically capturing DNA fragments at histone mark-enriched sites, which are generally less (hydroxy)methylated compared to the whole genome. Furthermore, when looking at individual histone marks, we found that those associated with active chromatin (H3K4me3, H3K27ac, H3K4me1) showed lower overall methylation and a higher 5hmC/5mC ratio than the bivalent/repressed chromatin (H3K27me3) mark, a pattern that was absent in whole-genome data (Supplementary Fig. 4b, Supplementary Fig. 8). Globally, individual CpG sites captured by 6B-C&T mapped to regions with high histone mark enrichment and showed differences in 5mC and 5hmC levels relative to the untargeted method (Supplementary Fig. 9).

We next investigated how 5mC and 5hmC levels in 6B-C&T differed from untargeted 6-base-seq data for the same CpG sites (Supplementary Fig. 10a). First, we identified CpG sites with statistically significant differences between both methods, i.e., differentially methylated CpGs (DMCs) and differentially hydroxymethylated CpGs (DHMCs). While most DHMCs between both methods were located in introns followed by intergenic and promoter regions, DMCs were primarily found in promoters followed by introns/exons (Supplementary Fig. 10b–e). Overall, DMCs predominantly showed lower 5mC levels in 6B-C&T relative to whole-genome (untargeted) data (Supplementary Fig. 10a, Supplementary Fig. 11). This suggests that DNA associated with histone modifications at these sites is mostly hypomethylated relative to the total pool of DNA fragments. This was particularly striking in DMCs co-occurring with H3K4me3 at imprinted gene loci39 such as Igfr2, Kcnq1ot1 and Peg3, suggesting that this active histone mark largely associates with the unmethylated allele (Supplementary Fig. 10a). For DHMCs, active chromatin marks (H3K4me3, H3K27ac and H3K4me1) showed similar numbers of sites with higher or lower 5hmC levels in 6B-C&T relative to the untargeted experiment; in contrast, for the repressive mark (H3K27me3) the majority of DHMCs showed higher 5hmC levels in 6B-C&T than in untargeted 6-base-seq (Supplementary Fig. 10a, Supplementary Fig. 11). Overall, these observations reveal important epigenetic differences between the DNA associated with target histone modifications and the bulk average at the same genomic sites. These quantitative relationships would be difficult to discern from comparative analysis of independent measurements of 5mC, 5hmC, and histone marks.

Histone mark-specific 5mC and 5hmC enrichment at enhancers
To demonstrate the applicability of 6B-C&T to determine the true base modification status at specific chromatin elements, we analysed a set of mESC enhancers that have been independently annotated by the presence or absence of specific histone marks by ChIP-seq40, including active (H3K27ac, H3K4me1), primed (H3K4me1), and poised (H3K4me1, H3K27me3) enhancers. We observed good agreement at these loci between 6B-C&T and the ChIP-seq data originally used for their annotation (Supplementary Fig. 12). Excellent agreement was also obtained for total DNA methylation levels (%modC) at all enhancer types between 6B-C&T and available C&T-BS data for the general enhancer mark, H3K4me1 (Supplementary Fig. 7c). We next found that different enhancer loci displayed distinct patterns of 5mC and 5hmC enrichment (Fig. 3a, Supplementary Fig. 13). Across all enhancer types, primed enhancers, defined by H3K4me1 alone, generally showed the highest levels of both 5mC (average ~13%) and 5hmC (average ~4%; Fig. 3b). It is notable that this increased enrichment at the H3K4me1 mark was specific to primed enhancers, since other H3K4me1-associated CpG sites showed much lower 5mC and 5hmC levels (Supplementary Fig. 14; Supplementary Data 1—Table 1). In contrast, active and poised enhancers were relatively depleted of both 5mC (~4% for active and ~2% for poised, on average) and 5hmC (~2% for active and poised, on average). For all enhancer types, H3K4me1 consistently showed the highest 5mC and 5hmC levels and the highest 5mC/5hmC ratio among all histone marks (Fig. 3b).

We then assessed whether the enrichment of 5mC and 5hmC at primed enhancers was specific to H3K4me1 by measuring these DNA modifications at other enhancer-associated histone marks potentially present at these sites, i.e., H3K27ac and H3K27me3. We first confirmed that annotated primed enhancers were enriched with the H3K4me1 mark but also showed appreciable levels of H3K27ac and H3K27me3 (Fig. 3c, Supplementary Fig. 15). Across all three histone marks, H3K4me1 displayed the highest fraction of both 5mC and 5hmC, and the highest 5mC/5hmC ratio (Fig. 3d, Supplementary Fig. 16). This pattern, which cannot be detected with previously reported methods like ChIP-BS41, suggests that H3K4me1-marked nucleosomes are sites for preferential 5mC retention at enhancers. Active (H3K27ac) and bivalent (H3K27me3) histone marks, which do not represent the primed enhancer state, showed a greater tendency towards demethylation (lower overall methylation levels, higher 5hmC/5mC ratio). These results suggest that 5mC and 5hmC are coupled to different extents with each chromatin mark and, more generally, that different enhancer states are defined by distinct combinations of epigenetic modifications.

5mC and 5hmC at H3K4me1 resolved distinct enhancer functional states
H3K4me1 is common to all enhancer functional states, so we next asked whether 6B-C&T methylation at this single feature better resolves different enhancer types compared to the more heterogeneous whole genome (untargeted) 6-base-seq data. For this, we used H3K4me1 6B-C&T data to train a machine learning model for classification of enhancer types (active, primed and poised) and compared its performance to a model with the same architecture trained on whole-genome 6B-seq data (Fig. 4a, Supplementary Fig. 17). For all enhancer types, the 6B-C&T-derived model showed superior performance when compared to the whole genome-derived model (Fig. 4b–d; Supplementary Fig. 18). This demonstrates that the H3K4me1-coupled 5mC and 5hmC signal more robustly resolves bona fide enhancer states than its unenriched equivalent. These results show that 6B-C&T can be used to evaluate co-occurring epigenetic features and retrieve improved information on individual molecular states (Fig. 1b).

Methods
Cell culture
E14TG2A mouse embryonic stem cells were cultured in serum/LIF conditions on gelatin-coated plates. Serum/LIF medium was prepared with high-glucose DMEM (Sigma, D6546-500mL) supplemented with 10% FBS (Gibco, 16141079), 1X GlutaMax-I (Gibco, 35050-038), 1X NEAA (Gibco, 11140-035), 0.1mM β-mercaptoethanol (Sigma, M3148) and 1000U/mL LIF (PeproTech, murine LIF, 250-02-25UG). For maintenance, cells were grown until 70–80% confluence and subcultured by dissociation with StemPro Accutase Cell Dissociation Reagent (Invitrogen A1110501). Cultures were routinely monitored and tested negative for Mycoplasma contamination. E14TG2A cells were a gift from Wolf Reik (Babraham Institute, Cambridge).

Preparation of tagmentation adaptors
Tn5ME-A, Tn5ME-B and Tn5MErev oligonucleotides were reconstituted to 200µM in annealing buffer (10mM Tris-HCl pH 8.0, 50mM NaCl, 1mM EDTA). Tn5ME-A and Tn5ME-B were mixed separately with equal amounts of Tn5MErev and annealed on a thermocycler (95°C for 2min, then cooled in 5°C steps with 5min incubation at each step until 25°C). The resulting mosaic-end adaptors are termed MEDS-A and MEDS-B, respectively. The ME2U hairpin adaptor was reconstituted to 25µM in Nuclease-Free Duplex Buffer (30mM HEPES, pH 7.5; 100mM potassium acetate, IDT 11-01-03-01) and annealed by heat denaturation at 95°C for 2min, followed by slow cooling at 0.1°C/second to 10°C. All oligonucleotide sequences (including chemical modifications) are listed in Supplementary Data 1 (Table 2).

pA-Tn5 transposome assembly
pA-Tn5 transposase was prepared in-house43 or obtained commercially (Active Motif, 53162) and loaded with application-specific DNA adaptors. For the standard CUT&Tag workflow, 50µL of pA-Tn5 (4µM) were mixed with a threefold excess of each pre-annealed mosaic-end adaptor, MEDS-A and MEDS-B, and rotated at room temperature for 1h. For 6-base-CUT&Tag, 50µL of pA-Tn5 (4µM) were mixed with a threefold excess of folded ME2U hairpin adaptor and similarly rotated at room temperature for 1h. Transposome samples were stored at −20°C.

Standard CUT&Tag (cleavage under targets and tagmentation)
Standard CUT&Tag was performed as in Kaya-Okur et al.44 with modifications. Briefly, Dynabeads MyOne Streptavidin T1 beads (Invitrogen 65601) were conjugated with biotinylated ConA45 (Sigma C2272) and activated in binding buffer (20mM HEPES pH 7.5, 10mM KCl, 1mM CaCl2, 1mM MnCl2) prior to cell binding. 500,000 cells per assay were lightly fixed for 2min with 0.1% formaldehyde in PBS and fixation was stopped by adding glycine to a final concentration of 120mM. Fixed cells were resuspended in wash buffer (20mM HEPES pH 7.5, 150mM KCl, 0.5mM spermidine, and Complete EDTA-Free Protease Inhibitor, Roche 11873580001) and immediately bound to ConA-conjugated MyOne T1 beads. Primary antibody binding was performed overnight at 4°C in wash buffer supplemented with 0.05% digitonin (Sigma 300410-250MG), 2mM EDTA, 0.1% BSA and 1:50 dilutions of antibody stocks (H3K27ac—Abcam, ab4729-100ug; H3K4me1—Abcam, ab8895; H3K4me3—Abcam, ab8580; H3K27me3—Cell Signalling Technology, 9733; normal rabbit IgG—Cell Signalling Technology, 2729). For secondary antibody binding, a guinea pig anti-rabbit IgG (Antibodies-Online ABIN101961: 2µg) was diluted 1:25 in dig-wash buffer (wash buffer + 0.05% digitonin) and binding performed at 25°C for one hour. pA-Tn5-MEDS-A/B transposome was diluted 1:250 in dig-300 buffer (20mM HEPES pH 7.5, 300mM KCl, 0.5mM spermidine, 0.01% digitonin, and Complete EDTA-Free Protease Inhibitor) and allowed to bind at 25°C for one hour. Feature-specific tagmentation was done in dig-300 buffer supplemented with 10mM MgCl2 at 37°C for one hour and stopped by washing cells with TAPS wash buffer (10mM TAPS pH 8.5, Thermo Scientific J63268.AE + 0.2mM EDTA). DNA was released and crosslinks were reversed by Proteinase K treatment at 55°C for one hour (0.5mg/mL Proteinase K Invitrogen EO0491, 0.5% SDS, 10mM Tris-HCl pH 8.0), followed by clean-up with Zymo DNA Clean and Concentrator (Zymo Research D4013). PCR was carried out using Q5 Hot Start High-Fidelity 2X Master Mix (M0494S) and indexed primers (Supplementary Data 1—Table 2) on a thermocycler (programme: 72°C for 5min, 98°C 30s; 12 cycles of 98°C 10s, 63°C 10s; final extension at 72°C 1min; 4°C hold). PCR products were purified by double-sided size selection with AMPure XP beads (Beckman Coulter A63880, 0.4X followed by 1.3X). QC was performed on a 4200 TapeStation (Agilent) and libraries were quantified using the NEBNext library quantification kit (E7630L). Libraries were sequenced on Illumina’s NextSeq 2000 sequencing system (paired-end 2×151-bp, dual index 8-bp).

6-base-CUT&Tag (6B-C&T)
Preparation of tagmented DNA for 6-base sequencing. Feature-specific tagmentation was performed as in standard CUT&Tag, with modifications. Since the size selection and 6-base enzymatic conversion steps reduce DNA recovery, we increased input material to maintain DNA yields and PCR duplication rates comparable to conventional CUT&Tag. Four CUT&Tag samples (500,000 cells each) were prepared in parallel per biological sample. pA-Tn5-ME2U transposome was diluted 1:25 in dig-300 buffer and bound to cells for one hour at 25°C, followed by tagmentation. Proteinase K digestion was performed at 37°C for one hour and DNA was purified on Zymo DNA Clean and Concentrator. Eluted DNA from all four samples (i.e., 4 × 500,000 cells) was pooled and subjected to two rounds of right-side size selection (0.4X) with SPRI magnetic beads (duet multiomics solution evoC kit, biomodal) and Zymo column purification to remove untagmented genomic DNA. 4pg of long DNA controls (tagmented unmethylated pUC19 and methylated Lambda DNA) were added to each sample at this stage. Gap repair was carried out with 8U Phusion DNA polymerase (M0530S) and 320U Taq DNA ligase (M0208S), in 1X Taq DNA ligase reaction buffer containing 0.8mM dNTPs (NEB N0447S) at 37°C for one hour. DNA was purified on Zymo columns and digested in exonuclease digestion mix (1X NEB buffer 4 NEB B7004S, 40U T7 Exonuclease NEB M0263S, 120U RecJf NEB M0264S, 80U Exonuclease I from E. coli NEB M0293S) for 4h at 37°C, followed by an additional round of Zymo column clean-up.

For optimisation of exonuclease digestion conditions, a range of enzymes were tested to specifically degrade linear dsDNA, in addition to the ones employed in the final workflow (see above). For tests with Exonuclease III (NEB M0206S), 1X rCutSMART (NEB B6004S) was used as the reaction buffer along with 100 U of Exo III per 500,000-cell equivalent of CUT&Tag material. For tests with Exonuclease VIII-truncated (NEB M0545S), 10U were used per 500,000-cell equivalent of CUT&Tag material, also in rCutSMART buffer.

6-base-seq library preparation. The following steps were performed with early access reagents from the duet multiomics solution evoC kit (biomodal). To digest ME2U hairpin termini, 23.7µL of exonuclease-treated DNA (enriched for circular dumbbell-like fragments) were treated with 3µL of HD enzyme (hairpin digestion module) and 3µL HD buffer for 30min at 37°C, followed by 1.8X SPRI clean-up (elution = 14µL water). DNA was denatured at 95°C for 2min and slowly cooled (−0.1°C/s) to 10°C to enable strand separation and hairpin re-folding. At this stage, 0.4pg short hmC control oligonucleotide (spike-in) was added to the sample. Strand copy synthesis (SC module) was setup with 13µL DNA, 2.5µL SC buffer, 2µL dNTP mix, 2µL SC enzyme 1, 2µL SC enzyme 2, 2.5µL hmCP additive and 1µL hmCP enzyme, incubated at 37°C for 30min, followed by denaturation at 95°C for 2min and slow cooling (−0.1°C/s) to 10°C. 2.5µL fork-head adaptor were then added to the DNA sample, along with 0.5µL hairpin ligation additive and 15µL hairpin ligation master mix, and the sample was incubated at 20°C for 15min and cleaned up with 1.2X SPRI beads (elution = 15µL water). After ligation, the methyltransferase (MT) module was used to methylate the copy strand at unprotected 5mCpG sites. DNA from the previous step was mixed with 11.6µL MT buffer, 0.4µL MT additive 1, 1µL MT additive 2 and 13µL freshly reconstituted MT enzyme, incubated at 23°C for one hour, then at 65°C for 15min, and cleaned up with 1.2X SPRI beads (elution = 31µL water). 30µL MT-treated DNA were treated with oxidation (Ox) mix (10µL Ox additive 1 solution, made by reconstituting solid Ox additive 1 in 400µL Ox buffer; 1µL DTT solution; 1µL hmCP additive; 1µL hmCP enzyme and 2µL Ox enzyme), followed by addition of 1µL of 1:1250 dilution of Ox additive 2. Samples were incubated at 37°C for one hour and cleaned up with 1.8X SPRI beads (elution = 31µL water). Deamination (DA module) was performed on 30µL DNA by adding 12.7µL water, 17.5µL DA buffer, 1.8µL ATP, 3.5µL MgCl2, 2µL DA enzyme 1 and 2.5µL DA enzyme 2 and incubating at 37°C for 90min. Samples were purified by adding 6.5µL DA clean-up reagent and concentration with 1.2X SPRI beads (elution = 21µL water). 20µL DNA were then PCR-amplified with 25µL PCR mastermix and 5µL UDI primers (programme: 98°C 30s; 12 cycles of 98°C 10s, 62°C 30s, 65°C 60s; final extension at 65°C 5min; 4°C hold). DNA was purified with 0.9X SPRI beads (elution = 14µL evoC dilution buffer), analysed on a 4200 TapeStation and quantified with the NEBNext library quantification kit. Sequencing was performed on the NextSeq 2000 platform (paired-end 2×151-bp, dual index 8-bp).

6B-C&T workflow without gap repair and exonuclease digestion
We initially developed 6B-C&T using a strategy for end repair of hairpin-containing tagmented DNA which slightly differs from the final optimized workflow (Supplementary Fig. 1). After two rounds of SPRI removal of untagmented material, DNA was end-repaired via strand displacement of the inserted hairpin and gaps were filled in using Phi29 DNA polymerase (NEB M0269S). The reaction (25µL) was setup with SPRI-eluted DNA, 1X Phi29 DNA polymerase reaction buffer, 0.2mM dNTPs, 0.3U/µL Phi29 DNA Pol and 0.1mg/mL albumin, and incubated at 30°C for 30min, then 65°C for 10min, followed by Zymo column clean-up. Eluted DNA was then subjected to USER treatment (HD enzyme from the duet multiomics solution evoC kit) to digest the original ME2U hairpin strands, followed by standard 6-base library construction.

Generation of ground-truth DNA methylation controls
For the 5hmC control, we used 80bp oligonucleotides (ATDbio, UK)26 that were either hydroxymethylated or unmodified at specific CpGs (Supplementary Data 1—Table 2). Here, short synthetic oligonucleotides were first annealed in appropriate pairs and quantified using a dsDNA Qubit Quantification Assay kit (Qiagen). Then, 100ng of the annealed oligonucleotide were ligated to synthetic ME2U hairpin adaptors. The construct was treated with USER enzyme and purified using SPRI beads.

For long genomic DNA controls, 40ng of enzymatically methylated (methylated using a CpG methyltransferase) bacteriophage λ DNA (EpigenDX) were added to an equal proportion of unmethylated pUC19 isolated from a methylation-negative strain of Escherichia coli (Dam–, Dcm–), and quantified using a dsDNA Qubit Quantification Assay kit (Qiagen). This mixture was tagmented using Tn5 preloaded with the ME2U hairpin. The tagmented DNA was purified using SPRI beads.

Whole-genome 6-base-seq using Tn5 tagmentation
Sample preparation and genomic DNA isolation. E14TG2A cells were harvested, spun down and washed once with PBS+0.1% BSA. After pelleting, the supernatant was removed and cells were stored at −20°C until genomic DNA extraction. Genomic DNA isolation was carried out using the Quick-DNA MicroPrep Kit (Zymo Research Cat No D3020) following manufacturer’s instructions. Purified DNA was quantified using Qubit and stored at −20°C.

Tagmentation. To perform whole-genome Tn5-based 6-base-seq, the ME1U hairpin adaptor (Supplementary Data 1—Table 2) was annealed in Nuclease-Free Duplex Buffer (30mM HEPES, pH 7.5; 100mM potassium acetate, IDT 11-01-03-01) at 25µM and then complexed in a 1:1 ratio with Tn5 in TPS buffer (Creative Enzymes Cat No NATE-1629). 12 pmol of loaded Tn5 were used to tagment 80ng of genomic DNA. Tagmentation was carried out in LM Buffer (Creative Enzymes Cat No NATE-1629) for 10min at 56°C. Before DNA purification, DNA was released by adding SDS to a final concentration of 0.05% and incubating at 55°C for 15min. After incubation, DNA was purified with 1.8X SPRI beads and eluted in 16.5µL of nuclease-free water. 1µL of eluted DNA was analysed on a TapeStation to determine the tagmentation efficiency.

6-base-seq library prep and sequencing. 15.5µL tagmented DNA were first gap-filled with 8U of Phi29 DNA polymerase (NEB Cat No M0269) in Phi29 reaction buffer (1X Phi29 buffer, 200µM dNTPs, 100µg/ml albumin) by incubating at 30°C for 30min, followed by denaturation at 65°C for 10min. DNA was then purified with 0.9X SPRI beads and subjected to USER (NEB Cat No 5505) digestion (in 1X rCutSmart buffer, NEB B6004S) to cleave the original ME1U hairpin strand. After USER reaction, DNA was purified with 1.8X SPRI beads, eluted with 12µL library dilution buffer (biomodal) and then denatured at 95°C for 2min followed by slow cooling (0.1°C/s) to 10°C. Resulting DNA was processed using the duet multiomics solution evoC kit (biomodal) according to the protocol described for 6B-C&T and amplified with six PCR cycles. After PCR, samples were purified with 0.9X SPRI beads, eluted in dilution buffer (biomodal), quantified with Qubit and analysed by Tapestation. Sequencing was performed on a NovaSeq 6000 platform (paired-end 2×151-bp, dual index 8-bp).

Bioinformatic analysis and data visualisation
Standard CUT&Tag data analysis. Libraries were demultiplexed using demuxFQ (flags: -c -d -i -e -t1 -r 0.01 -R -l9). FastQC-version 0.11.8 was used to assess the read quality. Bases with a quality score below 20 were trimmed from both reads using Cutadapt (cutadapt -q20). Reads were aligned to the GRCm38/mm10 reference genome using bwa-mem with default parameters. Duplicates were removed using Picard version 2.20.3 (Picard MarkDuplicates). BigWig files were created on the deduplicated BAM file, using deepTools bamCoverage (parameters: --normalizeUsing CPM, --binSize 10). Peaks were called using Seacr version 1.3 without input control reporting the top 1% by AUC regions, using both relaxed and stringent criteria. For consensus peak calling, overlap across biological replicates was calculated with intervene tools generating a series of .bed files. Peaks that were common to at least 2 of the 3 biological replicates were considered as consensus peaks. Metagene profiles were generated using deepTools.

6-base sequencing data pre-processing. FASTQs from 6-base-seq libraries were processed through a modified version of biomodal’s duet analysis pipeline. For 6B-C&T (ME2U hairpin), 9-nt scar sequences, AAGAGATAG (R1) and AAAAAACAA (R2), derived from the deaminated, USER-treated Tn5 mosaic end (ME) sequence, were trimmed from the 5′ ends of reads 1 and 2, respectively. Read pairs that lacked either scar sequence were discarded. Adaptor sequences containing the full 19 nt Tn5 ME sequence and 2 uracils (ME2U) hairpin were trimmed from the 3′ ends of reads 1 and 2. For whole-genome 6-base-seq (ME1U hairpin), NNNNNNNNAGATGTGTATAAGAGATAG (R1) & NNNNNNNNAAATATATATAAAAAACAA (R2) were trimmed from the 5′ ends of reads 1 and 2, respectively. Similarly, read pairs that lacked either scar sequence were discarded. Adaptor sequences containing the full 19 nt Tn5 ME sequence and 1 uracil (ME1U) hairpin were trimmed from the 3′ ends of reads 1 and 2.

In all cases, a further 9 bp were trimmed from the 3′ ends of reads 1 and 2 to remove the gaps generated during Tn5 transposition. Trimmed read pairs were then resolved using the bespoke resolution algorithm from Füllgrabe et al26. The resolved, single-end reads were aligned against a GRCm38/mm10 reference (gencode vM25) that contained additional sequences for pUC19, bacteriophage lambda and short spike-in oligonucleotides using bwa-mem2 (v2.2.1). Multi-mapping reads, secondary and supplementary alignments were filtered out. DNA (hydroxy)methylation was quantified per aligned position from MM tags as part of the duet analysis pipeline. Compiled sequencing data analysis reports for all datasets were obtained on NextFlow and are available in Supplementary Data 1 (Tables 3–7).

Call-rate matrix for 6B-C&T. A call-rate matrix was used (Fig. 1g) to measure the accuracy of 5mC and 5hmC calls. Each cell in the matrix m represents the rate at which the method calls a particular modification state X when the true modification is Y. For example, the cell m(5mC, C) represents the rate at which the method calls unmod C when the true modification status is 5mC. Each row of the matrix corresponds to the rate at which each modification status is called for a particular true modification status. The matrix was estimated using three different spike-in controls: fully unmethylated pUC19 DNA for the row corresponding to a true state of unmethylated C; fully methylated lambda for the row corresponding to a true state of 5mC; and synthetic oligonucleotide for the row corresponding to a true state of 5hmC. For a given row, the rate is calculated as the proportion of bases with each modification status in the set of bases for which the genetic base call is C and which are aligned to CpGs in the given control sequence.

Analysis and visualisation of 6B-C&T and WG 6B-seq data. Individual 5mC and 5hmC percentages were calculated at each CpG site for all replicates. For downstream analyses, only cytosines with a read coverage ≥20 were included, except for genome browser visualisation in IGV, where a minimum coverage of ≥10 was used. BigWig tracks for 6B-C&T–derived 5mC and 5hmC across the four histone marks were generated by merging individual replicates. Merged bedGraph files were created using bedtools (v2.31.0), and BigWig files for IGV visualisation were produced using bedGraphToBigWig.

Average 5mC and 5hmC percentages were calculated for each replicate of the 6B-C&T datasets (across the four histone marks) and the whole-genome 6-base-seq dataset using a custom R script. Differential methylation analyses, including generation of volcano plots, were performed for each 6B-C&T dataset relative to whole-genome 6-base-seq using methylKit (v1.34.0). Window-based analyses were conducted with the methylKit tiling function (win.size= 100; step.size= 100). Cytosines with read coverage <20 were excluded during methylKit processing. Differentially methylated positions were identified using the calculateDiffMeth function. Regions with p<0.05, |log2FC| ≥0.5, and an absolute methylation difference ≥15% between groups were defined as differentially methylated regions (DMRs).

Shared CpG enrichment profiles were generated using a custom bash workflow. The GRCm38/mm10 genome was partitioned into 1kb bins, and average enrichment values for each histone mark were computed across these bins using standard CUT&Tag datasets. For CpGs shared between whole-genome 6-base-seq and 6B-C&T datasets, cytosine modification percentages were plotted as a function of histone mark enrichment (log₂ scale) for each histone mark.

Average read profiles for 6B-C&T centred on each enhancer class were generated using deepTools. Briefly, bamCoverage was used to produce BigWig files from BAM files containing 6B-C&T–resolved reads, using a bin size of 20. Peak density distributions were computed with computeMatrix (100 bp bins) and visualised with plotProfile. The enhancer analysis for global distributions of 5mC and 5hmC was performed using a custom R script (v4.3.1). Briefly, all CpG sites overlapping with enhancer peaks (primed, poised, or active, defined in Cruz-Molina et al.40) were identified with BEDTools and levels of each CpG state (C, 5mC and 5hmC) at each enhancer were calculated and visualised using R and GraphPad Prism.

Publicly available CUT&Tag-BS data for mESCs were downloaded as raw FASTQ files from the European Nucleotide Archive (ENA), EMBL. Adaptor sequences were removed using Cutadapt with the parameters “-a CTGTCTCTTATACAC-ACTGTCTCTTATACAC -O 5 -q 0 -m 20 -p”. The trimmed reads were subsequently aligned to GRCm38/mm10 genome using Bismark v0.24.04446 with parameters “-X 1000 --non_bs_mm”. Publicly available ChIP-seq datasets for wild-type mESCs were downloaded from GEO. Raw reads were filtered to remove those with a mean base quality score <20 and aligned to the mouse genome using BWA. The resulting BAM files were processed, and peak calling was performed with MACS247 with parameters p≤0.01; Fold-enrichment ≥2 and Broad Region Calling ON.

Distributions of 5mC and 5hmC levels across genomic windows (corresponding to IGV snapshots) were visualized using a custom Python workflow that parsed TSV input files and removed invalid entries. Boxplots comparing methylation values across each enhancer mark and IgG were generated with matplotlib, using a fixed colour palette and ordering groups by their median values.

For differential methylation scatter plots, DSS-like estimates of mean CpG methylation for the two treatment groups (6B-C&T and WG 6B-seq), along with dispersion parameters, Wald statistics, and FDR-adjusted q-values, were derived from the calculateDiffMeth output of the methylkit package. These metrics were used to visualise concordance between treatment methylation means (μ₁ vs μ₂) by generating jittered, rasterized scatter plots using ggplot2.

Pearson correlations of modified CpG fractions between 6B-C&T and other libraries were computed for CpGs with ≥20× coverage. Pairwise correlations within samples were obtained using the getCorrelation function from the methylKit package.

Chemical structures in Fig. 1a were generated using ChemDraw®.
IGV48 was used for data visualisation in Fig. 1c, Fig. 2b and Fig. 3a. Plots in Fig. 3b, Fig. 3d, Supplementary Fig. 4 and Supplementary Fig. 16 were generated with Prism (GraphPad Software, Massachusetts USA).

Machine learning analysis
Data preparation. Raw counts were collected for each methylation state (C, 5mC, 5hmC and ambiguous modC, (i.e., modified C sites that cannot be resolved into 5mC or 5hmC) from 6-base-seq data (6B-C&T or WG 6B-seq). Raw counts were then normalised by the total counts at each CpG, and then the average was taken across experimental replicates. CpG sites mapping to enhancer regions were collected, and the sequential CpG methylation states in each enhancer region were organised as a sample. Enhancer annotations were obtained from Cruz-Molina et al.40

Model architecture. As H3K4me1 is the common feature of annotated active, poised, and primed enhancers, we used the samples curated from methylation profiles obtained from H3K4me1 6B-C&T to train a model for classification into these three enhancer classes. To further address the advantages of using H3K4me1-enriched 5mC and 5hmC information, we also trained an independent model with the same architecture based on the methylation profiles generated by whole genome 6B-seq. The model architecture (Supplementary Fig. 17) consists of a linear projection layer that projects the normalised CpG methylation states to a 32-dimension vector (d_model). This projection layer was followed by a self-attention layer, which consisted of four attention heads (num_head) and a two-layer feed-forward network with 256 hidden units (ff_dim), and was used to extract the features for the final classification layer. The classification layer is a two-layer fully connected neural network with 256 hidden units (mlp_h_dim) for each layer. The ‘elu’ activation function was used for non-linearity. A dropout probability of 0.2 (dropp) was set for avoiding overfitting. CpG methylation profiles from chromosome 1 were retained for model evaluation.

Enhancer classification. Enhancer regions were defined and classified into three categories (active, poised, and primed) according to a previous study40. Here, we designed a machine learning model to predict the type of each enhancer based on its 6-base methylation profile. Specifically, for a set of enhancers E = {e1, e2, ..., en}, each enhancer can be represented as a sequential vector of normalised signals {num_ambiguous_modC, mC, hmC, C} for its CpG loci. Specifically, ei = {mi1, mi2, ..., mik}, where mij is the j-th CpG locus in the i-th enhancer, given by mij = {x_ambiguous, x_modC, x_mC, x_hmC, x_C}, and x represents the normalised signal. The enhancer labels can be noted as Y ∈ {0, 1, 2}n, where 0 means active, 1 means poised, and 2 means primed. The classification task is then to find a classifier f : e → y to minimize the expected loss:

f* = argmin E[L(f(e), y)]
f∈F

Where L(·) is a loss function, we used the cross-entropy loss for the classification task.

Self-attention. To capture the associations among CpG methylation states in enhancers, we applied a self-attention module to the model49. After projection to a high-dimension space by the linear projector, the CpG methylation profile is simultaneously used as a query, key, and value matrix to calculate the attention matrix and the weighted value. The weighted value is then converted to informative embedding, which is further used for classification.

Statistics & reproducibility
Experiments were performed independently with different E14TG2A cultures (N=2 for H3K27me3 6B-C&T; N=3 for H3K27ac 6B-C&T; N=3 for H3K4me1 6B-C&T; N=3 for H3K4me3 6B-C&T; N=3 for IgG 6B-C&T and N=2 for WG 6B-seq), and no statistical method was used to predetermine sample size. Pearson correlation coefficient, representing effect size, was used to verify reproducibility. No data were excluded from any analyses. Descriptive statistics were computed using GraphPad Prism 10 (GraphPad) or R. Randomization or blinding were not required nor performed in this study. Differences between groups were assessed in Supplementary Fig. 10 using methylKit, which employs a logistic regression test for p value calculation, and in Supplementary Fig. 14 using a Wilcoxon rank sum test (results in Supplementary Data 1—Table 1).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Epigenetic modifications on DNA and histones act in concert to contribute and reflect changes in chromatin states. In this manuscript, Tavares et al. aimed to integrate the enrichment of histone modifications using the CUT&Tag method with their previously developed 6-letter sequencing, which enables simultaneous detection of 5-methylcytosine (5mC) and 5-hydroxymethylcytosine (5hmC) at single-base resolution.\nWhile detection of 5mC using ChIP or CUT&Tag combined with bisulfite sequencing has been previously reported, the 6-letter sequencing approach represents a novel and important advance. It enables simultaneous detection of both major DNA modifications without bisulfite treatment, thus offering a powerful tool to dissect coordinated transitions between epigenetic states.\nCombining CUT&Tag with 6-letter sequencing is not as straightforward as with bisulfite-based methods. To address this, the authors developed a custom hairpin transposase adapter and optimized the removal of DNA fragments lacking adapters on both ends. When comparing the performance of individual protocols with the combined 6-letter CUT&Tag (6L-C&T), they observed that both enrichment and DNA modification detection were retained at comparable efficiency.\nThe authors then applied 6L-C&T in mouse embryonic stem cells, using antibodies against H3K27ac, H3K4me3, H3K4me1, and H3K27me3. Enrichment revealed changes in DNA methylation and hydroxymethylation when compared to unenriched DNA, consistent with expected chromatin states. Notably, 5mC retention was observed at enhancers marked by H3K4me1, suggesting coordinated regulation at these loci. Furthermore, a machine learning model trained on 6L-C&T data from H3K4me1-marked regions outperformed whole-genome 6-letter sequencing in predicting enhancer states.\nIn summary, 6L-C&T is a novel technique that offers unique advantages for studying the interplay between DNA and histone modifications. It enables simultaneous, base-resolution profiling of two major DNA modifications on the same DNA molecules and benefits from increased sensitivity due to the absence of bisulfite treatment. The manuscript is well written and suitable for publication. I have following minor suggestions:\n1. The authors could discuss how 6L-C&T compares to conventional CUT&Tag in terms of cell input requirements. Does the incorporation of the hairpin adapter and nuclease treatment lead to a significant loss of DNA?\n2. It would be interesting to see locus-specific examples of differential methylation (e.g., at imprinted regions such as H19/Igf2), as these would highlight the method's potential to resolve coupling of allele-specific epigenetic states.", "reviewer_2": "Authors report 6-letter CUT&Tag, a multi-omics method for six-base sequencing (GATC/5mC/5hmC) at antibody targeted chromatin features. Approach used in E14 mESCs to profile the co-occurrence of 5mC and 5hmC at a range of histone PTMs (H3K4me1, H3K4me3, H3K27ac, H3K27me3) and relate to the functional state of enhancer elements (poised, primed, or active). In this example, successful application could offer a more effective means to classify a central transcriptional regulatory element. A potentially valuable approach but extensive revisions would be required to address numerous concerns.\nGeneral comments\nMuch essential background and context are missing, including: how is 5hmC created (e.g., TET enzyme identity, preference and kinetics); any proposed relationships between 5mC and 5hmC (including their relative abundance in cell types that includes mESCs); how this has been measured in prior studies (and how do the authors numbers compare); even the existence of 'third generation' sequencing approaches (ONT, PacBio) and their direct reading of 5mC / 5hmC (whole genome or immunotethered : an extremely active area of study. Such a level of omission on so many essential points is glaring.\nDirect Comments & Questions (in no particular order)\n'Sequencing-based methods have shaped our understanding of the interplay between these epigenetic marks [5mC and 5hmC]” . What are the current thoughts on this [partic on function gained by 5hmC vs. lost when 5mC is oxidised] ? And outstanding questions that could be addressed using 6L-CUT&Tag?\nScar-filtering is an innovative approach : congratulations to the authors\nDatasets exist for native 5mC/5hmC in a variety of sample types (including mESCs). The authors neither mention the alternate approaches used to create same (!), nor make any technical (sensitivity / potential biases) or data comparisons.\nFig 1F: How effective is the '6-letter' approach at distinguishing 5mC / 5hmC on gDNA (e.g., relative efficiency ? Sequence bias (? a 5hmC oligo would not fully inform on same)). How does approach used differ to that in Fig 4c / [PMID: 36747096] (and how has the improvement been achieved) ?\nSuppl Fig S1 is a flowchart of 6L-seq library preparation. This would benefit from a more detailed explanation of the key enzymes / steps. As one example, what is the role of DNMT5 ? (also source of this enzyme? Why was it chosen). It should not be necessary to return to the authors prior work [PMID: 36747096] for this background.\nEssential: CUT&Tag (and by extension 6L-CUT&Tag) requires an IgG background control for best practice. Whole genome 6L provides a very useful landscape (e.g., does heterogeneity exists at a genomic location?) but cannot identify the % signal at each location that is target-antibody mediated vs. approach background (e.g., as in figs 1D, 1E and 2A). Particularly matters here since 'open' locations (enhancers) are where most of the ATAC-like background can appear in a poorly controlled CUT&Tag experiment (of note H3K27me3 is known to 'dampen' this background, likely because of its sheer abundance).\nThe central biology explored in the manuscript is the 5mC/5hmC state at mouse enhancers classed as poised, primed, or active (by Cruz-Molina et al) : “we analysed a set of mESC enhancers that have been annotated by the presence or absence of specific histone marks”. What justifies this annotated set as a gold standard? It was created by ChIP-seq (a noisier method than CUT&Tag). Has it has been independently verified / widely accepted?\n> #187 [& fig 2a whisker plots]: The abundance of 5mC / 5hmC ranges dramatically across the enhancer locations .. which would give this reviewer pause before interpreting biology in differences at the bottom end of the abundance spectrum.\nFig 2a: Each IGV snapshot group should show the distribution of all histone PTMs under study (and an IgG control) : not just those proposed associated with each enhancer class. Windows are too narrow to demonstrate enrichment of H3K4me1/H3K27ac [relative to ? Would need to show flanking locations]. What is driving the scale choice in these locations for %5mC ? [e.g., 1-78 on far left] ?? As below: selected windows are great to make a point but should always be accompanied by whole genome row-linked ChAsE heatmaps).\nThe authors make extravagant use of violin plots to compare cytosine (hydroxy)methyl enrichment between 6L-CUT&Tag datasets. In many cases the large majority of the distribution is at or near zero, with a thin tail trailing upward : what percentage of the data is in the difference ? e.g., Extended Fig 4b what statistical tests were performed to confirm significance between each distribution?\nExtended Data Fig 3a. Authors’ choice to use volcano plots (to summarize the enrichment/depletion of cytosine methylation in 6L-CUT&Tag vs 6L-WGS) is non-intuitive. Could a simpler visualization be a dot plot (e.g., X-axis 5mC enrichment in 6L-WGS / Y-axis 5mC enrichment in 6L-CUT&Tag)? Expectation would be that a steep incline supports a greater enrichment of 5mC associated with the histone PTM (as by 6L-CUT&Tag) vs. general landscape (as by 6L-WGS).\nExtended Data Fig 6: Used to substantiate “that regions annotated as primed enhancers are enriched with the H3K4me1 mark but also show appreciable levels of H3K27ac and H3K27me3.”\nProfile plots seem to support that point, since the y-axes indicate more H3K4me1 than H3K27me3 or H3K27ac at these primed enhancers. Please also provide genomic heatmaps showing the three histone PTMs (row-linked; sorted descending by H3K4me1 signal) along with an IgG control (as above) to identify background signals.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "This manuscript by Tavares et al. is a study of the 5mC and 5hmC status of CpGs at pulled down histone modification sites in the genome. The authors use their previously developed 6-letter bisulfite free sequencing approach with CUT&Tag. The main advance is the use of the new 6-letter bisulfite free approach, over previous versions of this experiment that have been conducted using bisulfite sequencing. While this is an interesting advance, I do have some concerns that will need to be addressed.\n1. There is a lack of data comparison to the previous bisulfite methods, did they see similar things (or not, which is equally interesting).\n2. On page 4, the authors state “6L-C&T-derived 5mC and 5hmC levels at each histone mark were lower than the whole-genome average (Extended Data Fig. 1b)”, which is really interesting.\nHowever, they later say “active chromatin marks (H3K4me3, H3K27ac and H3K4me1) showed similar numbers of CpGs with higher 5hmC levels than in the whole-genome data and CpGs with lower 5hmC levels than in whole genome” and “the repressive mark (H3K27me3) showed the majority of CpG sites with higher 5hmC in 6L-C&T than in whole-genome”\nThese two sets of comments seem to be contradictory. The first that all 5mC/5hmC levels were below WG. Then the second two comments that some 5hmC level were similar to WG and some higher than WG. It didn’t help that the second two sentences didn’t have a figure attached to them.\n3. A general trend throughout is that the 5hmC levels always seem to be a very similar % of the 5mC levels (~25-35%). I’m surprised this isn’t noted and explained. Is this a factor of how the experiments were carried out? Is this identifying something important about 5hmC formation/removal? Is this the same in the WG samples?\n4. I don’t understand the point of the machine learning model that was developed, perhaps I misunderstand it. It seems they put the WG 5mC/5hmC and H3K4me1 5mC/5hmC (enhancer marker) data into the ML system and then show that the H3K4me1 version can predict enhancer regions better. But surely that is obvious, as you have preenriched the enhancer regions already?\n5. Most of the data is not in the main text, I think it would be useful to have a table(s) to summarise the findings as there is a tremendous amount of data that is skimmed over in the text.\n6. In plots Fig 2b, Ext data Fig 4, and Ext data Fig 6b, the violin plots show a very interesting trend. It’s almost as if the %5mC/5hmC is quantised. There are specific bumps of data at regular intervals. The authors should explain why this is, as I’m assuming this cannot be biological so will be a factor of how their experiments were conducts or data processed? What effect does this have on the conclusions?"}, "author_response": "Note from the authors: the method has now been renamed to 6-base-CUT&Tag (6B-C&T) since the term “base” is preferable to “letter” when referring to nucleobases. All instances of the text and figures, including the original reviewer’s comments below were updated for consistency and to avoid confusion.\nReviewer #1 (Remarks to the Author):\nEpigenetic modifications on DNA and histones act in concert to contribute and reflect changes in chromatin states. In this manuscript, Tavares et al. aimed to integrate the enrichment of histone modifications using the CUT&Tag method with their previously developed 6-base sequencing, which enables simultaneous detection of 5-methylcytosine (5mC) and 5-hydroxymethylcytosine (5hmC) at single-base resolution.\nWhile detection of 5mC using ChIP or CUT&Tag combined with bisulfite sequencing has been previously reported, the 6-base sequencing approach represents a novel and important advance. It enables simultaneous detection of both major DNA modifications without bisulfite treatment, thus offering a powerful tool to dissect coordinated transitions between epigenetic states.\nCombining CUT&Tag with 6-base sequencing is not as straightforward as with bisulfite-based methods. To address this, the authors developed a custom hairpin transposase adapter and optimized the removal of DNA fragments lacking adapters on both ends. When comparing the performance of individual protocols with the combined 6-base CUT&Tag (6B-C&T), they observed that both enrichment and DNA modification detection were retained at comparable efficiency.\nThe authors then applied 6B-C&T in mouse embryonic stem cells, using antibodies against H3K27ac, H3K4me3, H3K4me1, and H3K27me3. Enrichment revealed changes in DNA methylation and hydroxymethylation when compared to unenriched DNA, consistent with expected chromatin states. Notably, 5mC retention was observed at enhancers marked by H3K4me1, suggesting coordinated regulation at these loci. Furthermore, a machine learning model trained on 6B-C&T data from H3K4me1-marked regions outperformed whole-genome 6-base sequencing in predicting enhancer states.\nIn summary, 6B-C&T is a novel technique that offers unique advantages for studying the interplay between DNA and histone modifications. It enables simultaneous, base-resolution profiling of two major DNA modifications on the same DNA molecules and benefits from increased sensitivity due to the absence of bisulfite treatment. The manuscript is well written and suitable for publication. I have following minor suggestions:\n1. The authors could discuss how 6B-C&T compares to conventional CUT&Tag in terms of cell input requirements. Does the incorporation of the hairpin adapter and nuclease treatment lead to a significant loss of DNA?\nIn conventional CUT&Tag, theoretically half of all tagmented fragments are not amplified as they have the same adaptor at both ends (e.g., A/A and B/B), whereas fragments with A/B and B/A adaptor combinations are successfully amplified in PCR. In 6B-C&T, insertion of the same hairpin adaptor at either end of the DNA fragment means that in principle all double-sided tagmentation events (both ends of the fragment) should be recovered by PCR.\nIn practice, we have found that most fragments generated in 6B-C&T lack an adaptor at one or both ends. These fragments probably originate from DNA shearing or non-specific nuclease activity throughout the protocol and are efficiently repaired, ligated and amplified in the 6-base-seq workflow. Importantly, we developed a method to minimize the number of these fragments via exonuclease treatment and to computationally exclude them from downstream analysis. The increased number of steps and processing in the 6B-seq protocol does result in lower DNA yield relative to conventional CUT&Tag. However, we compensate for this by scaling up the starting number of cells four-fold, which results in library yields and read duplication rates similar to conventional CUT&Tag.\nThe Methods now include a statement justifying the input requirements of 6B-C&T in comparison to conventional CUT&Tag:\n“Since the size selection and 6-base enzymatic conversion steps reduce DNA recovery, we increased input material to maintain DNA yields and PCR duplication rates comparable to conventional CUT&Tag. Four CUT&Tag samples (500,000 cells each) were prepared in parallel per biological sample.”\n2. It would be interesting to see locus-specific examples of differential methylation (e.g., at imprinted regions such as H19/Igf2), as these would highlight the method's potential to resolve coupling of allele-specific epigenetic states.\nWe examined differential methylation at several genes reported to show allele-specific expression during mouse ESC differentiation (PMID: 30767785, PMID: 34226608). Please see below (Review-only Figure 1) for specific examples.\nThe Peg3/Usp29 bidirectional promoter exemplifies a general epigenetic pattern and harbours an imprinting control region (PMID: 27257070). It is characterized by a CpG-rich domain with global 5mC levels approximately half those of the surrounding regions, as determined by whole-genome 6-base sequencing, and by a pronounced enrichment of the promoter-associated histone mark (H3K4me3). This gene is paternally expressed in ESCs, with the paternal (active) allele associated with H3K4me3 and the maternal (repressed) allele associated with H3K9me3 (PMID: 34226608). Accordingly, %5mC at H3K4me3 (given by H3K4me3 6B-C&T) is ~0, indicating this histone mark is coupled with the unmethylated allele. No coupling could be determined with the other histone marks we profiled (H3K27me3, H3K4me1, H3K27ac), as they are not significantly enriched at this locus. Loci such as Kcnq1ot1 and Igf2r/Airn display similar patterns. H19 does not show appreciable levels of the histone marks profiled in our study, and we did not detect differential methylation at Igf2 in our data.\n6B-C&T is ideally suited to detect allele-specific coupling of histone marks with 5mC and 5hmC. We anticipate that exploratory analysis of histone modifications such as H3K9me3 and H3K20me3 will be highly informative in this context, since these marks are strongly associated with differential methylation at imprinted loci (PMID: 20180964).\nReviewer #2 (Remarks to the Author):\nAuthors report 6-base CUT&Tag, a multi-omics method for six-base sequencing (GATC/5mC/5hmC) at antibody targeted chromatin features. Approach used in E14 mESCs to profile the co-occurrence of 5mC and 5hmC at a range of histone PTMs (H3K4me1, H3K4me3, H3K27ac, H3K27me3) and relate to the functional state of enhancer elements (poised, primed, or active). In this example, successful application could offer a more effective means to classify a central transcriptional regulatory element. A potentially valuable approach but extensive revisions would be required to address numerous concerns.\nGeneral comments\n1. Much essential background and context are missing, including: how is 5hmC created (e.g., TET enzyme identity, preference and kinetics); any proposed relationships between 5mC and 5hmC (including their relative abundance in cell types that includes mESCs); how this has been measured in prior studies (and how do the authors numbers compare); even the existence of 'third generation' sequencing approaches (ONT, PacBio) and their direct reading of 5mC / 5hmC (whole genome or immunotethered : an extremely active area of study. Such a level of omission on so many essential points is glaring.\nAdditional background has been added (see pages 1-2 – Intro paragraphs 1-2), including the relative 5mC/5hmC abundance in mESCs:\n“Chemical modifications of DNA bases and histone tails regulate physical and functional properties of chromatin, including DNA accessibility, genomic architecture and binding of transcription regulators1-4. Methylation of cytosine bases on DNA is critical to establish the epigenetic landscape and governs gene activity and inheritance5. Cytosine bases can be methylated by DNA methyltransferases (DNMTs), that transfer a methyl group from the donor S-adenosylmethionine to form 5-methylcytosine (5mC)6, which can be further oxidised to form 5-hydroxymethylcytosine (5hmC) through the action of the Ten-Eleven Translocation family of proteins (TETs 1-3)7. TET enzymes can further oxidise 5hmC to other intermediates like 5-formylcytosine (5fC) and 5-carboxylcytosine (5caC)8, each of which can be excised by thymine DNA glycosylase, followed by restoration of C, completing a demethylation cycle (Figure 1). 5mC is the most abundant cytosine modification (~60% of all Cs in ESCs)9 and 5hmC is the next most abundant cytosine modification (less than 5% of all Cs in ESCs)9, with 5fC and 5caC being present at substantially lower levels (~0.002% and ~0.0003%, respectively)8. While 5hmC can be an intermediate in active DNA demethylation, it has also been shown to be a stable epigenetic mark that is associated with transcriptional activity and cell fate10-15. The functional roles of 5mC and 5hmC depend on their genomic context. Broadly, 5mC enrichment at promoters and enhancers is associated with transcriptional repression16, 17, whereas 5hmC is enriched at active regulatory regions – particularly enhancers – and also within the gene bodies of transcriptionally active genes18, 19”\n“The development of sequencing methods to detect 5mC and 5hmC has been an important focus in epigenomics. While early methods such as bisulfite sequencing20 cannot distinguish 5hmC from 5mC, later approaches – including oxidative bisulfite sequencing21, TET-assisted bisulfite sequencing22, APOBEC-coupled epigenetic sequencing23 and third-generation sequencing platforms (e.g., PacBio, Oxford Nanopore)24 – have enabled base-resolution mapping of both 5mC and 5hmC. More recent methods have now leveraged state-of-the-art enzymology and chemistry to enable bisulfite-free and simultaneous detection of 5mC and 5hmC on the same DNA fragment25-27.”\nOur data is consistent with previous reports showing the relative genome-wide abundance of 5mC corresponds to ~10-30X that of 5hmC. This information is now included in the text (page 3, Results paragraph 4):\n“As a reference, we also performed whole-genome 6-base-seq, which provides an untargeted picture representative of the total pool of DNA fragments. Using this approach, we confirmed that the relative 5mC/5hmC abundance was consistent with reported values in mESCs (Extended Data Fig. 1a).”\nDirect Comments & Questions (in no particular order)\n2. 'Sequencing-based methods have shaped our understanding of the interplay between these epigenetic marks [5mC and 5hmC]” . What are the current thoughts on this [partic on function gained by 5hmC vs. lost when 5mC is oxidised] ? And outstanding questions that could be addressed using 6B-CUT&Tag?\nSee response to Referee 2 point 1 above together with the following introductory text (page 2, Intro paragraph 2):\n“Despite such advancements, a key challenge in the field has been to determine how 5mC and 5hmC are coupled in space and time to other chromatin features such as histone marks, transcription factors and chromatin remodelers. Establishing this relationship is critical to appreciate how epigenetic features are coupled, when and where they co-occur, and how they determine causal relationships and regulatory roles. Current strategies that map DNA modifications and chromatin features largely rely on comparing profiling data from independent experiments for different features. This approach does not necessarily determine co-occurrence of features at the molecular level (Fig. 1b).”\nScar-filtering is an innovative approach : congratulations to the authors\n3. Datasets exist for native 5mC/5hmC in a variety of sample types (including mESCs). The authors neither mention the alternate approaches used to create same (!), nor make any technical (sensitivity / potential biases) or data comparisons.\n3.1 With respect the citation of approaches to measure 5hmC/5mC see response to Referee 2 point 1.\n3.2 To our knowledge, 6B-C&T is the first method to directly resolve 5mC and 5hmC bases on histone mark-containing DNA fragments. Previous 6-base sequencing methods have so far been performed solely on purified DNA for whole genome samples (e.g., PMID: 38846078, PMID: 38336903). To understand the main technical differences between the whole-genome vs targeted approaches, we have also performed extensive comparisons between whole-genome 6-base-seq and 6B-C&T throughout the manuscript (Extended Data Figs. 1-3).\n3.3 Our findings agree with reports that used other methods for simultaneous detection of DNA methylation at histone marks like ChIP-BS (PMID: 33844685), CUT&Tag-BS (PMID: 35028637) and scEpi-seq (PMID: 40999097). Specifically, we observe similar levels of total DNA methylation at the histone marks profiled in common with those studies, such as H3K4me1 and H3K27me3. In addition, our detection of different histone mark-associated fragments mapping to primed enhancers confirms and expands earlier characterisations of epigenetic heterogeneity at enhancer loci (PMID: 33844685).\nThe revised Discussion highlights the consistency between these previous reports and our study (page 5, Discussion paragraph 1):\n“6B-C&T detects 5mC and 5hmC bases that co-occur with target chromatin features. Profiling of major regulatory histone modifications, such as H3K4me1 and H3K27me3, reveals patterns that are consistent with other 5\nResponse to reviewers\nmeasurements of total DNA methylation at these marks31, 39, 40 and adds 5hmC as a new dimension to the multiomics analysis of individual DNA fragments.”\nWe have also directly compared our 6B-C&T data to publicly available data profiling total DNA methylation (modC) at H3K4me1 and H3K27me3 via CUT&Tag-BS (PMID: 35028637). We observed good agreement between the two methods both at the level of CUT&Tag fragment enrichments and at the level of modC detection (Supplementary Figure 6a-b, shown below). We also observed excellent agreement of total DNA methylation data from both methods at the enhancer regions studied throughout our manuscript (Supplementary Figure 6c).\nResults have been updated to include the above comparisons (page 3, Results paragraph 4; page 4, Results paragraph 7):\n“…we assessed how well 6B-C&T correlates with CUT&Tag-BS, a method that profiles total cytosine methylation (modC) on CUT&Tag DNA fragments via bisulfite conversion. Both genomic feature enrichment and CpG modification levels for 6B-C&T agreed well with publicly available CUT&Tag-BS data for H3K4me1 and H3K27me3 (Pearson r ≥ 0.83, Supplementary Fig. 6).”\n“Excellent agreement was also obtained for total DNA methylation levels (%modC) at all enhancer types between 6B-C&T and available C&T-BS data for the general enhancer mark, H3K4me1 (Supplementary Fig. 6c).”\nThis confirms 6B-C&T can retrieve published patterns of simultaneous histone mark and DNA methylation accumulation in addition to deconvolving 5mC from 5hmC on the same DNA fragment.\na\nb\n6B-C&T\nH3K4me1 CnT-BS\nrep 1 H3K27me3\n6B-C&T\nH3K4me1 rep 2\n6B-C&T\nH3K27me3 rep 1\n6B-C&T\nH3K27me3 rep 2\nc\nSupplementary Figure 6. Comparison between 6-base-CUT&Tag (this study) and available CUT&Tag-BS (GSE179266) data. (a) Pearson correlation matrices of genomic fragment enrichments for 6B-C&T compared to a publicly available C&T-BS dataset for two histone marks: H3K4me1 (left) and H3K27me3 (right). (b) Correlation matrices for genome-wide total DNA methylation (modC) fractions for the same datasets (left, H3K4me1; right, H3K27me3). Histograms of modified CpG fractions for 200-bp genomic tiles from each available experimental replicate are shown on the diagonal alongside scatter plots below the diagonal. Pearson correlation coefficients are displayed for each pairwise comparison above the diagonal. (c) Comparison of total DNA methylation (modC) from H3K4me1 6B-C&T and H3K4me1C&T-BS at enhancers. Stacked bar plots show % modC (green) and % C (grey) for different sets of annotated enhancer regions. Error bars represent standard deviation (3 independent experiments shown for 6B-C&T; 2 independent experiments shown for C&T-BS).\n7\n8.0\n4.0\n0.0\n0.0 0.4 0.8\n8.0\n4.0\n0.0\n0.0 0.4 0.8\n8.0\n4.0\n0.0\n8.0 4.0\n0.0\n0.0 0.4 0.8\n8.0 4.0\n0.0\n00..00 00..44 00..88\n8.0\n4.0\n0.0\n0.0 0.4 0.8\n0.94\n0.0 0.4 0.8\n8.0\n4.0\n0.0\n0.0 0.4 0.8\n8.0 4.0\n0.0\n0.0 0.4 0.8\n8.0\n4.0\n0.0\n0.93\n0.93\n6B-C&T H3K4me1\nrep 3\n0.0 0.4 0.8\n8.0 4.0\n0.0\n00..00 00..44 00..88\n8.0\n4.0\n0.0\n0.0 0.4 0.8 0.0 0.2 0.4 0.6 0.8 1.0\n0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.00.0 0.2 0.4 0.6 0.8 1.0\n0.1\n8.0\n6.0\n4.0\n2.0\n0.00.1\n8.0\n6.0\n4.0\n2.0\n0.0\n0.1\n8.0\n6.0\n4.0\n2.0\n0.0\n0.1\n8.0\n6.0\n4.0\n2.0\n0.0\n0.83\n0.83\n0.83\nCnT-BS H3K4me1 rep 1\nCnT-BS H3K4me1\nrep 2\n0.0 0.4 0.8\n8.0\n4.0\n0.0\n8.0\n4.0\n0.0\n0.86\n0.86\n8.0\n4.0\n0.0\n0.85\n0.93\n0.0 0.4 0.8\n8.0\n4.0\n0.0\nModified CpG Pearson correlation Modified CpG Pearson correlation\n0.89\n0.0 0.2 0.4 0.6 0.8 1.0 0.1\n8.0\n6.0\n4.0\n2.0\n0.0\nCnT-BS 6B-CnT CnT-BS 6B-CnT\nH3K4me1 H3K4me1 H3K27me3 H3K27me3\nRep 1Rep 2Rep 1Rep 2Rep 3\nRep 1 Rep 1 Rep 2\n0.84 0.84 1.00 1.00 1.00 6B-CnT H3K4me1 Rep 3\n0.91 0.99 1.00 6B-CnT H3K27me3 Rep 2\n0.84 0.84 1.00 1.00 1.00 6B-CnT H3K4me1 Rep 2\n0.83 0.83 1.00 1.00 1.00 6B-CnT H3K4me1 Rep 1 0.90 1.00 0.99 6B-CnT H3K27me3 Rep 1\n1.00 1.00 0.83 0.84 0.84 CnT-BS H3K4me1 Rep 2\n1.00 0.90 0.91 CnT-BS H3K27me3\n1.00 1.00 0.83 0.84 0.84 CnT-BS H3K4me1 Rep 1\n0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0\n0.89\n0.91\nH3K4me1 6B-CUT&Tag vs CUT&Tag-BS at enhancers\nmodC\nC\nnoitalyhtem\nGpC\n%\n100\n80\n60\n40\n20\n0\n6B-C&T active C&T-BS active 6B-C&T primed C&T-BS primed 6B-C&T poised C&T-BS poised\nResponse to reviewers\n4. Fig 1F: How effective is the '6-base' approach at distinguishing 5mC / 5hmC on gDNA (e.g., relative efficiency ? Sequence bias (? a 5hmC oligo would not fully inform on same)). How does approach used differ to that in Fig 4c / [PMID: 36747096] (and how has the improvement been achieved) ?\nFullgrabe et al. (PMID: 36747096) have described how the 6-base method has comparable/superior detection sensitivity and specificity when measured against other platforms, including bisulfite-based and enzymatic methods. It is important to note that 6-base-seq can only distinguish 5hmC from 5mC at CpG sites. Methylation outside of the CpG sequence context, such as CHH, is detected but cannot be deconvolved due to the method’s enzymology. Given the majority of (hydroxy)methylation in mESCs is present at CpG sites (PMID: 25042786), this was not a limitation in our study.\nWhile the ground-truth spike-in controls provide an accurate estimate of the method’s detection sensitivity and specificity, we also assessed whether 5hmC detection exhibits significant sequence preferences beyond the CpG context. Using chromosome 1 as an example (see Review-only Figure 2 below), we performed motif enrichment analysis using HOMER on 5hmC sites identified in 6-base-seq. The analysis revealed multiple distinct sequence motifs present at genomic 5hmC sites, which indicates that 5hmC is detected across diverse genomic contexts with no strong representation of a specific sequence bias other than the CpG context.\nRegarding how our 6-base-seq approach compares to Fullgrabe et al.: the main difference is that the previous study used sonication to shear gDNA, followed by end repair and hairpin ligation (3 steps), whereas we use Tn5 tagmentation to 8\nResponse to reviewers\nsimultaneously fragment and insert the hairpin adaptor (1 step). In our approach, we observe improved efficiency of 6-base-seq library generation, which is accompanied by slightly improved 5mC and 5hmC detection (Fig. 2d) when compared to the results shown in Fig 4c of Fullgrabe et al. (PMID: 36747096).\n5. Suppl Fig S1 is a flowchart of 6B-seq library preparation. This would benefit from a more detailed explanation of the key enzymes / steps. As one example, what is the role of DNMT5 ? (also source of this enzyme? Why was it chosen). It should not be necessary to return to the authors prior work [PMID: 36747096] for this background.\nWe have now explained the workflow in detail (Supplementary Figure 1a, caption pasted below for convenience). We also reference the original 6-base-seq paper (Fullgrabe et al. 2023, PMID: 36747096), where extensive details are provided, including the role of DNMT5, chosen for its high specificity for copy methylation over de novo methylation. Details of all reagents used for 6-base-seq library preparation, which are commercially available, are also provided in the Methods of the present manuscript.\n“Supplementary Fig. 1. Detailed 6-base-seq library preparation from tagmented DNA. (a) Workflow including DNA circularization and exonuclease treatment. For clarity, subsequent enzymatic steps are shown only for the bottom strand of the initial dsDNA fragment. DNA is first circularised by gap filling with Phusion DNA polymerase and nicks sealed with Taq DNA ligase. Subsequent exonuclease treatment specifically digests linear DNA while dumbbell-like fragments are protected. (Uracil-specific excision reagent) USER cleavage then generates a single nucleotide gap at both uracil sites present in each hairpin. This allows for strand separation of the original duplex via heat denaturation and hairpin re-folding on each strand via slow cooling. Klenow (3’ à 5’ exo–) next synthesises a copy strand on each individual strand. T4 polynucleotide kinase (PNK) is included for gap repair of USER-digested sites. A “fork-head” adaptor is then ligated to the free end of double-strand DNA fragments. DNMT5, chosen for its high specificity of copy methylation over de novo methylation, then specifically methylates the CpG on the recently synthesised copy strand only when the original strand contains a 5mCpG1. TET2 is next used to oxidise all 5mC groups to 5hmC followed by their protection in the same reaction by glycosylation using beta-glucosyltransferase (bGT). Finally, APOBEC3A (ssDNA deaminase) converts all unprotected cytosines to uracils, assisted by the dsDNA unwinding activity of the UvrD helicase. Primers complementary to the fork-head adaptor sequences and containing Illumina P5 and P7 sequences are then used to amplify all fragments by PCR to generate a library for Illumina sequencing. All cytosine states are schematically represented and colour-coded according to the legend. “5mC” = 5-methylcytosine; “5hmC” = 5-hydroxymethylcytosine; “5ghmC” = 5-glucosylhydroxymethylcytosine.”\n6. Essential: CUT&Tag (and by extension 6B-CUT&Tag) requires an IgG background control for best practice. Whole genome 6B provides a very useful landscape (e.g., does heterogeneity exists at a genomic location?) but cannot identify the % signal at each location that is target-antibody mediated vs. approach background (e.g., as in figs 1D, 1E and 2A). Particularly matters here since 'open' locations (enhancers) are where most of the ATAC-like background can appear in a poorly controlled CUT&Tag experiment (of note H3K27me3 is known to 'dampen' this background, likely because of its sheer abundance).\nWe have now performed 6B-C&T on E14TG2A mESCs using a non-specific IgG control to assess the amount of accessibility-driven background that is unrelated to target recognition by a specific antibody.\nTo estimate the contribution of non-specific background to the total signal detected with 6B-C&T genome-wide, we compared the percentage of USER-scarred read pairs for histone mark-targeted libraries, which are a proxy for total signal, and IgG control libraries, which serve as a proxy for the background:\nN = number of reads with scars on both ends (bona fide tagmentation events) and\ns\nN = number of reads with no scars or with single-end scars (random fragmentation, ns\nDNA shearing, partial tagmentation).\nF (fraction of scarred pairs) = !! = $ scars !! # !#! $# $#! $! à !! = %!%&'! = (scarred/unscarred fragment ratio).\n!#! ($'%!%&'!) \nWe can use the scarred/unscarred fragment ratio to estimate the contribution of non-specific, accessibility-driven background (given by IgG) to the total signal observed in the histone mark-targeted experiment:\n$! 5!%&'!()*+ %,#-',.) $#!()*+ %,#-',.) = [785!%&'!()*+ %,#-',.)] $! 5!%&'!(01!-,#2 3&'4) $#!(01!-,#2 3&'4) [785!%&'!(01!-,#2 3&'4)] \nOn average, F for IgG libraries = 1.6% and F for histone mark libraries ~ 55% scars scars (range = 50-60%):\n$! :.=:7< $#!()*+ %,#-',.) = [78:.:7<] = 0.013, or 1.3%\n$! :.== $#!(01!-,#2 3&'4) [78:.==] \nFrom the above figure, the IgG-estimated background represents, on average, ~1-2% of the total signal detected by histone mark 6B-C&T, which confirms the high specificity of this method.\nAlso see responses to points 9 and 12, which illustrate the above at specific loci and with genome-wide heatmaps.\n7. The central biology explored in the manuscript is the 5mC/5hmC state at mouse enhancers classed as poised, primed, or active (by Cruz-Molina et al) : “we analysed a set of mESC enhancers that have been annotated by the presence or absence of specific histone marks”. What justifies this annotated 10\nResponse to reviewers\nset as a gold standard? It was created by ChIP-seq (a noisier method than CUT&Tag). Has it has been independently verified / widely accepted?\nWe chose this annotation set as it was comprehensive and derived independently from CUT&Tag. This dataset covered the same histone modifications that we investigated using the same E14TG2A mESC cells and importantly the same culture conditions (serum/LIF medium) and so represents the same cell state of pluripotency/differentiation. Furthermore, a key feature in this dataset is the use of P300 occupancy as a criterion for enhancer type assignment, which enables the active enhancers to be confidently distinguished from other classes. This dataset has been widely used as a reference in other studies of mESCs (see PMID: 33620319, PMID: 32115407, PMID: 37456851, PMID: 40743391, PMID: 34792172). While we could have annotated enhancers de novo using our CUT&Tag data, we felt that it was essential to use comparative data independent of Tn5-based methods and their associated sequence/regional biases. This choice enabled us to cross-validate our genomic assignments against previously published, orthogonal data.\nSupporting this choice of enhancer annotation, we observed good agreement between 6B-C&T and the ChIP-seq data from Cruz-Molina et al. for all enhancer-associated histone marks at the reference loci used in our study (Supplementary Fig. 9, see below). This is now mentioned in the Results (page 4, Results paragraph 7) and rules out the noisier character of ChIP-seq as a potential confounder in the present study.\n“We observed good agreement at these loci between 6B-C&T and the ChIP-seq data originally used for their annotation (Supplementary Fig. 9).”\na b c\n6000\n5000\n4000\n3000\n2000\n1000\n0\n0 250 500 750 1000\nH3K27ac 6B-C&T H3K4me1 6B-C&T H3K27me3 6B-C&T\nH3K4me1 6B-C&T H3K4me1 6B-C&T\nd\nSupplementary Fig. 9. Comparison between 6-base-CUT&Tag (this study) and ChIP-seq data used for enhancer class annotation (GSE89211). Correlation scatter plots for ChIP-seq data from Cruz-Molina et al. (2017) and 6B-C&T data for enhancer-associated histone marks are plotted for different sets of annotated enhancer regions: (a) active enhancers (H3K27ac and H3K4me1), (b) primed enhancers (H3K4me1) and (c) poised enhancers (H3K27me3 and H3K4me1). In all cases, CPM-normalised genomic enrichments (bin size = 10 bp) are plotted for each experiment (y-axis = ChIP-seq and x-axis = 6B-C&T) and R-squared (R2) values are annotated on each plot. Least squares polynomial (first-degree) fitting was performed with numpy.polyfit to produce the fitting line. A single experimental replicate from each experiment was used as representative of the overall trend. (d) Pearson correlation matrix for data plotted in a-c. Pearson correlation coefficients (Pearson’s R) are shown for each pairwise comparison.\n12\n)aniloM-zurC(\nqes-PIhC\nca72K3H(\n)aniloM-zurC(\nqes-PIhC\n1em4K3H(\n)aniloM-zurC(\nqes-PIhC\n1em4K3H(\n)aniloM-zurC(\nqes-PIhC\n3em72K3H(\n)aniloM-zurC(\nqes-PIhC\n1em4K3H(\nActive enhancers Primed enhancers Poised enhancers\n300\n800\n250\n600 200\n400 150\n100\n200\n50\n0 0\n0 500 1000 0 100 200 300\n175\n150\n400\n125\n300 100\n75\n200\n50\n100\n25\n0 0\n0 500 1000\n0 100 200 300\nH3 K4 me H 1 a 3 c C K ti h 2 v I 7 e P a - H c s a 3 e C c K q t h i 4 v I m P e - e H s a 1 e 3 c q 6 t K i B v 2 e - 7 C a & a c H c T 3 6 t K i B v 4 - e C m & e p H T 1 r 3 i m 6 K B 4 e - m d C H e & p 1 3 T r K C i m 4 h m e I P d e p - H 1 s o 3 e i 6 s K q B e 4 - d m C & e H 1 p T 3 o K C is 2 h e 7 I P d m - p e s H e o 3 q i 3 s 6 K e B 2 d - 7 C m & e T p 3 o C is h e I d P-seq\n-0.05-0.04-0.05-0.05-0.06-0.08 0.68 0.75 0.79 1.00 H3K27me3 ChIP-seq poised (Cruz-Molina)\n-0.05-0.04-0.05-0.05-0.06-0.08 0.75 0.68 1.00 0.79 H3K27me3 6B-C&T poised\n-0.05-0.04-0.05-0.05-0.07-0.09 0.84 1.00 0.68 0.75 H3K4me1 ChIP-seq poised (Cruz-Molina)\n-0.05-0.04-0.05-0.05-0.06-0.08 1.00 0.84 0.75 0.68 H3K4me1 6B-C&T poised\n-0.26-0.20-0.26-0.27 0.87 1.00 -0.08-0.09-0.08-0.08 H3K4me1 ChIP-seq primed (Cruz-Molina)\n-0.20-0.16-0.20-0.21 1.00 0.87 -0.06-0.07-0.06-0.06 H3K4me1 6B-C&T primed\n0.77 0.85 0.92 1.00 -0.21-0.27-0.05-0.05-0.05-0.05 H3K27ac 6B-C&T active\n0.87 0.81 1.00 0.92 -0.20-0.26-0.05-0.05-0.05-0.05 H3K4me1 6B-C&T active\n0.73 1.00 0.81 0.85 -0.16-0.20-0.04-0.04-0.04-0.04 H3K27ac ChIP-seq active (Cruz-Molina)\n1.00 0.73 0.87 0.77 -0.20-0.26-0.05-0.05-0.05-0.05 H3K4me1 ChIP-seq active (Cruz-Molina)\n−1.00−0.75−0.50−0.25 0.00 0.25 0.50 0.75 1.00\nResponse to reviewers\n8. > #187 [& fig 2a whisker plots]: The abundance of 5mC / 5hmC ranges dramatically across the enhancer locations .. which would give this reviewer pause before interpreting biology in differences at the bottom end of the abundance spectrum.\nWe observed significantly higher 5mC and 5hmC levels at primed enhancers compared to active and poised enhancers. We note that within each enhancer type, 5mC and 5hmC levels can span a wide range. Coupled with intrinsic heterogeneity across the cell population, lower abundances of these marks will be harder to interpret biologically without other experimental data. We have therefore focused our main analysis on the higher levels of histone mark-specific 5mC and 5hmC found at primed enhancers.\n9. Fig 2a: Each IGV snapshot group should show the distribution of all histone PTMs under study (and an IgG control) : not just those proposed associated with each enhancer class. Windows are too narrow to demonstrate enrichment of H3K4me1/H3K27ac [relative to ? Would need to show flanking locations]. What is driving the scale choice in these locations for %5mC ? [e.g., 1-78 on far left] ?? As below: selected windows are great to make a point but should always be accompanied by whole genome row-linked ChAsE heatmaps).\nEach IGV snapshot now includes all enhancer histone marks along with the IgG control (see below). Windows were expanded to include flanking regions and the annotated enhancer under study is now indicated in each panel. We note that the poised enhancer example (Hoxa11 locus) is part of a bigger Hox gene cluster (~120 kb) which is entirely covered by H3K27me3 signal in E14TG2A mESCs. Further expanding the window in this case to capture the H3K27me3 signal drop would greatly reduce the amount of detail and interpretability of the snapshot.\nThe y-axis scales, previously presented as per IGV default, are now adjusted to fit the data range in each window and all tracks within each group now have the same scale (I – CUT&Tag enrichment; II – %5mC and III – %5hmC):\na\nActive enhancer Primed enhancer Poised enhancer\n(Cyrib locus) (Pak2 locus) (Hoxa11 locus)\nchr15:64,020,221-64,027,686 chr16:32,044,758-32,098,464 chr6:52,239,551-52,250,556\n7,467 bp 53 kb 11 kb\nH3K4me1[0 - 2.03] H3K4me1[0 - 1.3] H3K4me1 [0 - 4.75]\nH3K27ac[0 - 2.03] H3K27ac [0 - 1.3] H3K27ac [0 - 4.75]\nH3K27me3[0 - 2.03] H3K27me3[0 - 1.3] H3K27me3 [0 - 4.75]\nIgG[0 - 2.03] IgG[[00 -- 11..33]] IgG [0 - 4.75]\nH3K4me1 5mC%[0 - 59] H3K4me1 5mC%[0 - 50] H3K4me1 5mC% [0 - 60]\nH3K27ac 5mC%[0 - 59] H3K27ac 5mC%[0 - 50] H3K27ac 5mC% [0 - 60]\nH3K27me3 5mC%[0 - 59] H3K27me3 5mC%[0 - 50] H3K27me3 5mC% [0 - 60]\nIgG 5mC%[0 - 59] IgG 5mC%[0 - 50] IgG 5mC% [0 - 60]\nH3K4me1 5hmC%[0 - 24] H3K4me1 5hmC%[0 - 26] H3K4me1 5hmC% [0 - 16]\nH3K27ac 5hmC%[0 - 24] H3K27ac 5hmC%[0 - 26] H3K27ac 5hmC% [0 - 16]\nH3K27me3 5hmC%[0 - 24] H3K27me3 5hmC%[0 - 26] H3K27me3 5hmC% [0 - 16]\nIgG 5hmC%[0 - 24] IgG 5hmC%[0 - 26] IgG 5hmC% [0 - 16]\nEnhancer Enhancer Enhancer\nFigure 3a. Genome browser views of representative active, primed and poised enhancer loci.\nTracks show histone mark levels and percentages of 5mC and 5hmC obtained with 6B-C&T\nfor enhancer-associated histone marks, in addition to IgG control. The genomic location of each annotated enhancer is indicated with a horizonal bar and shaded area coloured according to enhancer type (active – blue, primed – green and poised – red). Tracks within each group (CUT&Tag fragment enrichments, % 5mC and % 5hmC) have the same y-axis scale.\nBoxplot quantification of the full windows on IGV snapshots above and of the annotated enhancers are now shown in Supplementary Figure 10:\na\n60\n50\n40\n30\n20\n10\nH3K27me3H3K27ac H3K4me1 IgG H3K27me3H3K27ac H3K4me1 IgG H3K27me3 H3K27ac H3K4me1 IgG\nb\nH3K27me3 H3K27ac H3K4me1 IgG H3K27me3 H3K27ac H3K4me1 IgG H3K27me3 H3K27ac H3K4me1 IgG\nH3K27me3 H3K27ac H3K4me1 IgG H3K27me3 H3K27ac H3K4me1 IgG\nSupplementary Fig. 10 Boxplots of %5mC and %5hmC in genomic windows (IGV\nsnapshots) shown in Main Figure 3. (a) Boxplots of %5mC from 6B-C&T at enhancer marks (H3K27me3 in red, H3K27ac in blue and H3K4me1 in green) and IgG control (black) for active (left column), primed (centre column) and poised (right column) enhancer loci. The top panels show the full genomic windows in Main Figure 3 and the bottom panels show only the\n%\nCm5\n%\nCmh5\n%\nCmh5\n%\nCm5\n%\nCm5\n%\nCmh5\n%\nCm5\n%\nCm5\n%\nCmh5\n%\nCmh5\n%\nActive enhancer window (Cyrib locus) Primed enhancer window (Pak2 locus) Poised enhancer window (Hoxa11 locus)\nchr15:64,020,221-64,027,686 (7,467 bp) chr16:32,044,758-32,098,464 (53 kb) chr6:52,239,551-52,250,556 (11 kb)\n30\n25\n20\n15\n10\n5\nchr15:64,023,886-64,024,353 (467 bp) chr16:32,066,114-32,073,757(7,643 bp) chr6:52,242,840-52,243,473 (633 bp)\n4 40\n3 30\n2 20\n1 10\n0 0\nActive enhancer window (Cyrib locus) Primed enhancer window (Pak2 locus) Poised enhancer window (Hoxa11 locus)\nchr15:64,020,221-64,027,686 (7,467 bp) chr16:32,044,758-32,098,464 (53 kb) chr6:52,239,551-52,250,556 (11 kb)\nAnnotated active enhancer Annotated poised enhancer\nchr15:64,023,886-64,024,353 (467 bp) chr6:52,242,840-52,243,473 (633 bp)\nH3K27me3H3K27ac H3K4me1 IgG\nCmh5\n%\n20\n15\n10\n5\n0\n14\n12\n10\n8\n6\n4\n2\n0\n14 10\n20 12 8\n10\n15 6\n8\n10 6 4\n4\n5 2\n2\n0 0 0\nAnnotated primed enhancer\nchr16:32,066,114-32,073,757(7,643 bp)\n25 5\n4\n20 4\n3 15 3\n2 10 2\n1 5 1\n0 0 0\nResponse to reviewers\nannotated enhancer window in each case. (b) Same as in a but for %5hmC. Whiskers represent Q3 + 1.5 x interquartile range (upper) and Q1 – 1.5 x interquartile range (lower).\nFor genome-wide heatmaps, please see response to point 12.\n10. The authors make extravagant use of violin plots to compare cytosine (hydroxy)methyl enrichment between 6B-CUT&Tag datasets. In many cases the large majority of the distribution is at or near zero, with a thin tail trailing upward : what percentage of the data is in the difference ? e.g., Extended Fig 4b what statistical tests were performed to confirm significance between each distribution?\nWe have now added boxplots to Extended Data Figure 4 to indicate the median, interquartile range, lower and upper bounds of each distribution in the violin plots (see below). In all cases, 80-95% of all data points are found within the distribution (< 3rd quartile + 1.5X interquartile range). The following were added to Supplementary Table 1: Wilcoxon rank sum test p-values that confirm significance for each comparison; total number of CpGs analysed; upper bound values for each distribution and the percentage of CpGs outside of the distribution (outliers).\na\nb\nc\n30\n20\n10\n16\nca72K3H\nta Cm5\n%\n3em72K3H\nta Cm5\n%\nca72K3H\nta Cmh5\n%\n3em72K3H\nta Cmh5\n%\nActive enhancers Active enhancers\n75\n30\n50 20\n25 10\n0 0\nPoised enhancers Poised enhancers\n75\n50\n25\n0\n1em4K3H\nta Cmh5\n%\n1em4K3H\nta Cmh5\n%\n100\nPoised enhancers\n75\n30\n50 20\n25 10\n0 0\n1em4K3H\nta Cm5\n%\n1em4K3H\nta Cm5\n%\nActive enhancers 40 Active enhancers\n75\n30\n50 20\n25 10\n0 0\nOther CpGs\n30\n20\n10\nPoised enhancers\n1em4K3H\nta\nCmh5\n%\nPrimed enhancers\n75\n50\n25\n0\n1em4K3H\nta\nCm5\n%\nResponse to reviewers\nThe scatter plot has now been added as suggested (Supplementary Fig. 8, see below) and closely matches the reviewer’s prediction: sites where 6B-C&T > WG follow a trend closer to the y-axis (steep incline), while sites with 6B-C&T < WG follow a trend closer to the x-axis (flat incline), deviating from the diagonal slope.\nThe volcano plots were kept in Extended Data as they provide a more complete view of differentially (hydroxy)methylated sites by integrating the fold-change and statistical significance for each site.\n1.00\n0.75\n0.50\n0.25\n0.00\n0.00 0.25 0.50 0.75 1.00\nSupplementary Figure 8. Differentially methylated CpGs (DMCs) and differentially hydroxymethylated CpGs (DHMCs) between 6B-C&T and whole-genome (untargeted) 6-base-seq. Scatter plots of %5mC (top row) and %5hmC (bottom row) for CpG sites (round dots) captured in both 6B-C&T and whole-genome (WG) 6-base-seq. In each plot, the y-axis coordinate of each CpG site is obtained from 6B-C&T data (replicate-averaged) and the x-axis coordinate is obtained from WG data (replicate-averaged). Differentially methylated sites (p-value < 0.05, |log FC| ³ 0.5) between both methods are coloured for each histone mark\n2\n(H3K27ac in blue, H3K4me3 in orange, H3K4me1 in green and H3K27me3 in red). Non-significant changes are shown in grey.\n12. Extended Data Fig 6: Used to substantiate “that regions annotated as primed enhancers are enriched with the H3K4me1 mark but also show appreciable levels of H3K27ac and H3K27me3.” Profile plots seem to support that point, since the y-axes indicate more H3K4me1 than H3K27me3 or H3K27ac at these primed enhancers. Please also provide genomic heatmaps showing the three histone PTMs (row-linked; sorted descending by H3K4me1 signal) along with an IgG control (as above) to identify background signals.\nWe have now added the suggested heatmaps in Extended Data Fig. 5 (see below).\n17\n)T&C-B6\nca72K3H(\nCm5\n%\n)T&C-B6\nca72K3H(\nCmh5\n%\n)T&C-B6\n3em4K3H(\nCm5\n%\n)T&C-B6\n3em4K3H(\nCmh5\n%\n)T&C-B6\n1em4K3H(\nCmh5\n%\n)T&C-B6\n3em72K3H(\nCmh5\n%\n% 5mC (WG 6B-seq) % 5mC (WG 6B-seq) % 5mC (WG 6B-seq) % 5mC (WG 6B-seq)\n5hmC\n% 5hmC (WG 6B-seq) % 5hmC (WG 6B-seq) % 5hmC (WG 6B-seq) % 5hmC (WG 6B-seq)\n)T&C-B6\n1em4K3H(\nCm5\n%\n)T&C-B6\n3em72K3H(\nCm5\n%\n5mC\n1.00 1.00 1.00\n0.75 0.75 0.75\n0.50 0.50 0.50\n0.25 0.25 0.25\n0.00 0.00 0.00\n0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00\n1.00 1.00 1.00 1.00\n0.75 0.75 0.75 0.75\n0.50 0.50 0.50 0.50\n0.25 0.25 0.25 0.25\n0.00 0.00 0.00 0.00\n0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00\nResponse to reviewers\na b\nActive enhancers Primed enhancers\n-6.0kbcentre6.0kb-6.0kbcentre6.0kb-6.0kbcentre6.0kb-6.0kbcentre6.0kb -6.0kbcentre6.0kb-6.0kbcentre6.0kb-6.0kbcentre6.0kb-6.0kbcentre6.0kb\ngene distance (bp) gene distance (bp)\nc\n-6.0kbcentre6.0kb-6.0kbcentre6.0kb-6.0kbcentre6.0kb-6.0kbcentre6.0kb\ngene distance (bp)\nExtended Data Figure 5. Genomic enrichments of 6B-C&T for enhancer-associated histone modifications at annotated enhancer types. Genomic heatmaps for each histone mark (H3K4me1 in green; H3K27ac in blue; H3K27me3 in red and IgG control in grey) are plotted in +/- 6 kilobase windows relative to the centre of each enhancer locus for a set of annotated (a) active, (b) primed and (c) poised enhancers.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nReviewer #4 (Remarks to the Author):\nThis manuscript by Tavares et al. is a study of the 5mC and 5hmC status of CpGs at\n18\nseneg\nseneg\nseneg\nH3K4me1 H3K27ac H3K27me3 IgG H3K4me1 H3K27ac H3K27me3 IgG\nPoised enhancers\nH3K4me1 H3K27me3 H3K27ac IgG\nResponse to reviewers\npulled down histone modification sites in the genome. The authors use their previously developed 6-base bisulfite free sequencing approach with CUT&Tag. The main advance is the use of the new 6-base bisulfite free approach, over previous versions of this experiment that have been conducted using bisulfite sequencing. While this is an interesting advance, I do have some concerns that will need to be addressed.\n1. There is a lack of data comparison to the previous bisulfite methods, did they see similar things (or not, which is equally interesting).\nSee response to reviewer 2, points 1 and 3.\n2. On page 4, the authors state “6B-C&T-derived 5mC and 5hmC levels at each histone mark were lower than the whole-genome average (Extended Data Fig. 1b)”, which is really interesting.\nHowever, they later say “active chromatin marks (H3K4me3, H3K27ac and H3K4me1) showed similar numbers of CpGs with higher 5hmC levels than in the whole-genome data and CpGs with lower 5hmC levels than in whole genome” and “the repressive mark (H3K27me3) showed the majority of CpG sites with higher 5hmC in 6B-C&T than in whole-genome”\nThese two sets of comments seem to be contradictory. The first that all 5mC/5hmC levels were below WG. Then the second two comments that some 5hmC level were similar to WG and some higher than WG. It didn’t help that the second two sentences didn’t have a figure attached to them.\nThe first comment compares the different genomic regions captured by 6B-C&T vs untargeted (whole-genome) 6-base-seq experiments. 6B-C&T captures a subset of the whole genome, i.e., histone mark-enriched sites. This subset is generally less (hydroxy)methylated than the rest of the genome, which means average 5mC and 5hmC levels obtained from 6B-C&T are lower than in the whole-genome case.\nThe second comment refers to a focused analysis comparing 6B-C&T data to untargeted (whole-genome) data at the same (shared) CpG sites. The observed differences reflect the different DNA fragments captured by the two methods: 6B-C&T fragments are physically associated with each histone mark whereas whole-genome fragments encompass all DNA fragments mapping to the same sites, which includes those not associated with each respective histone mark in the population. This results in the CpG-level differences shown in Extended Data Figure 3.\nWe have rewritten the relevant Results section (pages 3-4, Results paragraphs 5-6) and now mention each figure next to the relevant text. To remove any ambiguity between the two sets of comments, the CpG-level comparison between 6B-C&T and WG is now phrased in terms of differentially methylated CpGs (DMCs) and differentially hydroxymethylated CpGs (DHMCs). New Supplementary Figure 8 contains an alternative visualisation of DMCs and DHMCs in x-y scatter plots (see response to reviewer 2, point 11).\n“Methylation is known to be depleted in histone mark-enriched regulatory regions relative to the other regions of the in mESC genome. Consistent with this, we observed that average 5mC and 5hmC levels in 6B-C&T were lower than in the whole-genome (untargeted) 6-base-seq experiment (Extended Data Fig. 1a). This difference is due to the 6B-C&T method specifically capturing DNA fragments at histone mark-enriched sites, which are generally less (hydroxy)methylated compared to the whole genome.”\n“We next investigated how 5mC and 5hmC levels in 6B-C&T differed from untargeted 6-base-seq data for the same CpG sites (Extended Data Fig. 3a). First, we identified CpG sites with statistically significant differences between both methods, i.e., differentially methylated CpGs (DMCs) and differentially hydroxymethylated CpGs (DHMCs). While most DHMCs between both methods were located in introns followed by intergenic and promoter regions, DMCs were primarily found in promoters followed by introns/exons (Extended Data Fig. 3b-e). Overall, DMCs predominantly showed lower 5mC levels in 6B-C&T relative to whole-genome (untargeted) data (Extended Data Fig. 3a, Supplementary Fig. 8). This suggests that DNA associated with histone modifications at these sites is mostly hypomethylated relative to the total pool of DNA fragments. For DHMCs, active chromatin marks (H3K4me3, H3K27ac and H3K4me1) showed similar numbers of sites with higher or lower 5hmC levels in 6B-C&T relative to the untargeted experiment; in contrast, for the repressive mark (H3K27me3) the majority of DHMCs showed higher 5hmC levels in 6B-C&T than in untargeted 6-base-seq (Extended Data Fig. 3a, Supplementary Fig. 8).”\n3. A general trend throughout is that the 5hmC levels always seem to be a very similar % of the 5mC levels (~25-35%). I’m surprised this isn’t noted and explained. Is this a factor of how the experiments were carried out? Is this identifying something important about 5hmC formation/removal? Is this the same in the WG samples?\nEven though the 5hmC/5mC ratio varies over a relatively narrow range as pointed out by the reviewer, it reveals an interesting histone mark-dependent trend that is not captured by WG data. The following plots (see below), now included as Supplementary Figure 7, compare the relative abundances of 5hmC and 5mC from WG and 6B-C&T data at different histone mark sites.\n100%\n75%\n50%\n25%\n0%\nWG H3K27me3 WG H3K4me1\nWG H3K27ac WG H3K4me3\nSupplementary Figure 7. Relative abundances of 5mC and 5hmC from WG 6-base-seq and 6B-C&T at four histone modifications. In each graph, stacked bar plots show WG and 6B-C&T at CpG sites commonly detected by both experiments. 5mC is represented in dark and 5hmC in light colours. Percentages of each base are annotated for each condition.\n5hmC/5mC from 6B-C&T increases in the following order: H3K27me3 (repressed chromatin, 20%/80%) à H3K4me1 (enhancers, 24%/76%) à H3K27ac (active chromatin, 30%/70%) ~ H3K4me3 (promoters, 31%/69%). This suggests that 5mC removal and 5hmC accumulation might be physically coupled with each histone mark in a transcriptional activation-dependent manner. Conversely, 5hmC/5mC from WG is generally invariant (15%/85%), indicating that these individual trends are masked by the average landscape of DNA fragments at the same sites.\nIn the manuscript, we highlight this histone mark dependence of the 5hmC/5mC ratio in the following section (page 3, Results paragraph 5):\n“Furthermore, when looking at individual histone marks, we found that those associated with active chromatin (H3K4me3, H3K27ac, H3K4me1) showed lower overall methylation and a higher 5hmC/5mC ratio than the bivalent/repressed chromatin (H3K27me3) mark, a pattern that was absent in whole-genome data (Extended Data Fig. 1b, Supplementary Fig. 7). ”\n4. I don’t understand the point of the machine learning model that was developed, perhaps I misunderstand it. It seems they put the WG 5mC/5hmC and H3K4me1 5mC/5hmC (enhancer marker) data into the ML system and then show that the 21\nnoitalyhteM\nGpC\n%\n100%\n75%\n50%\n25%\n0%\nnoitalyhteM\nGpC\n%\n100%\n75%\n50%\n25%\n0%\nnoitalyhteM\nGpC\n%\n100%\n75%\n50%\n25%\n0%\nnoitalyhteM\nGpC\n%\n5mC 5mC 5mC 5mC\n5hmC 5hmC 5hmC 5hmC\n14.6% 20.3% 16.2% 23.8%\n85.4% 79.7% 83.8% 76.2%\n5mC 5mC\n5mC 5mC\n5hmC 5hmC\n5hmC 5hmC\n15.1% 15.5%\n29.5% 30.9%\n84.9% 84.5%\n70.5% 69.1%\nResponse to reviewers\nH3K4me1 version can predict enhancer regions better. But surely that is obvious, as you have preenriched the enhancer regions already?\nTo clarify, the purpose of the ML analysis was to test whether 5mC and 5hmC information for single histone mark (H3K4me1, common to all enhancer types) would robustly predict different enhancer states. It is important to mention that although all enhancer regions are marked by H3K4me1, whether 5mC and 5hmC physically co-occur with this histone modification on the same DNA fragment had never been measured. The results show that these two layers of epigenetic information (i.e., histone and DNA modifications) are coupled and intrinsically linked to different enhancer functional states. Furthermore, our approach demonstrates how informative methylation signatures can be extracted to define different enhancer states without the need for whole-genome sequencing or multiple histone mark enrichment profiles.\nWe have clarified the relevant text:\nResults (pages 4-5, Results paragraph 9):\n“H3K4me1 is common to all enhancer functional states, so we next asked whether 6B-C&T methylation at this single feature better resolves different enhancer types compared to the more heterogenous whole genome (untargeted) 6-base-seq data. For this, we used H3K4me1 6B-C&T data to train a machine learning model for classification of enhancer types (active, primed and poised) and compared its performance to a model with the same architecture trained on whole-genome 6B-seq data (Fig. 4a). For all enhancer types, the 6B-C&T-derived model showed superior performance when compared to the whole genome-derived model (Fig. 4b-d; Extended Data Fig. 7). This demonstrates that the H3K4me1-coupled 5mC and 5hmC signal more robustly resolves bona fide enhancer states than its unenriched equivalent. These results show that 6B-C&T can be used to evaluate co-occurring epigenetic features and retrieve improved information on individual molecular states (Fig. 1b).”\nDiscussion (page 5, Discussion paragraph 2):\n“…although all enhancer functional states are marked by H3K4me1, the co-occurrence of 5mC and 5hmC on the same DNA fragment with this histone modification had never been measured. Here, we show that both DNA marks co-occur at significantly higher levels with H3K4me1 than with other histone marks. Our machine learning analysis shows how these different layers of epigenetic information (i.e., histone and DNA modifications) are coupled and intrinsically linked with enhancer functional states. Importantly, the use of H3K4me1-specific methylation signatures enables the identification of functional enhancer states without the need for whole-genome sequencing or multiple histone mark enrichment profiles.”\n5. Most of the data is not in the main text, I think it would be useful to have a table(s) to summarise the findings as there is a tremendous amount of data that is skimmed over in the text.\nSupplementary Table 8 now contains a summary of all the findings and conclusions in Extended Data.\n6. In plots Fig 2b, Ext data Fig 4, and Ext data Fig 6b, the violin plots show a very interesting trend. It’s almost as if the %5mC/5hmC is quantised. There are specific a bumps of data at regular intervals. The authors should explain why this is, as I’m assuming this cannot be biological so will be a factor of how their experiments were conducts or data processed? What effect does this have on the conclusions?\nThe nature of 6-base-seq data is discrete since the fraction of (hydroxy)methylation at CpG sites is a function of the (hydroxy)methylation counts and the coverage at each site, both of which are discrete variables (fraction = counts/coverage).\nAt fixed coverage:\n1 methylation count covered 20 times = 5%\n2 counts covered 20 times = 10%\nb\nAt fixed counts:\n1 methylation count at a site covered 20 times = 5%\n1 methylation count at a site covered 40 times = 2.5%\nGiven the coverage threshold routinely employed for methylation calling (e.g., 10x, 20x) restricts the coverage parameter (denominator) to a very narrow range, the data values (%5hmC or %5mC) become a more direct function of discrete methylation counts (numerator). This produces a “banding” pattern in graphs showing values at individual CpG sites such as scatter plots (see example below) and violin plots. This pattern is absent in plots generated after data smoothing (e.g., window averaging) or descriptive statistics representations (e.g., boxplots).\nc\nd\nExtracted from Extended Data Fig. 2c – Illustration of methylation “bands” in 6-base-seq data.\nIn violin plots generated with ggplot2 in R, the violin shape is constructed by applying a kernel density estimate (KDE) to the data distribution. This smoothing process can introduce artifacts – such as spurious peaks or apparent gaps – around certain data\n23\n08\n06\n04\n02\nH3K27ac enrichment H3K27ac enrichment H3K27ac enrichment H3K27ac enrichment\nH3K4me1 enrichment H3K4me1 enrichment H3K4me1 enrichment H3K4me1 enrichment\nH3K4me3 enrichment H3K4me3 enrichment H3K4me3 enrichment H3K4me3 enrichment\nH3K27me3 enrichment H3K27me3 enrichment H3K27me3 enrichment H3K27me3 enrichment\n)GW(\n%Cm5\n)GW(\n%Cm5\n)GW(\n%Cm5\n)GW(\n%Cm5\n)T&C-L6\nca72K3H(\n%Cm5\n)T&C-L6\n1em4K3H(\n%Cm5\n)T&C-L6\n3em4K3H(\n%Cm5\n)T&C-L6\n3em72K3H(\n%Cm5\n)GW(\n%Cmh5\n)GW(\n%Cmh5\n)GW(\n%Cmh5\n)GW(\n%Cmh5\n)T&C-L6\nca72K3H(\n%Cmh5\n)T&C-L6\n1em4K3H(\n%Cmh5\n)T&C-L6\n3em4K3H(\n%Cmh5\n)T&C-L6\n3em72K3H(\n%Cmh5\n6L-CUT&Tag, shared CpGs WG 6-letter-seq, shared CpGs 6L-CUT&Tag, shared CpGs\nn=42313 n=42313 n=42313 n=42313\n−2 0 2 4 6 8 10 12 −2 0 2 4 6 8 10 12\nWG 6-letter-seq, shared CpGs 6L-CUT&Tag, shared CpGs WG 6-letter-seq, shared CpGs 6L-CUT&Tag, shared CpGs\nn=75292 n=75292 n=75292 n=75292\nWG 6-letter-seq, shared CpGs 6L-CUT&Tag, shared CpGs WG 6-letter-seq, shared CpGs 6L-CUT&Tag, shared CpGs\nn=412646 n=412646 n=412646 n=412646\nWG 6-letter-seq, shared CpGs 6L-CUT&Tag, shared CpGs WG 6-letter-seq, shared CpGs 6L-CUT&Tag, shared CpGs\nn=27381 n=27381 n=27381 n=27381\n001\n08\n06\n04\n02\n0\n−2 0 2 4 6 8 10 12\n53\n03\n52\n02\n51\n01\n5\n0\n−2 0 2 4 6 8 10 12\n04\n03\n02\n01\n0\n−2 0 2 4 6 8 10 12\n001\n08\n06\n04\n02\n0\n−2 0 2 4 6 8 10 12\n001\n08\n06\n04\n02\n0\n−2 0 2 4 6 8 10 12\n53\n03\n52\n02\n51\n01\n5\n0\n−2 0 2 4 6 8 10 12\n04\n03\n02\n01\n0\n0 2 4 6 8 10\n08\n06\n04\n02\n0\n0 2 4 6 8 10\n08\n06\n04\n02\n0\n0 2 4 6 8 10\n04\n03\n02\n01\n0\n0 2 4 6 8 10\n05\n04\n03\n02\n01\n0\n0 2 4 6 8 10\n08\n06\n04\n02\n0\n0 2 4 6 8 10\n001\n08\n06\n04\n02\n0\n0 2 4 6 8 10\n03\n52\n02\n51\n01\n5\n0\n0 2 4 6 8 10\n03\n02\n01\nvalues, particularly in discrete datasets. As a result, the density envelope may visually exaggerate or underrepresent the presence of data in specific intervals, potentially giving the impression that some intervals contain few or no observations when data points are in fact present. For this reason, we now use stacked bar plots to describe the data in Main Fig. 3b and moved the violin plots to Extended Data Fig. 6 to provide a means of alternative visualisation. Importantly, the above technical observations and visualisation changes do not affect the conclusions in the manuscript.\nReviewer #1 (Remarks to the Author):\nAuthors addressed concerns raised in the initial submission of the manuscript, with following remaining points:\n1. In imprinted gene loci authors demonstrated that e.g. 6B-C&T for H3K4me3 enriches for low 5mC loci, despite close to 50% 5mC present in WG dataset. These are predicted to be the most dramatic examples in the genome, where half of the alleles exist in two completely distinct chromatin and DNA methylation states. While authors decided not to put this figure in the manuscript, I suggest that they should highlight imprinted genes in volcano plots (for at least H3K4me3).\nWe appreciate the suggestion and have annotated CpG sites from the imprinted gene loci Igf2r, Kcnq1ot and Peg3 in the volcano plots in Supplementary Figure 10a and corresponding text:\nSupplementary Figure 10a. Differentially methylated and hydroxymethylated CpGs between 6B-C&T and whole-genome 6-base-seq. (a) Volcano plots of differentially methylated 5mCpGs (DMCs, top) and 5hmCpGs (DHMCs, bottom) for each 6B-C&T dataset relative to whole-genome 6-base-seq (WG). Differentially methylated sites (p-value < 0.05, |log FC| ³ 0.5) are coloured for each histone mark.\nThe percentage of significant differences where 6B-C&T values are lower (left) or higher (right) than in whole-genome data are indicated in each plot. CpG sites at example promoter regions showing allele-specific expression in E14TG2A mESCs (imprinted loci) are annotated in each plot (Igf2r = white filled circles; Kcnq1ot = white filled squares; Peg3 = white filled triangles).\nText:\n“Overall, DMCs predominantly showed lower 5mC levels in 6B-C&T relative to whole-genome (untargeted) data (Supplementary Figure 10a, Supplementary Figure 11). This suggests that DNA associated with histone modifications at these sites is mostly hypomethylated relative to the total pool of DNA fragments. This was\n)eulav−p(\n01gol− 5mC 5mC 5mC 5mC\nlog2FC log2FC\n5hmC 5hmC 5hmC 5hmC\nlog2FC\n)eulav−p(\n01gol−\nlog2FC\n)eulav−p(\n01gol−\nlog2FC\n)eulav−p(\n01gol−\n)eulav−p(\n01gol−\n)eulav−p(\n01gol−\nlog2FC\n)eulav−p(\n01gol−\n)eulav−p(\n01gol−\nH3K27me3 H3K4me3 H3K27ac H3K4me1\nvs vs vs vs\nWG WG WG WG\n50\n62% 38% 40 92% 8% 97% 3% 30 92% 8%\n20 20\n30\n20\n10 20 10\n10\n10\n0 0 0 0\n-50 -25 0 25 50 -50 0 -60 -30 0 30 -75 -50 -25 0 25\nlog2FC\n4 1199%% 8811%% 8 5544%% 4466%% 4 5500%% 5500%% 5 4477%% 5533%% 3 4\n6 3 2 3\n4 2 2\n1\n2 1 1\n0 0 0 0\n-20 0 20 40 -20 0 20 40 -30-20-10 0 10 20 30 -20 0 20\nlog2FC\nparticularly striking in DMCs co-occurring with H3K4me3 at imprinted gene loci39 such as Igfr2, Kcnq1ot1 and Peg3, suggesting that this active histone mark largely associates with the unmethylated allele (Supplementary Figure 10a).”\n2. Additional introduction is welcome, but it has a factual error. Statement “5mC is the most abundant cytosine modification (~60% of all Cs in ESCs)9 and 5hmC is the next most abundant cytosine modification (less than 5% of all Cs in ESCs)9 …” is incorrect. 5mC constitutes close to 2 - 4% of all Cs. Authors likely meant that 60% of Cs in CpG dinucleotides are methylated, which is closer the actual methylation value. The paper that was referenced does measure modification abundance in CpGs since they perform TLC on labelled nucleotides generated from DNA digestion with MspI, HpaII or Taqα1 – all of them are employed for labelling of Cs in CpGs. Tahiliani et al incorrectly labelled axis in some of their figures, which translated to confusion here.\nWe correct this oversight (see below). To avoid confusion, we now cite the study by Ito et al. (PMID: 21778364), where the authors quantified the genomic content of all these cytosine modifications using mass spectrometry.\n“5mC is the most abundant cytosine modification (2-4% of all Cs in ESCs) and 5hmC is the next most abundant cytosine modification (0.1-0.2% of all Cs in ESCs), with 5fC and 5caC being present at substantially lower levels (~0.002% and ~0.0003%, respectively)8.”\n3. Citations - while listing methods they should also list TAPS (PMID: 30804537) and CAPS (PMID: 33504799) as they are conceptually different from methods that relay on deamination of unmodified Cs (like bisulfite or EM-seq).\nWe have added the above pyridine borane sequencing methods to the Introduction:\n“While early methods such as bisulfite sequencing19 cannot distinguish 5hmC from 5mC, later approaches – including oxidative bisulfite sequencing20, TET-assisted bisulfite sequencing21, APOBEC-coupled epigenetic sequencing22, pyridine borane sequencing methods23, 24 and third-generation sequencing platforms (e.g., PacBio, Oxford Nanopore)25 – have enabled base-resolution mapping of both 5mC and 5hmC.”\nReviewer #2 (Remarks to the Author):\nGeneral Comments [Co-reviewers #2/3]\nThe authors have thoroughly addressed our original critiques by an improved introduction (with additional background and technology context), addition of essential IgG control reactions, and additional analyses / figures (including reduced claims of enhancer biology in areas of low PTM signal). The rebuttal is a thorough, good-faith effort to address concerns and describe the manuscript revisions. As such we support publication.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nReviewer #4 (Remarks to the Author):\nI am happy my comments have been addressed and supportive of publication"}]

Decision letter