Epigenetic reprogramming by TET enzymes impacts co-transcriptional R-loops

DNA oxidation by ten-eleven translocation (TET) family enzymes is essential for epigenetic reprogramming. The conversion of 5-methylcytosine (5mC) into 5-hydroxymethylcytosine (5hmC) initiates developmental and cell-type-specific transcriptional programs through mechanisms that include changes in the chromatin structure. Here, we show that the presence of 5hmC in the transcribed gene promotes the annealing of the nascent RNA to the template DNA strand, leading to the formation of an R-loop. Depletion of TET enzymes reduced global R-loops in the absence of gene expression changes, whereas CRISPR-mediated tethering of TET to an active gene promoted the formation of R-loops. The genome-wide distribution of 5hmC and R-loops shows a positive correlation in mouse and human stem cells and overlap in half of all active genes. Moreover, R-loop resolution leads to differential expression of a subset of genes that are involved in crucial events during stem cell proliferation. Altogether, our data reveal that epigenetic reprogramming via TET activity promotes co-transcriptional R-loop formation, disclosing new mechanisms of gene expression regulation.


Introduction
During transcription, the nascent RNA molecule can hybridize with the template DNA and form a DNA:RNA hybrid and a displaced DNA strand. These triple-stranded structures, called R-loops, are physiologically relevant intermediates of several processes, such as immunoglobulin class-switch recombination and gene expression (García-Muse and Aguilera, 2019). However, nonscheduled or persistent R-loops constitute an important source of DNA damage, namely, DNA double-strand targeting the enzyme to an active gene drives R-loop formation. Analysis of genome-wide distribution profiles shows a positive correlation between 5hmC and R-loops in mouse embryonic stem (mES) and in human embryonic kidney 293 (HEK293) cells, with a clear overlap of 5hmC and R-loops in approximately half of all active genes. We also show that 5hmC-rich regions are characterized by increased levels of phosphorylated histone H2AX (γH2AX), a marker of DNA damage. Finally, by determining the pathways more significantly affected by R-loops formed at 5hmC loci, we disclose novel links between R-loops and gene expression programs of stem cells.

Results
Transcription through 5hmC-rich DNA favors R-loop formation To assess the impact of cytosine methylation on R-loop formation, we performed in vitro T7 transcription of DNA fragments containing either native or modified cytosine deoxyribonucleotides (dCTPs). We synthesized three distinct DNA transcription templates, each composed of a T7 promoter followed by a 400 bp sequence containing a genomic region prone to form R-loops in vivo . Two of these sequences (ACTB P1 and ACTB P2) are from the transcription termination region of the human β-actin coding gene (ACTB); the third sequence is from the human APOE gene. The DNA templates for the in vitro transcription reactions were generated by PCR-amplification in the presence of dNTPs containing either native C, 5mC, or 5hmC ( Figure 1A). Successful incorporation of dCTP variants was confirmed by immunoblotting using specific antibodies against each variant ( Figure 1B). The formation of R-loops during the in vitro transcription reactions was inspected by blotting immobilized RNAs with the S9.6 antibody (S9.6 Ab), which binds DNA:RNA hybrids ( Figure 1C). To increase the specificity of hybrid detection, all samples were treated with RNase A in high-salt conditions in order to digest all RNA molecules except those engaged in R-loops. The specific detection of DNA:RNA hybrids was confirmed by blotting transcription reaction products previously digested with RNase H ( Figure 1C). In agreement with our hypothesis that 5hmC favors R-loops, increased amounts of DNA:RNA hybrids were detected in samples derived from in vitro transcription of 5hmC-rich ACTB P1, ACTB P2, and APOE DNA templates ( Figure 1D). To exclude the possibility that our results were biased by an inherent preference of the S9.6 Ab for hybrids containing 5hmC, we performed electrophoretic mobility shift assays (EMSAs) using the S9.6 Ab and DNA:RNA hybrid substrates of the same sequence but containing C, 5mC, or 5hmC. The S9.6 Ab was able to delay the run of the three substrates with similar kinetics, indicating that the Ab equally recognizes DNA:RNA hybrids formed with any of the three C variants (Figure 1-figure supplement 1).
We then performed atomic force microscopy (AFM) to directly visualize R-loop structures obtained in the in vitro transcription reactions ( Figure 1E). R-loops were identified as previously described (Carrasco-Salas et al., 2019;Klinov et al., 1998). Each individual DNA molecule establishing an R-loop structure in the AFM images was assigned manually. The frequency of these structures formed in the presence of C, 5hmC, or 5mC DNA templates was measured and normalized against the frequency formed in RNase H-treated samples ( Figure 1E). In agreement with the hypothesis that transcription of 5hmC-rich DNA templates results in increased R-loop formation, AFM data revealed that R-loop structures are more frequently formed in the presence of 5hmC.
To investigate if the DNA modification impacts in vitro transcription levels, we measured RNA synthesis from DNA templates containing unmodified C, 5hmC, or 5mC (Figure 1-figure supplement 2). These data show that the T7 polymerase is highly sensitive to DNA modifications since replacing C by either 5hmC or 5mC significantly decreased the transcript levels in vitro. On one side, detecting more R-loops on a lower-transcription levels setting (i.e., 5hmC-rich templates) further strengthens our hypothesis that 5hmC increases R-loop formation. However, we cannot draw any conclusions regarding the impact of 5mC on R-loop formation as a putative effect on R-loop levels could be masked by the significantly altered transcription. To clarify this aspect and further test our model, we continued our study with experiments performed in vivo.

TET enzymatic activity impacts endogenous R-loop levels
To test whether the 5hmC DNA modification induces R-loop formation in vivo, we quantified R-loop levels in mES cells after depletion of Tet1, Tet2, and Tet3. Despite the significant reduction in Tet enzymes ( Figure 2-figure supplement 1A), the levels of 5hmC were not significantly affected by Figure 1. 5-Hydroxymethylcytosine (5hmC) favors co-transcriptional R-loop formation. (A) Native or modified deoxyribonucleotides (dCTPs) were incorporated upon PCR amplification into DNA fragments with sequences from the transcription termination region of ACTB (ACTB P1 and ACTB P2) or APOE. (B) Incorporation of dCTP variants confirmed by immunoblotting using specific antibodies against 5-methylcytosine (5mC), 5hmC, and double-stranded DNA (dsDNA). (C) R-loops formed upon in vitro transcription reactions were detected by immunoblotting using the S9.6 antibody. RNase H-treated in vitro transcription reaction products (RH+) serve as negative controls. All data are representative of seven independent experiments with similar results. (D) S9.6 immunoblots were quantified and the R-loop levels normalized against the levels detected in the reaction products of DNA templates containing native C. Data represent the mean and standard deviation (SD) from seven independent experiments. *p<0.05, **p<0.01, and ***p<0.001, two-tailed Student's t-test. (E) In vitro transcription reaction products of ACTB P2 templates were visualized using atomic force microscopy (AFM). R-loop structures obtained from 5hmC-containing ACTB P2 transcription in the absence (RH-) or presence (RH+) of RNase H are shown. R-loops present in the transcription reaction products of C, 5mC, or 5hmC-containing ACTB P2 templates were counted in a minimum of 80 filaments observed in three individual AFM experiments.
The online version of this article includes the following figure supplement(s) for figure 1:  Tet1 or Tet2 depletion (Figure 2A, Figure 2-figure supplement 1B). In contrast, depletion of Tet3 resulted in a significant loss of 5hmC, an effect that was exacerbated by the simultaneous depletion of the three enzymes ( Figure 2A, Figure 2-figure supplement 1B). No significant changes were observed in 5mC levels ( Figure 2A, Figure 2-figure supplement 1B). This finding suggests that there is a partial redundancy in the activity of the three Tet enzymes in mES cells. The loss of Tet1 or Tet2 -but not of Tet3 -is compensated by the remaining Tets. In agreement with the hypothesis that 5hmC promotes R-loop formation, dot-blot hybridization of total cellular nucleic acids using the S9.6 Ab revealed reduced endogenous R-loop levels in mES cells after depletion of Tet3 and after co-depletion of the three Tet enzymes ( Figure 2B, Figure 2-figure supplement 1B). We also measured R-loop levels upon RNAi depletion of the Tet enzymes in NIH-3T3 mouse fibroblasts (Figure 2figure supplement 2A). As observed in mES cells, a significant reduction of 5hmC, but not 5mC, was obtained upon depletion of Tet3 and of the three Tets in mouse fibroblasts (Figure 2-figure  supplement 2B, C). The triple knockdown of the Tet enzymes significantly reduced the global levels of R-loops in mouse fibroblasts, whereas Tet3 depletion in these cells had a minor impact in R-loops ( Figure 2-figure supplement 2B and D). This effect was further confirmed by measuring R-loops formed at selected active genes by DNA:RNA immunoprecipitation (DRIP) in mES cells ( Figure 2C) and mouse fibroblasts (Figure 2-figure supplement 2E). The DRIP assays confirmed that R-loops are less abundant upon depletion of Tet enzymes. Importantly, simultaneous depletion of the three enzymes did not affect the expression levels of the analyzed genes in mES cells and mouse fibroblasts ( Figure 2D, Figure 2-figure supplement 2F). These data suggest that the activity of Tet enzymes promotes the formation of R-loops in the absence of changes in transcription levels.
Next, we employed a modified CRISPR-based system to target TET enzymatic activity to specific loci (Liu et al., 2016). We used a pool of three specific guide RNAs (gRNAs) to direct a catalytically inactive Cas9 nuclease fused to the catalytic domain of TET1 (dCas9-TET1) to the last exon of the APOE gene in human osteosarcoma (U-2 OS) cells. As a control, dCas9 was fused to an inactive mutant version of the TET1 catalytic domain (dCas9-dTET1). Local enrichment of 5hmC following dCas9-TET1 targeting at the APOE locus was confirmed by DNA immunoprecipitation using antibodies specific for 5mC or 5hmC-modified nucleotides ( Figure 2E). The highest levels of 5hmC were detected at the gene segment adjacent to the gRNAs-target region. R-loop levels detected by DRIP peaked significantly at the gRNAs-target and in the downstream region, upon tethering of dCas9-TET1 but not of dCas9-dTET1 ( Figure 2F). These differences were not caused by changes in APOE gene expression levels ( Figure 2G). The increased levels of R-loops detected far from the dCas9-TET1 target site are consistent with the view that R-loops have the capacity to extend from their inception locus. Accordingly, R-loops can be up to several hundred base pairs long and may extend over the entire gene body of shorter and/or highly transcribed genes (Sanz et al., 2016;Chen et al., 2017). Collectively, these data suggest that editing 5hmC density by changing the expression levels or the genomic distribution of TET enzymes influences R-loop formation in cells.
5hmC and R-loops overlap genome-wide at transcriptionally active genes To further inspect the link between 5hmC and R-loops, we performed computational analyses of 5hmC antibody-based DNA immunoprecipitation (hMeDIP-seq) and DNA:RNA immunoprecipitation (DRIP-seq) datasets from mES and HEK293 cells Matarese et al., 2011;Jin et al., 2014;Nadel et al., 2015). To assess individual genome-wide distribution profiles, R-loops density was probed over fixed windows of ±10 kbp around the 5hmC peaks ( Figure 3A, Figure 3-figure supplement 1A). The resulting metagene plots and heatmaps revealed a marked overlap between 5hmC-rich loci and R-loops. This overlap is also evident in the individual distribution profiles of 5hmC and R-loops along two long regions of chromosome 17 (Figure 3-figure supplement 2). Despite the distinct distribution patterns of 5hmC (well-defined peaks) and R-loops (reads spanning genomic regions with highly heterogeneous lengths, ranging between a few dozen to over 1 kb Chen et al., 2015), we could obtain a statistically significant Pearson correlation coefficient between both (p<0.05) ( Figure 3B  Metagene profiles revealed very similar patterns of intragenic distribution, with both 5hmC and R-loops increasing towards the transcription termination site (TTS), where they reached maximum levels ( Figure 3F). At the TSS, however, the 5hmC DNA modification was mostly absent, whereas R-loops were abundant. The detection of R-loop peaks at TSS regions is in agreement with previous studies (Ginno et al., 2013;Ginno et al., 2012) and implies that 5hmC is not necessary for co-transcriptional DNA:RNA hybridization and R-loop formation.
The observed overlap between 5hmC and R-loop peaks at the TTS raises the hypothesis that Tet activity may be involved in transcription termination by directing the formation of R-loops. Defects in transcription termination result in the accumulation of readthrough transcripts extending beyond the TTS (Nojima and Proudfoot, 2022). In agreement with a role in transcription termination, TET1-KO human ES cells displayed significantly higher levels of readthrough transcripts genome wide when compared to wt human ES cells ( Figure 3G).
We then sought to simultaneously detect 5hmC and R-loops at the same loci in individual mES cells. We performed proximity ligation assays (PLAs) using S9.6 and anti-5hmC antibodies ( Figure 4A). Control reactions without primary antibodies and with each antibody alone did not produce a significant signal. Staining of mES cells with S9.6 and anti-5hmC antibodies gave rise to a robust PLA signal scattered throughout the nucleus, which was mostly lost after digestion of cells with RNase H ( Figure 4B).

5hmC-rich loci are prone to DNA damage
Disruption of R-loop homeostasis is a well-described source of genomic instability (García-Muse and Aguilera, 2019). For instance, co-transcriptional R-loops increase conflicts between transcription and replication machineries by creating an additional barrier to fork progression (Hamperl et al., 2017;Helmrich et al., 2013). Such conflicts may cause DNA damage, including DSBs, which can be revealed using antibodies against γH2AX. Indeed, R-loops overlap with γH2AX-decorated chromatin at different locations such as TTS (Hatchi et al., 2015). We then sought to investigate if 5hmC creates conditions for DNA damage by promoting R-loop formation. We analyzed the genomic distribution of γH2AX by interrogating chromatin immunoprecipitation followed by sequencing (ChIP-seq) data from HEK293 cells (Bunch et al., 2015). The individual distribution profiles of γH2AX were analyzed over fixed windows of ±10 kbp around the 5hmC peaks detected in the same cells ( Figure 5A). The resulting metagene plots revealed marked enrichment of γH2AX at 5hmC-rich loci. The genic distribution of 5hmC and R-loops along three different genes further showed co-localization of the two marks with γH2AX ( Figure 5B). Analysis of γH2AX and 5hmC distribution within active genes revealed a low yet statistically significant Pearson correlation coefficient (p<0.05) ( Figure 5C).

R-loops formed at 5hmC-rich regions impact gene expression in mES cells
To gather insights into the functional impact of R-loops at 5hmC-rich DNA regions, we analyzed wholetranscriptome (RNA-seq) of mES cells overexpressing RNase H, a condition resulting in genome-wide loss of R-loops . Amongst the genes that were differentially expressed, we found that 64 and 48% of all downregulated and upregulated genes, respectively, displayed R-loops overlapping with 5hmC ( Figure 6A). Pathway analysis revealed that these differentially expressed genes (Supplementary file 1) are involved in the mechanistic target of rapamycin (mTOR) (downregulated) and MYC (upregulated) signaling pathways ( Figure 6B and C). mTOR and MYC are known to play dTET1 to the last exon of APOE in U-2 OS cells. R-loop data were normalized against RNase H-treated samples. *p<0.05, **p<0.01, two-tailed Student's t-test. (G) APOE transcription levels upon targeting dCas9-TET1 or dCas9-dTET1 to the last exon of the gene in U-2 OS cells. Data shown are the mean and SD from at least three independent experiments.
The online version of this article includes the following figure supplement(s) for figure 2: Figure supplement 1. Impact of Tet depletion in mouse embryonic stem (mES) cells in R-loops, 5-hydroxymethylcytosine (5hmC), and 5-methylcytosine (5mC).   opposite roles in establishing diapause, the temporary suspension of embryonic development driven by adverse environmental conditions (Fenelon et al., 2014), a stage that ES cells mimic when cultured in vitro. mTOR, a major nutrient sensor, acts as a rheostat during ES cell differentiation and reductions in mTOR activity trigger diapause (Bulut-Karslioglu et al., 2016). While overexpression of RNase H in mES cells did not reveal any significant changes in the cell cycle progression (Figure 6-figure supplement 1A and B), we observed a significantly decreased expression of genes related to pluripotency (Pou5f1) and germ layer commitment (Sox17, Sox6, Dll1) pathways ( Figure 6D). These data support the view that R-loops formed upon TET epigenetic reprogramming regulate gene expression in stem cells.

Discussion
In this study, we probed the hypothesis that 5hmC facilitates the co-transcriptional formation of noncanonical DNA secondary structures, known as R-loops. Data from in vitro transcription reactions and AFM provide direct evidence showing that transcription through 5hmC-rich DNA favors R-loop formation. By depleting TET enzymes in mES cells and fibroblasts, we demonstrate that TET activity increases cellular R-loop levels. Notably, the diminished levels of R-loops observed in TET-depleted cells did not result from impaired transcription, suggesting that 5hmC directly promotes R-loop formation. In agreement, tethering TET enzymes to a specific genomic locus using a CRISPR/Cas9-based system increases the levels of R-loops at the target locus. As 5hmC is mostly absent from the TSS, other chromatin and DNA features (e.g., histone modifications, DNA-supercoiling or G-quadruplex structures; García-Muse and Aguilera, 2019) known to induce R-loop formation are likely to operate in these regions. In contrast, the robust overlap between R-loops and 5hmC at the TTS of active genes suggests a putative causal link. Mechanistically, 5hmC may impact R-loop formation by either destabilizing the DNA duplex or altering RNA polymerase II elongation rate. Indeed, 5hmC modifies the DNA helix structure by favoring DNA-end breathing motion, diminishes the thermodynamic stability of the DNA duplex, and destabilizes GC pairing (Mendonca et al., 2014;Wanunu et al., 2011). It also weakens the interaction between DNA and nucleosomal histones (Mendonca et al., 2014), which is thought to accelerate RNA polymerase II elongation but can also facilitate nascent RNA annealing with the template DNA strand favoring R-loop formation. Future studies will clarify which one of these mechanisms, if not all, contribute to the observed impact of 5hmC on R-loops.
Acting as a promoter of R-loops, well-established drivers of DNA damage (García-Muse and Aguilera, 2019), 5hmC may indirectly harm genome integrity. Indeed, we found that 5hmC-rich loci are hotspots for DNA damage genome-wide. While such unscheduled R-loops formed at 5hmC-rich loci may threaten genomic integrity, regulated formation of R-loops at specific 5hmC-decorated loci may exert important regulatory roles. Indeed, R-loops play diverse physiological functions (García-Muse and Aguilera, 2019), such as the regulation of gene expression. Our findings that genomewide 5hmC and R-loops overlap robustly at the TTS of active genes and that TET-deficiency drives transcription readthrough support a model whereby TET enzymes act upstream of R-loop formation during transcription termination (Skourti-Stathaki et al., 2011).
TETs may play dual roles as both oncogenic and tumour suppressor genes, with the former arising as the consequence of altered expression levels or function, as observed in several cancers, such as triple-negative breast cancer (Bray et al., 2021;Good et al., 2018). In addition to altering the distribution along the Acadl, Cr1l, and Srr genes. Density signals are represented as reads per kilobase (RPKMs). (F) Metagene profiles of 5hmC and R-loops distribution in active genes. The gene body region was scaled to 60 equally sized bins, and ±10 kbp gene-flanking regions were averaged in 200 bp windows. TSS: transcription start site; TTS: transcription termination site. Density signals are represented as RPKMs, and error bars (gray) represent standard error of the mean. (G) Metagene profiles of genes showing transcription readthrough in wild-type and TET1 KO human ES cells. All gene regions were scaled to 2000 bp (gene body) and divided in equal bins of 100 bp. 1000 bp regions averaged in 100 bp bins were added upstream the TSS and downstream the TTS region. *p<0.05, Mann-Whitney rank test.
The online version of this article includes the following figure supplement(s) for figure 3:    expression levels of tumour suppressors or oncogenes (Bray et al., 2021), our findings suggest that TET-driven changes in the DNA methylation landscape may as well drive transcription-dependent genome damaging events that could facilitate cancer development and progression. In agreement with this view, a TET1 isoform that lacks regulatory domains, including its DNA-binding domain, but retains its catalytic activity, is enriched in cancer cells (Good et al., 2017), suggesting that mistargeted TET activity may drive oncogenic events, such as genomic instability. Conversely, TET activity deposits 5hmC at DNA damage sites induced by aphidicolin or microirradiation in HeLa cells and prevents chromosome segregation defects in response to replication stress (Kafer et al., 2016). While the role that TETs play during carcinogenesis is not yet clear, the impact of 5hmC on stem cell differentiation and development has been extensively studied (Ficz et al., 2011). By driving the developmental DNA methylome reprogramming, TETs carry out numerous functions related to early developmental processes. Here, we disclose a putative new role for R-loops as mediators of 5hmC-driven gene expression programs in stem cells. Our gene ontology analysis revealed that R-loops formed at 5hmCrich regions impact the expression of genes involved in establishing diapause. This stage of temporary suspension of embryonic development is triggered by adverse environmental conditions (Fenelon et al., 2014). Accordingly, changes in the activity of mTOR, a major nutrient sensor, control ES cell commitment to trigger diapause (Bulut-Karslioglu et al., 2016). The mTOR signaling pathway was significantly downregulated upon global R-loop suppression by RNase H. Conversely, MYC targets,  which prevent ES cells from entering the state of dormancy that characterizes diapause (Scognamiglio et al., 2016), were amongst the genes more significantly upregulated upon RNase H overexpression in mES cells. These cells also displayed reduced expression levels of genes related to pluripotency (Pou5f1) and germ layer commitment (Sox17, Sox6, Dll1) pathways. Whether the controlled 5hmCdriven formation of R-loops at specific genes drives stem cells fate and how do TET enzymes capture the environmental cues to target R-loop formation at selected genes are important questions that emerge from our findings. Our study sets the ground for further research aimed at investigating the role of R-loops in ES cells.

RNA isolation and quantitative RT-PCR
Total RNA was isolated using TRIzol reagent (15596018, Invitrogen). cDNA was prepared through reverse transcriptase activity (MB125, NZYTech). RT-qPCR was performed in the ViiA 7 Real-Time PCR system (Applied Biosystems) using PowerUp SYBR Green Master Mix (A25918, Applied Biosystems). Relative RNA expression was estimated as follows: 2 (Ct reference -Ct sample) , where Ct reference and Ct sample are mean threshold cycles of RT-qPCR done in duplicate for U6 snRNA or Gapdh mRNA and for the gene of interest, respectively. Primer sequences are presented in Supplementary file 2B.

Dot blot of genomic R-loops, 5mC, and 5hmC
Cells were lysed in lysis buffer (100 mM NaCl, 10 mM Tris pH 8.0, 25 mM EDTA pH 8.0, 0,5% SDS, 50 μg/mL Proteinase K) overnight at 37°C. Nucleic acids were extracted using standard phenolchloroform extraction protocol and resuspended in DNase/RNase-free water. Nucleic acids were then fragmented using a restriction enzyme cocktail (20 U each of EcoRI, BamHI, HindIII, BsrgI, and XhoI). Half of the sample was digested with 40 U RNase H (MB085, NZYTech) for 48 hr at 37°C to be used as a negative control in R-loops blotting. Digested nucleic acids were cleaned with standard phenolchloroform extraction and resuspended in DNase/RNase-free water. Nucleic acids samples were quantified in a NanoDrop 2000 spectrophotometer (Thermo Scientific), and equal amounts of DNA were deposited into a positively charged nylon membrane (RPN203B, GE Healthcare). Membranes were UV-crosslinked using UV Stratalinker 2400 (Stratagene), blocked in 5% (m/v) milk in PBSt (PBS 1× containing 0.05% [v/v] Tween 20) for 1 hr at room temperature, and immunoblotted with specific antibodies. For the loading control, membranes were stripped in 0.5% SDS for 1 hr at 60°C, followed by blocking and re-probing. Details of the antibodies used are included in Supplementary file 2C.

Proximity ligation assay (PLA)
E14 mES cells were grown on coverslips and fixed/permeabilized with methanol for 10 min on ice, followed by 1 min acetone on ice. Cells were then incubated with primary antibodies for 1 hr at 37°C, followed by a pre-mixed solution of PLA probe anti-mouse minus (DUO92004, Sigma-Aldrich) and PLA probe anti-rabbit plus (DUO92002, Sigma-Aldrich) for 1 hr at 37°C. Localized rolling circle amplification was performed using Detection Reagents Red (DUO92008, Sigma-Aldrich), according to the manufacturer's instructions. Slides were mounted in 1:1000 DAPI in Vectashield. For the RNase H control, fixed cells were treated with 3 U/μL RNase H (MB085, NZYTech) for 1 hr at 37°C prior to incubation with the antibodies. Images were acquired using the Point Scanning Confocal Microscope Zeiss LSM 880, 63×/1.4 oil immersion, with stacking acquisition and generation of maximum intensity projection images. PLA foci per nucleus were quantified using ImageJ. Details of the antibodies used are mentioned in Supplementary file 2C. In vitro transcription PCR products were subject to in vitro transcription using the HiScribe T7 High Yield RNA Synthesis Kit (E2040S, NEB), which relies on the T7 RNA polymerase to initiate transcription from a T7 promoter sequence (present in our fragments). Reactions were performed for 2 hr at 37°C, using 1 μg of DNA as template, according to the manufacturer's instructions. The resulting RNA was column-purified with NucleoSpin RNA isolation kit (740955.250, Macherey-Nagel) and quantified in a NanoDrop 2000 spectrophotometer (Thermo Scientific).

Dot blot of R-loops formed in in vitro
Half of each in vitro transcription product was treated with 10 U RNase H (MB085, NZYTech) at 37°C overnight to serve as negative control. Then, all samples were treated with 0.05 U RNase A (10109142001, Roche) at 350 mM salt concentration for 15 min at 37°C and ran on agarose gel. Nucleic acids were transferred overnight to a nylon membrane through capillary transfer. The membrane was then UV-crosslinked twice, blocked in 5% milk in PBSt for 1 hr at room temperature, and incubated with the primary antibody at 4°C overnight. Signal quantification was performed using ImageJ. Details of the antibodies used are included in Supplementary file 2C.

Atomic force microscopy
RNase A-treated in vitro transcription products, treated or not with RNase H, were purified through phenol-chloroform extraction method and resuspended in DNase/RNase-free water. DNA solution was diluted 1:10 in Sigma ultrapure water (with final 10 mM MgCl 2 ) and briefly mixed to ensure even dispersal in solution. A 10 μL droplet was deposited at the center of a freshly cleaved mica disc, ensuring that the pipette tip did not contact the mica substrate. The solution was let to adsorb on mica surface for 1-2 min to ensure adequate coverage. The mica surface was carefully rinsed with Sigma ultrapure water, so that excess of poorly bound DNA to mica is removed from the mica substrate. Afterward, the mica substrate was dried under a gentle stream of argon gas for approximately 2 min, making sure that any excess water is removed. DNA imaging was performed using a JPK Nanowizard IV atomic force microscope, mounted on a Zeiss Axiovert 200 inverted optical microscope. Measurements were carried out in tapping mode using commercially available ACT cantilevers (AppNano). After selecting a region of interest, the DNA was scanned in air, with scan rates between 0.5 and 0.9 Hz. The setpoint selected was close to 0.3 V. Several images from different areas of the same sample were performed and at least three independent samples for each condition were imaged. All images were of 512 × 512 pixels and analyzed with JPK data processing software.

DNA:RNA immunoprecipitation (DRIP)
Cells were collected and lysed in 100 mM NaCl, 10 mM Tris pH 8.0, 25 mM EDTA, 0.5% SDS, 50 μg/ mL Proteinase K overnight at 37°C. Nucleic acids were extracted using standard phenol-chloroform extraction protocol and resuspended in DNase/RNase-free water. Nucleic acids were then fragmented using a restriction enzyme cocktail (20 U each of EcoRI, BamHI, HindIII, BsrgI, and XhoI), and 10% of the digested sample was kept aside to use later as input. Half of the remaining volume was digested with 40 U RNase H (MB085, NZYTech) to serve as negative control, for 72 hr at 37°C. Digested nucleic acids were cleaned with standard phenol-chloroform extraction and resuspended in DNase/ RNase-free water. DNA:RNA hybrids were immunoprecipitated from total nucleic acids using 5 µg of S9.6 antibody (MABE1095, Merck Millipore) in binding buffer (10 mM Na 2 HPO 4 pH 7.0, 140 mM NaCl, 0.05% Triton X-100), overnight at 4°C. 50 µL protein G magnetic beads (10004D, Invitrogen) were used to pull down the immune complexes at 4°C for 2-3 hr. Isolated complexes were washed five times (for 1 min on ice) with binding buffer and once with Tris-EDTA (TE) buffer (10 mM Tris pH 8.1, 1 mM EDTA). Elution was performed in two steps, for 15 min at 55°C each, using elution buffer (50 mM Tris pH 8.0, 10 mM EDTA, 0.5% SDS, 60 µg/mL Proteinase K). The relative occupancy of DNA:RNA hybrids was estimated by RT-qPCR as follows: 2 (Ct Input−Ct IP) , where Ct Input and Ct IP are mean threshold cycles of RT-qPCR done in duplicate for input samples and specific immunoprecipitations, respectively. Data were normalized against the corresponding RNase H-treated samples and plotted as absolute numbers or as fold change over control. Primer sequences are shown in Supplementary file 2B.

5-(Hydroxy)methylated DNA immunoprecipitation ((h)MeDIP)
Cells were collected and lysed in 100 mM NaCl, 10 mM Tris pH 8.0, 25 mM EDTA, 0.5% SDS, 50 μg/ mL Proteinase K overnight at 37°C. Samples were sonicated with four pulses of 15 s at 10 mA intensity using a Soniprep150 sonicator (keeping tubes for at least 1 min on ice between pulses). Fragmented nucleic acids were cleaned with standard phenol-chloroform extraction protocol and resuspended in DNase/RNase-free water. 10% of sample was kept aside to use later as input. The remaining volume was denatured by boiling the samples at 100°C for 10 min, followed by immediate chilling on ice and quick spin. Samples were divided in half, and 5 µg of anti-5mC antibody (61255, Active Motif) or 5 µg of anti-5hmC antibody (39791, Active Motif) were used to immunoprecipitate 5mC and 5hmC, respectively, in binding buffer (10 mM Na 2 HPO 4 pH 7.0, 140 mM NaCl, 0.05% Triton X-100), overnight at 4°C. 50 µL protein G magnetic beads (10004D, Invitrogen) were used to pull down the immune complexes at 4°C for 2-3 hr. Isolated complexes were washed five times (for 1 min on ice) with binding buffer and once with TE buffer (10 mM Tris pH 8.1, 1 mM EDTA). Elution was performed in two steps, for 15 min at 55°C each, using elution buffer (50 mM Tris pH 8.0, 10 mM EDTA, 0.5% SDS, 60 µg/mL Proteinase K). The relative occupancy of 5mC and 5hmC was estimated by RT-qPCR as follows: 2 (Ct Input−Ct IP) , where Ct Input and Ct IP are mean threshold cycles of RT-qPCR done in duplicate for input samples and specific immunoprecipitations, respectively. Primer sequences are presented in Supplementary file 2B.
Cell cycle analysis pEGFP-N1 (GFP coding plasmid used as control) was purchased from Addgene, and pEGFP-RNaseH1 (GFP-tagged RNase H1 coding plasmid) was kindly provided by Robert J. Crouch (NIH, USA). Seeded mES cells were transfected with GFP (control) or GFP-tagged RNase H coding plasmids using Lipofectamine 3000 Transfection Reagent (L3000015, Invitrogen). 24 or 48 hr later, cells were trypsinized and pelleted by centrifugation at 500 × g for 5 min. Cells were fixed in cold 1% PFA for 20 min at 4°C, followed by permeabilization in 70% ethanol for 1 hr at 4°C. Cells were then treated with 25 μg/ mL RNase A (10109142001, Roche) in PBS 1× at 37°C for 20 min, followed by staining with 20 μg/mL propidium iodide (P4864, Sigma-Aldrich) in PBS 1× for 10 min at 4°C. Flow cytometry was performed on a BD Accuri C6 (BD Biosciences), and data were analyzed using FlowJo software.

Electrophoretic mobility shift assay (EMSA)
DNA:RNA hybrids formed with either C-, 5hmC-, or 5mC-containing DNA were obtained by incubating ssDNA with the complementary ssRNA in annealing buffer (100 mM KAc, 30 mM HEPES pH 7.5). Native and C-modified oligonucleotides were ordered from IDT (Supplementary file 2E). Hybrid formation was confirmed in a native polyacrylamide gel. Increasing amounts of S9.6 antibody (MABE1095, Merck Millipore) were added to the DNA:RNA hybrids, and the complexes were ran in a native polyacrylamide gel to assess the S9.6 capacity to bind hybrids containing each of the three C variants. The amount of free probe was quantified using ImageJ.

5hmC, R-loop, and ɣH2AX genome-wide characterization
The HTS datasets produced by immunoprecipitation (DRIP-seq, ChIP-seq, and hMeDIP-seq) were analyzed through the same workflow. First, the reads were aligned to the reference mouse and human genome (mm10 and GRCh38/hg38 assemblies, respectively) with Bowtie (Langmead et al., 2009), and filtering for uniquely aligned reads. Enriched regions were identified relative to the input samples using MACS (Zhang et al., 2008), with a false discovery rate of 0.05. Finally, enriched regions were assigned to annotated genes, including a 4-kilobase region upstream the TSS and downstream the TTS. Gene annotations were obtained from mouse and human GENCODE annotations (M11 and v23 versions, respectively) and merged into a single transcript model per gene using BEDTools (Quinlan and Hall, 2010). For individual and metaprofiles, uniquely mapped reads were extended in the 3′ direction to reach 150 nt with the Pyicos (Althammer et al., 2011). Individual profiles were produced using a 20 bp window. For the metaprofiles centered around 5hmC peaks, 5hmc-enriched regions were aligned by the peak summit (maximum of the peak) and the read density for the flanking 10 kbp were averaged in a 200 bp window. For the metagene profiles, the gene body region was scaled to 60 equally sized bins and ±10 kbp gene-flanking regions were averaged in 200 bp windows. All profiles were plotted as normalized reads per kilobase per million mapped reads (RPKMs). A set of in-house scripts for data processing and graphical visualization were written in bash and in the R environmental language (https://www.r-project.org/; R Core Team, 2018). SAMtools (Li et al., 2009) and BEDTools were used for alignment manipulation, filtering steps, file format conversion, and comparison of genomic features. Statistical significance of the overlap between 5hmC and R-loops was assessed with enriched regions and permutation analysis. Briefly, random 5hmC and R-loops-enriched regions were generated 1000 times from annotated genes using the shuffle BEDTools function (maintaining the number and length of the originally datasets). The p-value was determined as the frequency of overlapping regions between the random datasets as extreme as the observed.

Transcriptome analysis
Expression levels (transcripts per million [TPMs]) from RNA-seq and GRO-seq datasets were obtained using Kallisto (Bray et al., 2016), where reads were pseudo-aligned to mouse and human GENCODE transcriptomes (M11 and v23, respectively). Transcriptionally active genes for 5hmC and R-loops annotation were defined as those with expression levels higher than the 25th percentile. Differential expression in mES cells overexpressing RNase H was assessed using edgeR (v3.20.9) and limma (v3.34.9) R packages (Robinson et al., 2010;Ritchie et al., 2015). Briefly, sample comparison was performed using voom transformed values, linear modeling, and moderated t-test as implemented in limma R package, selecting significantly differentially expressed genes with B-statistics higher than zero. Significantly enriched pathways of up-and downregulated genes (with overlapping R-loops/5hmC regions) were selected using Fisher's exact test and all expressed genes as background gene list. Evaluated pathways were obtained from the hallmark gene sets of Molecular Signatures Database (MSigDB) (Liberzon et al., 2015) and filtered using false discovery rate-corrected p-values<0.05.
For the analysis of transcription readthrough, transcriptome profiles from human embryonic stem cells (WT and TET1 KO) were obtained from a GEO (GSE169209). RNA-seq data were mapped to the reference human genome (GRCh38) with the STAR v2.7.8a using default parameters (Dobin et al., 2013). Transcription readthrough levels were evaluated by counting the number of reads mapping downstream the TTS using ARTDeco (Roth et al., 2020) and human genome annotation from GENCODE project (GENCODE release 37). Genes with an enrichment in transcriptional readthrough in TET1 KO samples relative to the control were identified. Metagene profiles were built using the computeMatrix tool from the deepTools v3.5.1 (Ramírez et al., 2016) and default packages from Python language. Genes were scaled to equally sized bins of 100 bp so that all annotated TSSs and TTSs were aligned. Regions of 1 kb were added upstream of TSS and downstream of TTS and also averaged in 100 bp bins. All read counts were normalized by the number of mapped reads (RPKM). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Additional files • Transparent reporting form

Author contributions
• Source data 1. Original, uncropped images of blots and gels.

Data availability
All data generated or analysed during this study are included in the manuscript and supporting files.
The following previously published datasets were used: Continued on next page