The chromatin network helps prevent cancer-associated mutagenesis at transcription-replication conflicts

Article https://doi.org/10.1038/s41467-023-42653-0 The chromatin network helps prevent cancer-associated mutagenesis at transcription-replication conflicts Aleix Bayona-Feliu 1,2,3 , Emilia Herrera-Moyano 1,2 , Nibal Badra-Fajardo 1 , Iván Galván-Femenía 3 , María Eugenia Soler-Oliva 1,2 &AndrésAguilera 1,2 Genome instability is a feature of cancer cells, transcription being an important source of DNA damage. This is in large part associated with R-loops, which hamper replication, especially at head-on transcription-replication conflicts (TRCs). Here we show that TRCs trigger a DNA Damage Response (DDR) involving the chromatin network to prevent genome instability. Depletion of the key chromatin factors INO80, SMARCA5 and MTA2 results in TRCs, fork stalling and R-loop-mediated DNA damage which mostly accumulates at S/G2, while histone H3 Ser10 phosphorylation, a mark of chromatin compaction, is enriched at TRCs. Strikingly, TRC regions show increased mutagenesis in cancer cells with signatures of homologous recombination deficiency, transcription-coupled nucleotide excision repair (TC-NER) and of the AID/ APOBEC cytidine deaminases, being predominant at head-on collisions. Thus, our results support that the chromatin network prevents R-loops and TRCs from genomic instability and mutagenic signatures frequently associated with cancer. Cancer is a multifactorial disease triggered by the dysfunction of several cellular processes that fuel uncontrolled cell proliferation. Genome instability (GIN) is a hallmark of cancer cells and one of the main causes of cell transformation and cancer evolution1.Unstablegenomes favor mutagenesis of tumor suppressors and oncogenes thus promoting malignant cell transformation and facilitating tumor adaptation to a wide range of scenarios, including therapeutic treatments. Remarkably, chromatin-regulating enzymes are emerging as key factors in maintaining genome integrity2,andfine-tuning of chromatin structure and function has been linked to prevent transcriptionassociated DNA damage thus ensuring genome fidelity. Consistently, chromatin-modulating activities are frequently found altered in tumors3. The transcription machinery may pose a roadblock to replication fork progression, promoting replication stress and DNA damage4that can potentially be enhanced by the occurrence of non-B DNA structures such as R-loops, a three-stranded nucleic acid structure composed by a DNA–RNA hybrid and a displaced single strand DNA (ssDNA) that form co-transcriptionally. Although R-loops have physiological roles in particular processes such as class switch recombination, they may pose a threat to genome integrity. Consequently, cells possess several mechanisms to counteract the pathological accumulation of R-loops, either by preventing their formation via RNA binding or processing factors, such as the THO complex, SRSF1 and TDP43 among others5–10, or by resolving factors like DNA:RNA helicases such as SETX, UAP56/DDX39b or DDX511–15 and RNases H, such as RNH1 or RNH216. In addition, R-loops may indirectly be resolved by the action of DDR factors, such as BRCA1, BRCA2, FANCD2 or ATR, among others17–21. The Encyclopedia of DNA Elements (ENCODE) project22,23 allows access to a large number of genome-wide data for comparative studies, in particular from K562 cells, used mostly as a standard for Received: 7 February 2023 Accepted: 9 October 2023 Check for updates 1 Centro Andaluz de Biología Molecular y Medicina Regenerativa CABIMER, Universidad de Sevilla-CSIC-Universidad Pablo de Olavide, 41092 Seville, Spain. 2 Departamento de Genética, Facultad de Biología, Universidad de Sevilla, 41012 Seville, Spain. 3 Institute for Research in Biomedicine (IRB Barcelona), The Barcelona Institute of Science and Technology (BIST), Barcelona, Spain. e-mail: [email protected];[email protected] Nature Communications | (2023) 14:6890 1 1234567890():,; 1234567890():,; genome-wide studies. Furthermore, cancer-associated mutagenesis is also accessible through the Catalogue of Somatic Mutations in Cancer (COSMIC) database24,enablingidentification of mutational signatures associated with a specific etiology. Recently, the availability of DRIPcseq (DNA:RNA hybrid immunoprecipitation followed by sequencing of the cDNA derived from the RNA moiety of hybrids) and OK-seq data (Okazaki fragment sequencing)25 in K562 cells allowed us to predict genomic sites prone to TRCs26. Here, we show that in addition to the SWI/SNF complex, other chromatin remodeler factors such as INO80, SMARCA5 and MTA2 contribute to prevent R-loop accumulation in cells. Then, we used the ENCODE genome-wide data in K562 cells and COSMIC database to cross ChIP-seq and cancer-associated mutagenesis with our genomic TRC dataset to unveil possible structural and functional features of TRC sites. We show that a wide range of chromatin remodelers, modifiers and epigenetic marks, as well as transcription and DDR factors are enriched at TRC sites, in particular at those occurring in head-on orientation. This correlates with a high incidence of single nucleotide variations (SNVs) and insertions and deletions (indels) mutations at such sites in tumoral cells, as found in the COSMIC database. Strikingly, such mutagenic signatures are found highly enriched at sites where R-loops and head-on TRCs occur preferentially in cells depleted of SWI/SNF.Altogether,these results establish a direct link between the DDR and epigenetic factors at TRCs and their contribution to prevent mutagenesis and genomic instability mediated by DSB repair and translesion synthesis (TLS) associated with fork stalling that is over-represented in cancer cells. This study may open new perspectives in the design and development of new cancer therapies. Results Enrichment of chromatin, DNA, and RNA metabolism factors at TRCs Different studies have addressed the contribution of DNA repair factors on DNA–RNA hybrid accumulation and its impact on genome integrity. The results support the view that DNA repair factors regulate DNA–RNA hybrid homeostasis upon fork stalling, DSBs formation and post-replicative repair (PRR)17–21,27,28. The relevance of chromatin to this process has been provided by different studies26,29–32, but it seems yet incomplete considering the high number of factors influencing chromatin function. Indeed, the frequency of mutations in chromatin remodeler and modifier genes is high in cancer, genome instability being a hallmark of tumor cells. Accordingly, the SWI/SNF complex, one of the most frequently altered chromatin complexes in cancer, plays a key role helping solve TRCs26,30. Consequently, we did a global search for additional chromatin and functionally interacting factors. Despite the large number of genes demonstrated or proposed to regulate R-loop homeostasis and therefore to impact on TRCs and genome integrity, our knowledge on specific factors associated with TRC sites is still limited. The ENCODE Project provides massive, curated genome-wide data that can be used for this purpose, particularly regarding factors involved in chromatin organization. Thus, we analyzed the abundance of the whole catalog of proteins for which ChIP-seq data are available in ENCODE at R-loop- enriched TRC sites to unveil the factors and molecular processes contributing to preserve genome integrity at TRCs. We first gathered the ENCODE ChIP-seq data available for the K562 cell line and integrated them with the R-loop sites obtained by DRIPc-seq, where replication fork directionality (RFD) is mostly homogenous among cell population (RFD > |0.75|), to ensure accurate and directional TRC analysis as previously shown26. As control, we also included FANCD2 ChIP-seq data33 to the analysis, a standard marker for replication fork stalling34 previously described to accumulate at these sites26. Then, we plotted the coverage of proteins along R-loop-enriched TRC sites according to the RFD and measured protein abundance at TRCs (Fig. 1a). We extended the plots to +/−1MbfrompeaksastheDNA damage foci highlighted by γH2AXmayextendupto1–2Mb in mammals and we were analyzing DDR factors. Next, we performed gene set enrichment analysis (GSEA) to determine whether specific chromatin functional categories could be retrieved either enriched or depleted at the R-loop-rich TRC sites. To do so, we selected a subset of Gene Ontology (GO) functional categories including chromatin, but also transcription, DNA replication and DNA damage GO termsas controls, and we searched for significant enrichments using the GSEA software (UC San Diego and Broad Institute)35,36. First, we screened for factors whose abundance was significantly changed at the R-loop-enriched TRC sites in respect of the surrounding region by filtering those factors with Rank Metric Scores (RMS) > | 0.25| at TRCs. Our analysis revealed 209 factors that were enriched and 12 depleted (Fig. 1b). In accordance with published data, several factors previously described to participate in R-loop metabolism including SMARCA4, FANCD2, histone acetylation, BRD4, EP400, YY1 and EWSR1 were found enriched at TRCs17,18,26,32,37–40.Consistently, proteins encompassing transcription by RNA polymerase II (RNAPII), as well as factors involved in the DNA replication and DNA Damage Response (DDR) were also found enriched at TRC sites (Fig. 1c, Supplementary Fig. 1a–f, and Supplementary Table 1). Indeed, the FANCD2 Fanconi Anemia factor, the RAD51 recombinase (RAD51), the nibrin (NBN) component of the MRN complex and the RAD21 cohesin component emerged in our analysis as DNA repair factors peaking at TRCs. Strikingly however, canonical DNA replication components were observed dropping significantly at these regions (Fig. 1c, Supplementary Fig. 1e, and Supplementary Table 2), with a significant reduction in replicative helicases as indicated by the helicase GO term. These results are consistent with fork collapse and DNA breakage at TRCs as largely reported17,18,21,28,41, suggesting that these events are concomitant to disengaging of the replicative MCM helicase followed by the action of DSB repair factors. In addition to the above-mentioned factors harboring a functional role during DNA and RNA metabolism, an important number of epigenetic factors were also found associated with TRC sites, including chromatin remodelers and modifiers as well as specificepigenetic marks (Fig. 1d, Supplementary Fig. 1d, and Supplementary Table 1). Thus, components of the 4 main chromatin remodeling families (SWI/ SNF, ISWI, INO80 and CHD) were enriched at TRC sites. More specifically, these comprised the factors of the SWItch/Sucrose Non Fermentable (SWI/SNF) complex SMARCA4, SMARCC2/BAF170, ARID1B and ARID2, the members of the Inositol 80 (INO80) family EP400 and YY1, the component of the Imitation SWItch (ISWI) complex SMARCA5, and the MTA members of the Nucleosome Remodeling Deacetylase (NuRD) complex, which belongs to the Chromodomain helicase DNA- binding (CHD) remodeling family. In addition, chromatin modifiers mostly involved in histone acetylation and methylation were also overrepresented at TRCs (Fig. 1d, Supplementary Fig. 1d, and Supplementary Table 1). The first class included the HDAC1, HDAC3 and SIN3B histone deacetylase complexes (HDACs), the GATAD2B component of the MeCP1 HDAC complex, as well as the bromo-domain factor BRD4, whereas the second included the PHF8, KDM1A histone lysine demethylase (KDM) and the L3MBTL2 histone methyl-lysine binding protein of the Polycomb group (PcG). Consistently, histones H3 and H4 monomethylated H3K9me1 and H4K20me1, di-methylated H3K4me2 and H3K79me2, and tri-methylated H3K4me3 and H3K36me3, as well as acetylated H3K9ac and H3K27ac were increased at TRC sites. Moreover, components regulating DNA methylation and higher order chromatin structure such as DNMT1 and TRIP13 or FOXA1, respectively, were also found enriched at TRCs. Interestingly, similar to the underrepresentation of the replicative MCM helicases, some specific epigenetic marks and factors were significantly reduced at TRCs (Fig. 1c, Supplementary Fig. 1d, and Supplementary Table 2). These included the H3K27me3 epigenetic mark, and the CBX8 chromo-box member of the PcG PRC1-like complex, the Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 2 CTCFL testis-specific insulator binding factor or the PRMT5 arginine N-methyltransferase 5, among others. Next, we investigated whether the role of the identified factors was more relevant in head-on versus co-directional TRCs. A detailed comparison of factor abundance revealed that 196 showed higher enrichment at head-on collisions, while only 5 were more abundant at co-directional TRCs (Supplementary Fig. 2a, b and Supplementary Tables 3 and 4). Factors enriched at head-on TRCs covered all GO categories previously identified (Supplementary Fig. 2b and Supplementary Table 3). In contrast, TRC-depleted proteins were found more abundant in co-directional TRCs (Supplementary Fig. 2c, d and Supplementary Table 4). These included the replicative helicases MCM2/3/5/7 and factors involved in chromatin organization like H3K27me3, PRMT5, CTCFL or SUZ12. Altogether, our data supports the view that TRCs, besides triggering a strong DDR associated with changes in the presence of replication and DSB repair factors, are linked to chromatin modification and remodeling, especially at head-on TRCs. Replication might be PKNOX1 -4-224 Upstream Enrichment Score Downstream Enrichment Score -3 -2 -1 1 2 3 MCM3 H3K27me3 MCM5 MCM2 MCM7 SUZ12 NRF1 H3K79me2 GMEB1 SUPT5H PRDM10 SM A RCA4 FANCD2 En r iched facto r sDepleted f acto r s 209 12 DRIPc-seq R - l o o p OK-seq TRCs TRC aEnrichment Score bEnrichment Score nEnrichment Score RF fa c to r a f a c to r b f a c to r n ... A TP-dependent chromatin remodelling RN A polyme r ase complex Histone modifications Site of DN A damage Replication fo r k processing Transc r iption by RN A polyme r ase II Transc r iption termination Chromatin organization DN A damage checkpoint signalling DN A replication checkpoint signalling Replication fo r k protection Chromatin Chromatin remodelling RN A binding Transc r iption initiation DN A replication Double st r and break r epair Replication fo r k DN A repai r Transc r itption elongation MCM complex DN A -helicase activity TRC-En r iched TRC-depleted F actors per GO term (%) 0 0 8 0 0 1 0 4 0 6 20 a bc Upstream Downstream MOST RELEVANT TRC-ENRICHED EPIGENETIC FACTORS Histone acetylation SMARCA4 ARID1B ARID2 SMARCC2 DPF2 HLTF IKZF1 NR3C1 SMARCA5 SWI/SNF family ISWI family INO80 family CHD family (NuRD complex) YY1 NFRKB EP400 MTA1 MTA2 MTA3 TRIM28 GATAD2B IKZF1 KDM1A ASH1L PHF8 L3MBTL2 DNMT1 TRIP13 FOXA1 Histone acetylation Histone methylation DNA methylation SIN3B HDAC1 HDAC3 KAT8 BRD4 NCOA6 NCOR1 TRIM24 HCFC1 Higher-order structures Chromatin Remodelling Chromatin Modification H3K79me2 H4K20me1 H3K4me2 H3K4me3 H3K9me1 H3K36me3 H3K9ac H3K27ac Histone PTMs Histone methylation d R - l o o p RF ENCODE Project Fig. 1 | Factors enriched and depleted at transcription-replication conflicts (TRCs). a Schematic summary of the ChIP-seq data downstream processing for enrichment/depletion analysis at TRCs. bScatter plot indicating the fold change in content for each factor when comparing the TRC site versus its upstream and downstream sequences. cPercentage of factors enriched and depleted at TRCs for each of the GO categories analyzed. dSchematic chart with the most relevant chromatin factors found enriched at TRCs. Bold text is used to highlight gene subcategories. TRC transcription-replication conflict, RF replication fork. The arrow under RF indicate replicationfork directionality. Source data are provided as aSourcedatafile. See also Supplementary Figs. 1 and 2. Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 3 hindered leading to alternative mechanisms as suggested by the decrease of MCM helicases. Histone 3 serine 10 phosphorylation extends ~1 Mb upstream of TRCs and concentrates in head-on orientation Our results support the idea that fine tuning of chromatin is key to prevent TRCs from resulting in DNA breaks, since chromatin structure needs to be dynamic to facilitate DNA replication and transcription and prevent DNA damage. In this context, harmful R-loops have been found associated with histone H3 Ser10 phosphorylation (H3S10pho) from yeast to mammals42, which suggest a functional link between H3S10pho and genetic instability caused by R-loops. This conclusion is supported by the observation that yeast histone H3 mutants unable to generate H3S10pho show high levels of R-loops not accompanied by an increase in DNA damage43. However, how R-loop-associated H3S10pho determines R-loop harmfulness is unclear. Given the strong abundance of chromatin-modulating factors at TRCs, particularly head-on TRCs, we compared H3S10pho at co-directional versus head-on TRCs. Since MYCN was recently shown to activate on chromatin Aurora- A, which phosphorylates H3S10pho in S phase and promotes R-loop suppression44, we used the available genome-wide data of H3S10pho profiling during S phase from this study to investigate the coverage of H3S10pho at TRCs. First, we compared the H3S10pho genomic distribution in the neuroblastoma cell line IMR-5 with the DRIPc-seq data in control K562 cells (Fig. 2a, b). Consistent with published data from yeast42,43, the genomic profile of H3S10pho during S phase correlated with that of R-loops in K562 cells. Metadata analysis showed a clear accumulation of H3S10pho at these sites: 1951 out of 2771 (70.4%) of H3S10pho-accumulating genes were also R-loop prone (Fig. 2c, d). Furthermore, a massive and asymmetric increase of H3S10pho was detected around head-on TRCs at the Mb scale that was not observed at co-directional TRCs (Fig. 2e). H3S10pho extended from ~0.2 Mb to ~1 Mb upstream the TRC, and certain accumulation was also appreciable by ~0.5 Mb downstream. This is consistent with H3S10pho foci accumulation detected in IF experiments under conditions triggering R-loop accumulation in human cells42. Notably, chromatin interactions are more intense at head-on rather than co-directional Fig. 2 | Evaluation of H3S10pho signal over TRCs. a Representative screenshot of a genome region showing co-localization of R-loops (blue) and H3S10pho (orange). bExamples of two R-loop-prone genes (RPL13A and H1-2) showing co-localization with H3S10pho. cMetapeak analysis. H3S10pho mean coverage around +/−5kbof R-loop peaks. dVenn diagram showing co-occurrence of R-loop-prone genes (DRIPc; blue) and H3S10pho (orange) mark in control K562 cells. eMetanalysis at megabase scale. H3S10pho mean coverage around +/−1Mbofhead-onandco- directional R-loop peaks. Replication fork directionality is indicated. fAnalysis of chromatin interactions around TRCs. Replication fork directionality is indicated. TRC transcription-replication conflict, HO head-on, CD co-directional, RF replication fork. Arrows under RF indicate replication fork directionality. Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 4 TRCs, especially upstream of the TRCs where H3S10pho accumulates, as determined by HiC data metanalysis (Fig. 2f), consistent with the broadly reported connection between H3S10pho and chromatin condensation. These results further confirm the connection between R-loops and H3S10pho genome-wide in human cells and reveal a major occurrence of the modification at head-on collisions. Depletion of TRC-enriched chromatin factors leads to R-loop- mediated DNA damage Alteration of epigenetic factors is common in cancer cells, and their deficiencies have been associated with genome instability, TRCs being a potential source of DNA damage and genome instability. Given that the GO term with the highest enrichment scores at TRCs was ATP- dependent chromatin remodelers (Fig. 1and Supplementary Fig. 1), we next assayed whether silencing the expression of components from the ISWI, INO80 and CHD families, whose impact on TRCs had not been analyzed previously had an impact on R-loop metabolism, as observed for the SWI/SNF family26. Thus, we analyzed the consequences of depleting SMARCA5 and INO80, the core subunits of chromatin remodelers ISWI and INO80, and MTA2, a subunit of the Nucleosome Remodeling and Deacetylase (NuRD) complex that belongs to the CHD remodeling family, on genome integrity. We depleted SMARCA5, INO80 and MTA2 in HeLa cells via small interference RNA (siRNA) that caused >80% depletion efficiency as determined by Western Blot (WB) (Supplementary Fig. 3a, b). Then, we analyzed DNA break accumulation by immunofluorescence (IF) using the γH2AX antibody as a marker. In most cases, gene silencing resulted in an accumulation of γH2AX foci consistent with an increase in DNA breaks (Fig. 3a and Supplementary Fig. 3c, d), as reported in siSMARCA4 cells26. The percentage (%) of cells with more than 5 γH2AX foci raised significantly by ~1.6× in siSMARCA5, siINO80 and siMTA2 cells compared to siRNA control cells (siC). Next, to unveil whether the increase on DNA damage formation was due to abnormal R-loop metabolism, we overexpressed RNase H1 (RNH1), which degrades the RNA moiety of DNA-RNA hybrids, and measured the impact on DNA break occurrence. The results showed a significant decrease in the % of cells containing more than 5 γH2AX foci, close to the levels of siC cells in all the three conditions showing increased DNA damage upon gene depletion, supporting an R-loop- dependent formation of DNA breaks (Fig. 3a and Supplementary Fig. 3c, d). Next, we assayed whether the levels of R-loops were increased, as determined by IF using the anti-DNA-RNA S9.6 antibody in siINO80 and siSMARCA5 cells that showed the highest γH2AX foci accumulation. Consistently, the mean nuclear S9.6 IF signal was significantly higher in siSMARCA5 (26.3) and siINO80 (34.5) cells compared to siC control cells (19.6), with a nucleolar increase of S9.6 signal also detectable, especially in siINO80 cells (Fig. 3b and Supplementary Fig. 3e, f). Importantly, the S9.6 signal was significantly reduced by overexpressing RNH1, (Fig. 3b and Supplementary Fig. 3e, f), confirming the specificity of detection of DNA–RNA hybrids. Concomitantly, R-loops were also assayed by DNA–RNA immunoprecipitation (DRIP) followed by quantitative PCR (qPCR). In accordance with the results coming from IF experiments, the precipitated material measured as % of input was significantly increased in RNAPII- transcribed genes (RPL13A,TAF9B and FOXP4 in siSMARCA5 cells and RPL13A and FOXP4 in siINO80 cells) when compared to siC cells (Fig. 3c). Moreover, the quantity of immunoprecipitated material was also significantly higher for rDNA genes (28 S in siSMARCA5 and 5’ region and 28 S in siINO80), coincident with nucleolar IF results (Supplementary Fig. 3f). In all cases, the DNA-RNA hybrid signal was reduced dramatically by treating the samples with ribonuclease H (RNH) prior to immunoprecipitation, confirming again that results referred specifically to DNA-RNA hybrids (Fig. 3c). Altogether, these results show that in addition to SWI/SNF, chromatin remodeler factors such as INO80, SMARCA5 and MTA2 control harmful R-loop accumulation in cells, implying a global role of the chromatin remodeling network in R-loop homeostasis and associated genome instability. Chromatin alteration causes R-loop-associated DNA replication impairment and S-phase DNA damage R-loops toxicity is in most cases linked to their ability to block DNA replication during S phase4. Therefore, we studied whether unscheduled R-loop formation in siSMARCA5 and siINO80 cells resulted in TRCs, stalled forks and S phase-associated DNA damage. We first investigated TRC occurrence in SMARCA5 and INO80 depleted cells by Proximity Ligation Assay (PLA) using antibodies against PCNA and the elongating form of RNA polymerase II (RNAPII) phosphorylated at S2 (RNAPII-S2P), as previously described26.RNAPII- S2P + PCNA PLA foci median values increased significantly by 1.75-fold in both siSMARCA5 and siINO80 cells (Fig. 4a). Consistently, increased fork stalling and DNA damage during S phase were also detected in both conditions, as determined by IF of FANCD2 foci using an anti-FANCD2 antibody, and the analysis of DNA break accumulation along the cell cycle by IF of γH2AX foci, using DAPI staining to infer cell cycle phase by the DNA content, as previously established26. Indeed, a significant 2-fold increase in the % of cells with more than 5 FANCD2 foci and median values of γH2AXfocipercellinSphasewereobservedinSMARCA5and INO80-depleted cells (Fig. 4b and Supplementary Fig. 4a). Notably, all the observed phenotypes were RNH1-sensitive;theincreaseinRNAPII- S2P + PCNA PLAs, FANCD2 and S-phase γH2AX foci were significantly suppressed by RNH1 overexpression. These results support that unscheduled R-loop formation boosts TRCs leading to increased fork stalling and DNA breaks (Fig. 4a, b and Supplementary Fig. 4a). To further investigate the function of these chromatin factors at R-loops and TRCs, we analyzed the genome-wide distribution of SMARCA5 and YY1 subunit of INO80 and observed a clear overlap with R-loops and FANCD2 enriched regions, previously determined by DRIPc-seq and ChIP-seq in K562 cells (Fig. 4c, d). Indeed, 3339, 1950 and 816 R-loop peaks colocalized with SMARCA4, SMARCA5 and YY1, respectively, being FANCD2 clearly enriched over these sites, as shown by metaplot analysis of FANCD2 ChIP-seq (Fig. 4e, f). At the gene level, 68.0% target genes of SMARCA5 and 64.0% target genes of YY1 are R- loop-positive genes. Of them, 39.4% and 38.3% are also enriched in FANCD2 (Fig. 4g, h). Remarkably, 81.3% FANCD2-target genes are enriched in SMARCA5 and INO80 proteins, as deduced from ChIP-seq (Supplementary Fig. 4b). Strikingly, despite their general association with R-loops and TRCs, chromatin remodelers mostly associate with specificand mutually exclusive R-loop subsets. Indeed, only 247 R-loops showed enrichment in all three remodeling activities, while 1927, 660 and 299 R-loops colocalized only with SMARCA4, SMARCA5 and YY1, respectively (Fig. 4e). Similarly, metapeak analysis unveiled clearly different profiles of each chromatin remodeler over the R-loop regions. While SMARCA4 shows a clear asymmetric distribution, consistent with previous data, SMARCA5 and YY1 accumulate with a sharp profile at the R-loop sites (Supplementary Fig. 4c). In addition, YY1 is enriched asymmetrically with a gradual increase along the transcription unit until reaching the R-loop prone site. Peak annotations also revealed that YY1-enriched R-loops are at promoter and 5’UTR regions, while R-loops enriched in SMARCA4 and SMARCA5 are enhanced at intronic and intergenic areas (Supplementary Fig. 4d). The R-loop-prone genes enriched in chromatin remodelers do not overlap and display specific features. Among these genes, 248 are bound by all three remodelers, whereas 1096, 371 and 215 only present peaks of SMARCA4, SMARCA5 or YY1, respectively (Supplementary Fig. 4e). Such genes present high expression values when compared to the whole genome, consistent with a higher chance of R-loop Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 5 formation. Nevertheless, the properties of each gene subset are distinct. Thus, R-loop-prone genes enriched in SMARCA4 are significantly longer, SMARCA5-enriched R-loop genes display higher GC content, and YY1-enriched genes show higher expression levels (Supplementary Fig. 4f–h). These results support the view that chromatin remodelers prevent TRCs from resulting in DNA damage during S phase, as reported for the BAF subcomplex of SWI/SNF26,30, under clearly distinct genome contexts. R-loop-prone sites correlate with somatic mutation hotspots in cancer genomes Cancer is a multifactorial disorder where genomic alterations are recurrent, chromatin malfunctioning having a direct impact on the evolution and prognosis of the disease as evidenced by the high number of alterations found at chromatin-modulating genes in malignant cells45. The variety of proteins enriched at genomic regions with high potential to accumulate TRCs, known to be stimulated by R- loops, opens the possibility that such regions are potential hotspots of Fig. 3 | Assessment of DNA damage in siSMARCA5, siINO80, and siMTA2 cells. aPercentage of cells with >5 γH2AX foci in control (siC) and SMARCA5, INO80, MTA2-depleted cells that overexpress (+) or not (−) RNH1. Data expressed as relative to siC. Mean + SEM are plotted (n= 4 (siSMARCA5 and siINO80) and n = 6 (siC and siMTA2) independent experiments). (Unpaired Student’sttest, one-tailed). bQuantification of nuclear S9.6 mean signal intensity in siSMARCA5 and siINO80 cells treated as in (a). Data presented as scatter plot (n> 100 cells examined over 3 independent experiments). Median values are indicated. (Mann–Whitney U-test, two-tailed). cDRIP-qPCR analysis of RNAPII (RPL13A,FOXP4,TAF9B)andRNAPI- transcribed genes (5’and 28S rDNA) of siC, siSMARCA5 (left) and siINO80 (right) cells. Signal values normalized with respect to the siC control and plotted as mean ± SEM (n=5(RPL13A,FOXP4 and 5’rDNA)andn=4(TAF9B and 28S rDNA) independent experiments). (Unpaired Student’sttest, one-tailed). Representative images, with nuclear perimeter highlighted (yellow dashed line), are shown. Gene regions amplified by qPCR in DRIP-qPCR experiments are indicated with a red line on drawings of the genes tested. Scale bars and pvalues are indicated. Source data are provided as a Source data file. See also Supplementary Fig. 3. Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 6 mutations which could be specifically manifested in cancer cells, consistent with the specific mutation signatures and their potential instigators identified in transformed cells46. Therefore, an analysis of themutationpatternsofR-loop-pronesitesmightprovideinsightinto the possible mechanisms of R-loop-induced genetic instability, where chromatin dynamics might play relevant roles. We crossed the mutation data from the COSMIC24,adatabase holding extensive curated somatic mutation information related to human cancers, with genome-wide R-loops obtained by DRIPc-seq26. We first gathered coding and non-coding mutation data from the COSMIC database and divided mutations into SNVs, deletions and insertions and analyzed their abundance at R-loop-hot regions versus a RNH1 Ce lls with > 5 F ANCD2 foci (relative) FF - + - + - + p=0.0004 p=0.0016 p=0.0080 p=0.0040 si C siS MAR CA5 siIN O80 0 1 2 3 RNH1 siC siSMARCA5 FANCD2 siINO80 - + 4 p<0.0001 p<0.0001 p=0.0008 p<0.0001 p<0.0001 0 1 0 20 30 RN AP II -S 2P +P CN A P L A A f oc i/ nu cl eu s A siC 4 7 1 7 2 s iSMAR C C A 5 s iIN O O8 0 - + - + - + RNH1 RNH1 siC siSMARCA5 RNAPII-S2P + PCNA siINO80 - + b b DRIPc SMARCA4 SMARCA5 YY1 FANCD2 RFD genes - 1 c SMARCA4 R-loops SMARCA5 R-loops YY1 (INO80) R-loops 1927 660 299 196 247 74 969 fgh DRIPc SMARCA4 SMARCA5 YY1 FANCD2 RFD genes DRIPc SMARCA4 SMARCA5 YY1 FANCD2 RFD genes N ALF 1 LRR C 8 B L RR C 8C de FANCD2 Mean Coverage R-loop SMARCA4 R-loops 0 . 9 - 1 . 2 2 5 -5kb -5kb SMARCA5 R-loops 0 . 8 5 - 1 . 1 2 5 -5kb -5kb YY1 (INO80) R-loops 0 . 8 - 1 . 1 5 -5kb -5kb R-loop R-loop trx trx 514 285 1166 177 3791 4388 6107 689 937 406 3562 5485 1311 3314 SMARCA5 R-loops FANCD2 YY1 (INO80) R-loops FANCD2 trx c h r 1 : 4 9 , 0 0 2 , 9 2 5 - 8 0 , 3 2 2 , 2 1 7 5 M b 0 - 5 0 - 1 0 0 - 4 0 - 4 1 - 2 . 5 + 1 - 1 0 - 1 0 0 - 1 0 0 - 4 0 - 4 1 - 3 . 5 + 1 1 1 1 1 1 - 1 0 - 1 0 0 - 1 5 0 - 5 0 - 5 1 - 5 + 1 2 5 0 k b 5 0 0 k b 5.5 m 5.5 m Peak correlation Gene correlation Gene correlation Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 7 other regions (Fig. 5a–c). Metadata analysis resulted in noticeable SNV, deletion and insertion metapeaks at R-loop sites (Fig. 5d), suggesting that R-loops are an important source of mutagenesis in cancer cells. Indeed, the mutation burden was significantly higher over R-loop areas compared to the mutation load of the same genomic ranges being randomized along the same genes (Supplementary Fig. 5a). Consistently, mutations were clearly detectable at R-loop-prone genes, such as those routinely tested by DRIP-qPCR (RPL13A,EGR1)(Fig.5b, c). Furthermore, the mutation profiles showed clear enrichments over 5’ and 3’end of R-loop-prone genes (Fig. 5e), as also reported for R-loops. Given the tight connection observed between R-loops and mutagenesis, we also measured the impact of R-loop formation on mutagenesis genome-wide to gain further insight into the mutagenic consequences of these structures in cancer. We first determined mutation hotspots genome-wide by using MACS2 algorithm47,acomputational method commonly used to detect genome regions with coverage enrichments (Supplementary Fig. 5b). Then, we crossed mutation-enriched sites with genes and classified these according to expression and R-loop accumulation in K562 cells. Strikingly, 53.7% of mutation hotspots colocalize with R-loop-prone genes, while only 25.7% and 12.6% coincide with R-loop-reluctant and silenced genes, respectively (Fig. 5f). Further statistical analysis revealed a significant association between R-loop-prone genes and mutation hotspots (Supplemental Fig. 5c). We next compared the impact of mutagenesis in head-on and codirectional R-loop-associated TRCs. Either SNVs, deletions or insertions were found significantly higher at head-on TRCs, consistent with these being the most important compromising genome integrity26,28,41,48–50. Indeed, the number of R-loop-prone sites with more than 1 SNV/kb significantlyrises from 76.9% in co-directional to 81.1% in head-on TRCs, while R-loops with >0.1 indel/kb increases from 20.9–23,6% in co-directional to 27.4–28.2% in head-on TRCs (Fig. 5g). Consistently, SNVs and indels were much more abundant at sites accumulating unscheduled R-loops after depletion of SMARCA4 than of the UAP56/DDX39B DNA-RNA helicase (Supplementary Fig. 5d), in accordance with SMARCA4 but not UAP56 playing a key role helping resolve TRCs26. These results are consistent with a major genotoxicity of head-on TRCs, as widely reported. Finally, we analyzed the impact of chromatin remodeling deficiencies on R-loop mutation burden. We selected data from tumors with at least 1000 mutations (9,046 tumors with 23,013,925 SNVs and 1,636,758 indels) and sorted COSMIC tumor samples based on the status of chromatin remodeler genes; then we compared the mutation load of samples carrying disruptive chromatin remodeler mutations to those with proficient chromatin remodeler genes. Strikingly, R-loop mutation burden was significantly higher in tumors holding SMARCA4,SMARCA5 and INO80 gene deficiencies (Fig. 5h), consistent with a direct role of these activities preventing R-loop mutagenesis in cancer. Altogether, these results are consistent with a major role of chromatin remodelers helping prevent mutation burden at R-loop- prone regions in cancer cells. Chromatin patterns associated with R-loop-prone mutational signatures Mutational processes might be driven by different etiologies during cancer development, each one causing specificandidentifiable mutational patterns commonly referred as mutational signatures. Therefore, we next deepened into the potential causes of R-loop mutagenesis boosts observed in cancer cells and their association to the epigenetic context. First, we analyzed the single-nucleotide variant (SNV) and insertion–deletion (indel) mutational signatures associated to these sites. We used SigProlifer51 to identify the mutational patterns and its similarity to previously described signatures in the mutation dataset. The result revealed 17 SNV and 13 indel mutational signatures (Supplementary Figs. 6 and 7 and Supplementary Tables 5 and 6) matching a wide range of mutation etiologies. Then, wecomputed the numberof mutations corresponding to each signature present at R-loop sites and performed metaplot analysis of them around +/−50 kb from the R-loop peak centers. Notably, mutation signatures peaked at R-loop sites, being those associated with dysfunction of DNA repair processes the most abundant. Indeed, mutations resulting from defects in homologous recombination (HR; BRCA1 and BRCA2 mutations), transcription-coupled Nucleotide Excision Repair (NER-TCR), DNA mismatch repair (MMR), translesion synthesis (TLS; polymerase η), base excision repair (BER), topoisomerases (topoisomerase 1 and 2A), non-homologous end joining (NHEJ), replication slippage and AID/ APOBED-mediated cytidine deamination scored high at these sites (Fig. 6a and Supplementary Fig. 8a, b). Consistently, deregulation of several factors involved in these processes, including BRCA1, BRCA2 and TOP1 induces R-loop-dependent genome instability9,19,20,52. Remarkably, mutations associated with HR-defects and replication slippage were among the most abundant (Fig. 6a and Supplementary Fig. 8a, b), suggesting that fork stalling at TRCs may be a major source of mutagenesis at these regions. Importantly, mutational signature comparisons unveiled a major impact of HR dysfunction, NER-TCR and the AID/APOBEC family of cytidine deaminases on head-on TRCs in cancer cells (Fig. 6aand SupplementaryFigs.9and10).Similarly, signatures analogous to those induced by reactive oxygen species (ROS) or aristolochic acid exposure, as well as one with unknown etiology characterized by increased >5 bp insertion mutations were also higher at head-on TRCs. In addition, despite not reaching significant differences, a higher abundance of mutational footprints related to translesion synthesis (TLS) polymerase ηor topoisomerase was also appreciable at head-on TRCs. We next investigated the association of chromatin factors enriched at TRCs in our screening with the different mutational patterns peaking at R-loop-prone regions. To this aim, we first measured the frequency of mutational signatures at R-loop-prone sequences enriched for each of the most relevant TRC-enhanced epigenetic factors identified (Fig. 1d) and compared it to the frequency at all R-loop- prone genomic regions. The results unveiled differential associations that were grouped into 7 clusters (Supplementary Fig. 11). Consistent with previous data26, SWI/SNF complex factors SMARCA4 and ARID1B Fig. 4 | Analysis of DNA damage along the cell cycle in siINO80 and siSMARCA5 cells. a Quantification of RNAPII-S2P + PCNA PLA (TRCs) foci in control (siC), siSMARCA5 and siINO80 cells. Data presented as scatter plot (n>50cells examined over 3 independent experiments). Median values are indicated. (Mann- Whitney U test, two-tailed). bPercentage of cells with more than 5 FANCD2 foci in siC, siSMARCA5 and siINO80 cells that overexpress (+) or not (−) RNH1. Data expressed as relative to siC. Mean + SEM are plotted (n=5(RNH1−)andn=6 (RNH1+) independent experiments). (Unpaired Student’sttest, one-tailed). cRepresentative screenshot of a genome region showing co-localization of R-loops (DRIPc; light blue), SMARCA5 (yellow), YY1 (INO80; green), SMARCA4 (red) and FANCD2 (purple) at sites with high R-loop abundance. Replication fork directionality (RFD) is also shown. dZoomed-in examples of R-loops (DRIPc; light blue), SMARCA5 (yellow), YY1 (INO80; green), SMARCA4 (red) and FANCD2 (purple) colocalization. eVenn diagram showing correlation between SMARCA4 (red), SMARCA5 (yellow) and YY1 (INO80; green) peaks colocalizing with R-loops. fMetapeak analysis of FANCD2 at R-loops colocalizing with SMARCA4, SMARCA5 or YY1. Transcription direction is indicated. gVenn diagram showing R-loop (DRIPc; light blue), FANCD2 (purple) and SMARCA5 (yellow) genome-wide co-occurrence in control K562 cells. hVenn diagram showing R-loop (DRIPc; light blue), FANCD2 (purple) and YY1 (INO80; green) genome-wide co-occurrence in control K562 cells. Representative images, with nuclear perimeter highlighted (yellow dashed line), are shown. Scale bars and pvalues are indicated. RFD replication fork directionality, trx transcription. Arrows under trx indicate transcription direction. Source data are provided as a Source data file. See also Supplementary Fig. 4. Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 8 were clustered together and enriched at sequences with high SNV signatures driven by APOBEC mutagenesis and HR deficiencies and sequences with indel patterns observed at sites of head-on TRCs. In contrast, SMARCA5 and MTA1/2/3 strongly associated with APOBEC mutagenesis, MMR deficiencies and topoisomerase 1/2A-dependent mutagenesis, as also observed for H4K20me1, a histone modification tightly associated with DNA repair and replication53.Otherfactors including YY1, an INO80 subunit, or SIN3B HDACs correlated with APOBEC-induced mutagenesis, MMR deficiencies, DNA replication slippage and different patterns of unknown etiology. Strikingly, DRIPc SNVs deletions insertions genes EGR1 0-230 0-0.26 0-0.15 0-0.02 DRIPc SNVs deletions insertions genes -5.0kb +5.0kb 0.8 1.0 1.2 SNVs trx 1.5 2.0 2.5 3.0 3.5 deletions trx 2.5 3.0 3.5 4.0 insertions trx RPL13A RPS11 DRIPc SNVs deletions insertions genes 0-75 0-0.15 0-0.1 0-0.03 a bcd 50kb Mean Cov. (x10-2) h Mean Cov. (x10-3) Mean Cov. (x10-4) SMARCA4 proficient R-loop Mutation Burden (R-loop mutations / tumor sample) SMARCA4 deficient 3000 2000 1000 0 p < 2.22e−16 p < 2.22e−16 SMARCA5 proficient R-loop Mutation Burden (R-loop mutations / tumor sample) SMARCA5 deficient 3000 2000 1000 0 p < 2.22e−16 INO80 proficient R-loop Mutation Burden (R-loop mutations / tumor sample) INO80 deficient 3000 2000 1000 0 p<0.0001p<0.0001 p<0.0001 e 0.2kb -0.2kb gene body 0.5 1 3 2.5 2 1.5 TSS TES Mutations (Mean Coverage) (x10-2) R-loop-prone genes f p=0.0429 p=0.0411 p=0.0055 SNVs (>1SNV/kb) deletions (>0.1del/kb) insertions (>0.1ins/kb) 280/364 (76.9%) 1442/1779 (81.1%) 86/364 (23.6%) 502/1779 (28.2%) 76/364 (20.9%) 487/1779 (27.4%) CD HO CD HO CD HO g 60 40 20 0 Mutation hotspots (%) Expression: R-loops: - -++ -+ 12.6 25.7 53.7 p<0.0001 p<0.0001 p<0.0001 R-loop -5.0kb +5.0kb R-loop -5.0kb +5.0kb R-loop ZNF707 0-75 0-0.15 MAPK15 NR_136279.1 SCRIB PUF60 EPKK1 PLEC PARP10 SPATC1 SMPD5 EXOSC4 CYC1 WDR97 HGH1 0-0.1 0-0.03 chr8: 143,676,271-144,143,192 Fig. 5 | Analysis of the mutational landscape at TRCs in cancer. a Representative screenshot of a genome region showing co-localization of R-loops (DRIPc; blue) with SNVs (orange), deletions (green) and insertions (yellow). bExample of an R- loop-prone gene (EGR1) showing co-localization with SNVs (orange), deletions (green) and insertions (yellow). cExample of an R-loop-prone gene (RPL13A) showing co-localization with SNVs (orange), deletions (green) and insertions (yellow). dMutation metanalysis. Mean coverage of SNVs (orange), deletions (green) and insertions (yellow) around +/−5 kb of R-loop peaks. Transcription direction is indicated. eMetagene analysis of mutations over R-loop-prone genes. The arrow below the graph indicates transcription direction. fPercentage of mutation hotspots in cancer colocalizing with R-loop-prone genes, expressed genes reluctant to R-loop formation and silenced genes. (Chi-square with Yates’correction test, twotailed). gQuantification of the number of sites with >1SNV/kb, >0.1 deletions/kb and >0.1 insertions/kb at head-on and co-directional TRCs. (Fisher’s exact test, onetailed). hComparison of R-loop mutation burden between tumor samples proficient and deficient in SMARCA4 (left panel), SMARCA5 (middle panel) and INO80 (right panel). Data presented as box plot. The lower and upper hinges of the box plots correspond to the first and third quartiles (the 25th and 75th percentiles), while the upper and lower whiskers extend from the hinge to the largest and smaller value no further than 1.5× inter-quartile range, respectively. Data beyond the end of the whiskers are considered outliers and plotted individually. n=9046 biologically independent tumors sorted according to gene status. (Wilcoxon test, two-tailed). Scales and pvalues are indicated. SNV single-nucleotide variant, HO head-on, CD co-directional, trx transcription. Arrows under trx indicate transcription direction. Source data are provided as a Source data file. See also Supplementary Fig. 5. Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 9 74. Rice, P., Longden, L. & Bleasby, A. EMBOSS: The European Molecular Biology Open Software Suite. Trends Genet. 16, 276–277 (2000). 75. Ramírez, F. et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 44, W160–W165 (2016). 76. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9,357–359 (2012). 77. Yu, G., Wang, L. G. & He, Q. Y. ChIP seeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 31,2382–2383 (2015). 78. Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38,576–589 (2010). 79. Kuhn, R. M., Haussler, D. & James Kent, W. The UCSC genome browser and associated tools. Brief. Bioinformatics 14, 144–161 (2013). 80. Perez-Llamas, C. & Lopez-Bigas, N. Gitools: analysis and visualisation of genomic data using interactive heat-maps. PLoS ONE 6, 19541 (2011). 81. Neph, S. et al. BEDOPS: high-performance genomic feature operations. Bioinformatics 28,1919–1920 (2012). 82. Quinlan,A.R.&Hall,I.M.BEDTools:aflexible suite of utilities for comparing genomic features. Bioinformatics 26,841–842 (2010). 83. Islam, S. et al. Uncovering novel mutational signatures by de novo extraction with SigProfilerExtractor. Cell Genomics https://doi.org/ 10.1016/j.xgen.2022.100179 (2022). 84. Aragon, T. J. epitools: Epidemiology Tools R package version 0.5-7. https://CRAN.R-project.org/package=epitools (2012). 85. Galili, T. dendextend: An R package for visualizing, adjusting and comparing trees of hierarchical clustering. Bioinformatics 31, 3718–3720 (2015). 86. RStudio Team. RStudio: integrated development environment for R (RStudio Team, 2019). 87. Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 29,24–26 (2011). 88. Bayona-Feliu, A. et al. The chromatin network helps prevent cancerassociated mutagenesis at transcription-replication conflicts. github/zenodo https://doi.org/10.5281/zenodo.8359279 (2023). 89. Bayona-Feliu, A. et al. The chromatin network helps prevent cancerassociated mutagenesis at transcription-replication conflicts. github/zenodo https://doi.org/10.5281/zenodo.8390817 (2023). Acknowledgements We thank F. Supek (Institute for Research in Biomedicine, Barcelona) for computational resources provided and discussions. Research was part of the project I + D + i PID2019-104270GB-I00/BMC, funded by MCIN/ AEI/10.13039/501100011033/, as a main funding source, as well as grants from the European Research Council (ERC2014 AdG669898 TARLOOP) and Foundation “Vencer el Cancer”. A.B.-F. was supported in part by Juan de la Cierva postdoctoral contracts (FJCI-2017-34536/ IJC2020-044963-I) from the Spanish Ministry of Science and Innovation, E.H.-M. by a Marie Sklodowska-Curie Action postdoctoral contract (H2020-MSCA-IF-2020-RteRloop –101027467) and M.E.S.-O. by a predoctoral FPI contract of the Ministry of Science and Innovation. Author contributions A.B.-F. and A.A. designed the study and the experiments. A.B.-F. performed most of the experiments and the bioinformatic analysis. E.H.-M. and N.B.-F. contributed with specific experiments. I.G.-F. and M.E.S.-O. performed part of the bioinformatic analysis. A.B.-F. and A.A. wrote the manuscript. All authors read, discussed, and agreed with the final version of this manuscript. Competing interests The authors declare no competing interests. Additional information Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-42653-0. Correspondence and requests for materials should be addressed to Aleix Bayona-Feliu or Andrés. Aguilera. Peer review information Nature Communications thanks Zhijun Duan and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available. Reprints and permissions information is available at http://www.nature.com/reprints Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2023 Article https://doi.org/10.1038/s41467-023-42653-0 Nature Communications | (2023) 14:6890 16