Mouse embryonic stem cells were provided by Matthias Stadtfeld. Briefly, KH2 embryonic stem cells (ESCs)  were cultured on irradiated feeder cells in KO-DMEM (Invitrogen) supplemented with l-glutamine, penicillin/streptomycin, nonessential amino acids, β-mercaptoethanol, 1000 U/mL LIF, and 15% FBS (ESC medium). To remove feeder cells from ESCs, cells were trypsin digested and pre-plated in ESC medium for 30 min. Supernatant containing ESCs was used for further experiments.
Mouse KLF4 has been cloned into pHAGE2-tetO-MCS-ires-tdTomato vector (obtained from Matthias Stadfeld’s lab, ) for the production of lentiviruses, using the following primers:
Fwd: 5′– gcggccgcATGGCTGTCAGCGACGCTCT
Rev: 5′– ggatccTTAAAAGTGCCTCTTCATGTGTAAGG
KLF4 R462A mutation has been generated using the site-directed mutagenesis kit from Agilent #210518. HEK 293T cells were used for the production of lentiviruses, obtained from ATCC (cat. No. CRL 3216). Lentiviral infection of KLF4 knockout mESC  was performed by spin-infection, and the cells were transferred to feeders and expanded with puromycin. After selection, KLF4 expression was induced with doxycycline (1μg/ml) for 2 days. Finally, the cells were pre-seeded (30 min) to remove the feeders, and the ES cells were processed as described in the “Cell culture” section. KLF4 protein expression has been checked by western blot using an antibody from Santa Cruz (#sc-20691, now discontinued) and using H3 as a loading control (anti-H3, Abcam, ab1791).
Assembly of the transposase
Tn5 transposase was assembled with methylated adaptors as per the T-WGBS protocol . Ten microliters of each adapter with incorporated methylated cytosines (Tn5mC-Apt1 and Tn5mC1.1-A1block; 100 μM each; Additional file 2: Table S1) was added to 80 μl of water and annealed in a thermomixer with the following program: 95 °C for 3 min, 70 °C for 3 min, 45 cycles of 30 s with a ramp at − 1 °C per cycle to reach 26 °C. Fifty microliters of annealed adapters was incubated with 50 μl of hot glycerol and 10 μl of this mixture was incubated with 10 μl of Ez-Tn5 transposase (from the EZ-Tn5 insertion kit) at room temperature for 30 min to assemble the transposome.
ATAC-seq and M-ATAC
ATAC-seq and M-ATAC were performed with 50,000 mESC as per the original ATAC-seq protocol . Cells were washed in cold PBS and resuspended in 50 μl of cold lysis buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% IGEPAL CA-630). The tagmentation reaction was performed in 25 μl of TD buffer (Illumina Cat #FC-121-1030), 2.5 μl transposase (either the Nextera transposase (ATAC-seq), or the transposase containing the methylated adaptors (M-ATAC, see section “Assembly of the transposase” for details), and 22.5 μl of nuclease-free H2O at 37 °C for 30 min. Purified DNA (on column with the Qiagen Mini Elute kit) either bisulfite converted (M-ATAC, see section “Bisulfite conversion” for details) or directly amplified (ATAC-seq, see “Amplification of ATAC-seq and ChIP-seq libraries” for details).
ChIP-seq and M-ChIP
ChIP-seq and M-ChIP were performed on mESC as per the original ChIPmentation protocol . Five microliters of CTCF antibody (Millipore 07-729) or 25 μl of KLF4 antibody (R&D AF3158) was combined to protein A (for CTCF) or G (for KLF4) magnetic beads and added to sonicated chromatin (from 200 to 700 bp, checked on agarose gel) from 10 million mESC, for 3 to 6 h rotating in the cold room. Beads were washed as per the original ChIPmentation protocol : twice with TF-WBI (20 mM Tris-HCl/pH 7.4, 150 mM NaCl, 0.1% SDS, 1% Triton X − 100, 2 mM EDTA), twice with TF-WBIII (250 mM LiCl, 1% Triton X-100, 0.7% DOC, and 10 mM Tris-HCl, 1 mM EDTA), and twice with cold Tris-Cl pH 8.0 to remove detergent, salts, and EDTA. During the second wash, the whole reaction was transferred to a new tube to decrease tagmentation of unspecific chromatin fragments sticking to the tube wall. Beads were resuspended in 25 μl of the tagmentation reaction mix (10 mM Tris pH 8.0, 5 mM MgCl2, and 10% v/v dimethylformamide), and tagmentation was performed for 1 min at 37 °C with either 1 μl of the Nextera transposase (ChIP-seq) or the transposase containing the methylated adaptors (M-ChIP, see section “Assembly of the transposase” for details). Then, beads were washed twice with TF-WBI (20 mM Tris-HCl/pH 7.4, 150 mM NaCl, 0.1% SDS, 1% Triton X − 100, and 2 mM EDTA) and twice with TET (0.2% Tween − 20, 10 mM Tris-HCl/pH 8.0, 1 mM EDTA). During the last wash, the whole reaction was transferred to a new tube to decrease carry-over of tagmented unspecific fragments stuck to the tube wall. Chromatin was eluted and decrosslinked by 70 μl of elution buffer (0.5% SDS, 300 mM NaCl, 5 mM EDTA, 10 mM Tris-HCl pH 8.0) containing 20 μg of proteinase K for 2 h at 55 °C and overnight incubation at 65 °C. Eluted and purified DNA was either bisulfite converted (CTCF M-ChIP, see section “Bisulfite conversion” for details) or directly amplified (CTCF ChIP-seq, see “Amplification of ATAC-seq and ChIP-seq libraries” for details).
Purified DNA was bisulfite converted following the T-WGBS protocol  with the EZ DNA methylation kit (Zymo). Oligonucleotide replacement was performed by incubating 9 μl of tagmented M-ATAC or M-ChIP purified DNA with 2 ng of phage lambda DNA as carrier, 2 μl of dNTP mix (2.5 mM each, 10 mM), 2 μl of 10× Ampligase buffer, and 2 μl of replacement oligo (Tn5mC-ReplO1, 10 μM; Additional file 2: Table S1) in a thermomixer with the following program: 50 °C for 1 min, 45 °C for 10 min, ramp at − 0.1 °C per second to reach 37 °C. One microliter of T4 DNA polymerase and 2.5 μl of Ampligase were added, and the gap repair reaction was performed at 37 °C for 30 min. DNA was purified using SPRI AMPure XP beads with a bead-to-sample ratio of 1.8:1 and eluted in 50 μl of H2O. Five microliters was kept as an unconverted control sample, and 45 μl was bisulfite converted using the EZ DNA methylation kit (Zymo). Briefly, the gap repair reaction was performed by adding 5 μl of M-dilution buffer and 15 min incubation at 37 °C, and bisulfite treatment was performed by adding 100 μl of liquid CT-conversion reagent in a thermomixer with the following program: 16 cycles of 95 °C for 15 s followed by 50 °C for 1 h. Converted DNA was purified on a column and amplified (see section “Amplification of M-ATAC and M-ChIP libraries” for details).
Amplification of ATAC-seq and ChIP-seq libraries
Purified DNA (20 μl) was combined with 2.5 μl of each primer and 25 μl of NEB Next PCR master mix as per the original ATAC-seq protocol . For ATAC-seq, DNA was amplified for 5 cycles and a monitored quantitative PCR was performed to determine the number of extra cycles needed not exceeding 12 cycles in total to limit the percentage of duplicated reads. DNA was purified on a column with the Qiagen Mini Elute kit. For ChIP-seq, DNA was amplified as per the ChIPmentation protocol  in a thermomixer with the following program: 72 °C for 5 min; 98 °C for 30 s; 14 cycles of 98 °C for 10 s, 63 °C for 30 s and 72 °C 30 s; and a final elongation at 72 °C for 1 min. DNA was purified using SPRI AMPure XP beads with a bead-to-sample ratio of 1:1 and eluted in 20 μl of H2O.
Amplification of M-ATAC and M-ChIP libraries
Purified converted DNA was amplified as per the original T-WGBS protocol . Briefly, 10 μl of DNA was combined with 1.25 μl of each primer (25 μM each) and 12.5 μl of high-fidelity system KAPA HiFi uracil+ PCR master mix. DNA was amplified for 5 cycles, and a monitored quantitative PCR was performed to determine the number of extra cycles needed, not exceeding 12 cycles in total to limit the percentage of duplicated reads.
Sequencing of the libraries and data processing
For ATAC-seq, ChIP-seq, M-ATAC, and M-ChIP, libraries were quantified using Kapa qPCR kit and sequenced using the HiSeq 2500 for paired-end 50-bp reads. ChIP-seq for histone modifications in mESC were downloaded from GEO (H3K4me1: GSM1000121, H3K27ac: GSM1000126, H3K4me3: GSM1000124). Data processing was performed as per the pipeline available on Github (https://github.com/skoklab/EpiMethylTag). Briefly, reads were trimmed using trim-galore/0.4.4, and aligned to the mm10 assembly of mouse genome using bowtie2  for ChIP-seq and ATAC-seq, and using Bismark/0.18.1 (bowtie2)  for M-ChIP and M-ATAC to account for bisulfite conversion. Reads with quality < 30 and duplicates were removed using Samtools/1.3 . Peaks were called using Macs/2.1.0  with the following parameters: --qvalue 0.01 --nomodel --shift 0 -B --call-summits. Narrow peaks were considered for further analysis. Bigwigs were generated from bam files with RPKM normalization using Deeptools  for visualization on IGV.
Bioinformatic analysis of data
The distribution of fragment lengths was assessed with Deeptools/2.3.3 with option “--maxFragmentLength 1000”, and Pearson correlations of read counts with Deeptools/2.3.3 and default parameters. Heatmaps and average profiles were performed on merged bigwig files using Deeptools/2.3.3. Default parameters from Bismark/0.18.1 (Bowtie2)  were used to generate coverage files containing methylation information. Only cytosines in a CpG context were used for subsequent analysis. For Fig. 3d and Additional file 1: Figure S5d, e, the plots were centered on CpGs in M-ATAC peaks from the different groups highlighted in Fig. 3a. For Fig. 4a, lists of CpGs were subsampled using BEDTools  to consider only the CpGs inside CTCF motifs, and the average plots were centered on those CpGs. Genomic annotations were performed using ChIPseeker . CTCF motif locations in CTCF M-ChIP/ChIP and M-ATAC, and KLF4 motifs in M-ChIP peaks were determined using the FIMO tool from MEME , with the motif PWM from Jaspar database (MA0139.1 for CTCF and MA0039.1 and MA0039.2 for KLF4). PWM was manually modified to look at methylation frequency at different combinations of C2 and C12 dinucleotides of CTCF motif. Scripts are available on Github (https://github.com/skoklab/EpiMethylTag). In order to account for possible lack of specificity of the anti-KLF4 antibody, we filtered out ChIP-seq peaks present in Klf4−/− cells. Peaks shared or specific to either WT or mutant KLF4 were identified using BEDTools . For the ChIP enrichment versus CpG methylation plots, we plotted the peak score versus the beta values of the CpG probes within the peaks, using peaks called via MACS2 for CTCF (Additional file 1: Figure S6b) and via PeaKDEck for KLF4 (Fig. 4c).
To quantify the probability of clustering CpG probes into low, medium, and highly methylated groups, we assumed that beta values (i.e., the sampling mean) are normally distributed with the mean beta value (b) and variance (b (1 − b))/((n − 1)) where n is the total number of reads. This allows us to quantify the probability that each probe belongs to its designated cluster as P(b < Ch) − P(b < Cl) where Ch and Cl are the high and low thresholds of the cluster respectively. In Additional file 1: Figure S6g, the points and corresponding contours are colored based on their designated cluster. The x-axis is the beta value and the y-axis is the probability that beta lies within the cluster limits. For all clusters, more than half of the CpGs have a high probability of being in the assigned group (> 72%).