Repeat RNAs associate with replication forks and post-replicative DNA

Repeat RNAs associate with replication forks and post-replicative DNA HELENE M. GYLLING, 1 CRISTINA GONZALEZ-AGUILERA, 1,5 , 6 MARTIN A. SMITH, 2,3,7 , 8 DOMINIK C. KACZOROWSKI, 2 ANJA GROTH, 1,4 and ANDERS H. LUND 1 1 Biotech Research and Innovation Centre, University of Copenhagen, Copenhagen, 2200, Denmark 2 Garvan Institute of Medical Research, Darlinghurst, New South Wales 2010, Australia 3 St-Vincent’s Clinical School, Faculty of Medicine, UNSW Sydney, Darlinghurst, New South Wales 2010, Australia 4 The Novo Nordisk Center for Protein Research (CPR), University of Copenhagen, Faculty of Health Sciences, University of Copenhagen, 2200 Copenhagen, Denmark ABSTRACT Noncoding RNA has a proven ability to direct and regulate chromatin modifications by acting as scaffolds between DNA and histone-modifying complexes. However, it is unknown if ncRNA plays any role in DNA replication and epigenome maintenance, including histone eviction and reinstallment of histone modifications after genome duplication. Isolation of nascent chromatin has identified a large number of RNA-binding proteins in addition to unknown components of the replication and epigenetic maintenance machinery. Here, we isolated and characterized long and short RNAs associated with nascent chromatin at active replication forks and track RNA composition during chromatin maturation across the cell cycle. Shortly after fork passage, GA-rich-, alpha- and TElomeric Repeat-containing RNAs (TERRA) are associated with replicated DNA. These repeat containing RNAs arise from loci undergoing replication, suggesting an interaction in cis. Postreplication during chromatin maturation, and even after mitosis in G1, the repeats remain enriched on DNA. This suggests that specific types of repeat RNAs are transcribed shortly after DNA replication and stably associate with their loci of origin throughout the cell cycle. The presented method and data enable studies of RNA interactions with replication forks and post-replicative chromatin and provide insights into how repeat RNAs and their engagement with chromatin are regulated with respect to DNA replication and across the cell cycle. Keywords: RNA; DNA replication; chromatin; repeats INTRODUCTION Mammalian genome duplication relies on correct DNA replication and reestablishment of the local chromatin environment on the two newly synthesized sister chromatids. DNA replication is initiated by sequential activation of several origins of replication in S-phase (Bell and Dutta 2002; Mechali 2010). The replication program is carefully orchestrated according to genome organization, with a trend for transcriptionally active euchromatin to replicate early in S-phase (Marchal et al. 2019). The replication machinery disrupts chromatin ahead of the replication fork and must correctly restore the epigenetic landscape on both new daughter strands (Alabert and Groth 2012). Parental histones with their large variety of post-translational modifications (PTMs) are segregated onto the daughter strands and mixed with nucleosomes assembled from new histones (Alabert et al. 2015; Annunziato 2015). Parental histone recycling is remarkably accurate, allowing the histone PTM landscape to be reproduced after replication but with PTM levels twofold reduced due to dilution with new histones (Alabert et al. 2015). Post-replication, in a process termed chromatin restoration, new histones acquire modifications identical to those of the nearby parental histones in order to maintain epigenetic regulation and preserve 5 Present address: Departamento de Biología Celular, Universidad de Sevilla, 41012, Seville, Spain 6 Present address: Centro Andaluz de Biología Molecular y Medicina regenerativa (CABIMER), Universidad de Sevilla-CSIC- Universidad Pablo de Olavide, 41092, Seville, Spain 7 Present address: CHU Sainte-Justine Research Centre, Montreal, H3T 1C5, Canada 8 Present address: Department of Biochemistry and Molecular Medicine, Faculty of Medicine, University of Montreal, H3C 3T5, Canada Corresponding author: [email protected] Article is online at http://www.rnajournal.org/cgi/doi/10.1261/rna. 074757.120. © 2020 Gylling et al. This article is distributed exclusively by the RNA Society for the first 12 months after the full-issue publication date (see http://rnajournal.cshlp.org/site/misc/terms.xhtml). After 12 months, it is available under a Creative Commons License (Attribution- NonCommercial 4.0 International), as described at http://creativecommons.org/licenses/by-nc/4.0/. 1104 RNA (2020) 26:1104–1117; Published by Cold Spring Harbor Laboratory Press for the RNA Society cell identity. This is a highly heterogeneous process that can take minutes or hours depending on the nature and location of histone marks (Alabert et al. 2015; Reveron- Gomez et al. 2018), the regulation and mechanism of which is largely unknown. Generally, PTMs associated with active transcription are restored rapidly, with recent evidence in mouse embryonic stem cells (mESC) showing that RNA Polymerase II transcription recommences within 30 min after DNA replication (Stewart-Morgan et al. 2019). The epigenetic landscape is largely shaped by histone modifications and DNA methylation, which are deposited and removed by histone and DNA modifiers (Rothbart and Strahl 2014). Recruitment of these factors to precise genomic sites is facilitated by transcription and recognition of chromatin features, including long noncoding RNAs (lncRNA) (Kouzarides 2007; Johnson and Straight 2017). Over the last decade, lncRNAs have broadened our understanding of the roles of RNA from encoding proteins to having important functions in almost all cellular pathways via a number of different functionalities, including binding to DNA, other RNAs and proteins (Quinn and Chang 2016). With the exception of Y-RNAs, which regulate replication initiation (Christov et al. 2006; Ge and Lin 2014), the role of RNA at ongoing replication forks and in chromatin maturation remains elusive. Nascent chromatin capture (NCC) has been used to identify proteins associated with active replication forks and newly replicated DNA (Alabert et al. 2014). Surprisingly, a substantial number of these proteins have a described function in RNA binding and processing. Based on this and the many functional roles of RNA, such as the scaffolding of macromolecular complexes, we speculated that RNA could play a role in chromatin replication. To address this question, we developed an NCC- based method to detect RNA at active replication forks and maturing chromatin. The method utilizes pulse-label- ing of DNA in replicating synchronized cells with biotin dUTP, enabling isolation of chromatin at specific time points post replication. We find a number of RNA transcripts arising from genes to be moderately enriched specifically at nascent chromatin. These genes are not replicated at the time of labeling, suggesting that the RNA associate with nascent chromatin in trans. In addition, we find considerable enrichment of specific groups of repeat containing RNAs at mid/late S-phase replicated DNA shortly after replication. These RNAs arise from DNA loci that are being replicated at the time of labeling, implying that these RNAs associate with nascent chromatin in cis. Notably, the repeat RNA species remain associated with the replicated chromatin loci for at least 10 h, suggesting that these RNAs are stably bound to chromatin. We propose that repeat RNAs are transcribed shortly after fork passage as a consequence of increased accessibility of DNA upon incorporation of new acetylated histones (Annunziato 2012) and that they remain associated with and influence chromatin structure at the loci across the cell cycle. RESULTS Development of NCC-RNA-seq method With the aim to identify specific RNAs associated with replication forks and maturing chromatin, we optimized the NCC technique for RNA isolation (Fig. 1A). This protocol allows for the isolation of active replication forks and nascent chromatin by incorporation of biotin-dUTP during DNA synthesis. Maturation of the pulse-labeled newly replicated chromatin can then be followed across the cell cycle (Alabert et al. 2014, 2015). To allow comparison of RNA- seq data with published NCC proteomics and DNA-seq data from HeLa cells (Alabert et al. 2014; Reveron-Gomez et al. 2018), we adapted the same synchronization strategy—that is, pulse labeling cells in mid-S phase 3 h after release from a G1/S block (Fig. 1B). We isolated RNA associated with nascent chromatin (N), and mature chromatin harvested 2 (M2), and 10 h (M10) after DNA replication (Fig. 1B,C; Alabert et al. 2015). We incorporated several controls, including a negative control without b-dUTP labeling (no b-dUTP) to account for binding of nonspecific RNAs and chromatin immunoprecipitation (ChIP) of histone H3, representing total chromatin-bound RNA. Input controls were included to account for bias between highly and lowly expressed RNAs. Importantly, we pull down PCNA, a key DNA polymerase processivity factor, in the nascent sample only, ensuring isolation of active replication forks (Supplemental Fig. S1B). We also observed an increase in the level of H3K9me3 with increasing chase time (Supplemental Fig. S1B), consistent with the restoration marks on new histones after chromatin replication. Because isolated RNAs from the NCC samples were 20–4000 nt in length (Supplemental Fig. S1A), we performed sequencing of both long and short RNAs to obtain a comprehensive overview of RNAs associated with the replication complex and maturing chromatin. Long RNA- seq experiments were performed in quadruplicates with input and H3 controls taken at the time of labeling (Fig. 1C). Initial RNA-seq analysis using genome-guided reference transcriptome quantification revealed that more than 50% of mapped reads were lost due to reads mapping outside of known transcripts. However, mapping to the genome resulted in few unmapped reads. Instead of relying on annotated genes and transcripts, we performed ab initio transcriptome assembly using Trinity, which assembles transcripts based on the raw sequencing reads without a reference sequence (Grabherr et al. 2011). In addition, we mapped reads to the human genome version 38 (hg38) and used another annotation-free algorithm— DERfinder—to identify expressed regions based on clustering of reads in the genome, for example, unannotated RNAs at replication forks www.rnajournal.org 1105 exons (Frazee et al. 2014; Collado-Torres et al. 2017a,b). Both methods were analyzed with EdgeR for enrichment analysis (Jurka 2000; Jurka et al. 2005). Based on initial principal component analysis (PCA) and clustering of samples after EdgeR processing, one nascent replicate (N_112) was removed from the analysis since it clustered with the input samples instead of the other nascent samples (Supplemental Fig. S1C). We prepared the short RNA-seq in triplicates using a capture and amplification by tailing and switching (CATS) strategy, which generates cDNA libraries ranging from 3 to 300 nt in length. In these experiments, we included cell cycle-matched input controls for all time points, H3 IP and a negative control without b-dUTP incorporation (Fig. 1C). No rRNA depletion or size selection was performed to reduce any potential enrichment bias (excluding the removal of primers during cDNA library preparation). This resulted in libraries that, in some instances, had more than 60% of reads mapping to repeat-con- taining RNAs. To mitigate the impact of multimapping repetitive reads and allow for meaningful quantification and statistical analysis of the data, we first mapped reads to human repeat sequences to uniquely map repeat-containing reads, and subsequently aligned the unmapped reads to the hg38 genome (Fig. 1D). Expressed regions were identified and quantified using DERfinder and the analysis of the repeat mapping and hg38 mapping was done with EdgeR (Supplemental Fig. S1D,E). Long RNA sequencing identifies genes that are lowly yet specifically enriched at nascent chromatin, whereas GA-rich and telomeric RNAs are highly enriched at nascent, 2 and 10 h mature chromatin We first investigated the diversity of RNAs associated with replication forks. Interestingly, we found the nascent samples to cluster separately in the PCA plot, and when comparing data from the nascent samples (N) against all other samples (M2, M10, H3, and input), we found 737 RNAs specifically associated with the replication fork (nascent hits) (Fig. 2A; Supplemental Fig. S1C;Supplemental Table S1). Overlapping the RNAs with annotation files for repeats and genes showed that all nascent hits (except 15 unannotated regions) mapped to the exons of 309 different protein-coding genes, spread across all chromosomes. These RNAs displayed low enrichment between the nascent and the M2 samples with a log fold-change (logFC) < 2.5 (except one transcript with logFC 3.8). Based on this low enrichment, it is unlikely that the observed RNAs associate with all actively replicating forks, as this would likely have resulted in a much higher logFC between the nascent and mature chromatin. Similarly, we analyzed specific enrichment for transcripts in the M2 and M10 samples, respectively, but found only a subtle enrichment in terms of logFC. All three gene lists (N, M2, and M10 hits, Supplemental Table S1) were analyzed for characteristics via B A C D FIGURE 1. Developing nascent chromatin capture RNA-seq (NCC RNA-seq). (A) Workflow of NCC RNA-seq protocol. (B) Experimental setup illustrating synchronization, labeling and sample harvesting. FACS diagrams showing distribution of cells in cell cycle phases at harvest time points (also indicated in blue cell cycle illustration), N = nascent chromatin, M2 = 2 h mature chromatin, M10 = 10 h mature chromatin. (C) Schematic representation of sample setup for long RNA and short RNA sequencing. Orange stars indicate samples that were included in short but not long RNA-seq. Three out of four samples were labeled 3 h into S phase with a pulse of biotin dUTP. Nascent (N) and unlabeled samples were harvested 15 min after labeling, whereas mature samples were collected 2 h (M2) and 10 h (M10) later. The unlabeled sample in short RNA-seq was split in two and used for histone 3 (H3) immunoprecipitation and as a control for unspecific binding to Streptavidin T1 beads (no b-dUTP). Input controls were taken from indicated samples before cross-linking. (D) Bioinformatics pipeline for long (red) and short (blue) NCC RNA-seq. Gylling et al. 1106 RNA (2020) Vol. 26, No. 9 Enrichr, and although pathways, ontologies, etc. showed statistical enrichment, no consistent pattern of enrichment was found within each data set or in common (Chen et al. 2013; Kuleshov et al. 2016). We also compared nascent, M2 and M10 samples against input and H3, which revealed a significant enrichment for a large number of RNAs (Fig. 2B). Based on this we created a hit list (long hits), composed of transcripts that were commonly enriched between N, M2 and M10 when individually compared to input and H3 (FDR ≤0.05 and logFC > 0) (Supplemental Table S1). Long hits consisted of 405 transcribed sequences annotated as being part E BA F CD G FIGURE 2. GA-rich repeats are enriched in N, M, and M10 samples from long RNA-seq. (A) RNAs in long RNA-seq plotted according to log foldchange (logFC) versus average counts per million (CPM) from nascent (N) and 2 h mature (M2) sample comparison. Red and gray points indicate significantly enriched (FDR = <0.05) RNAs in N sample compared to M2. (B) Enriched RNAs (FDR= <0.05) from N versus input comparison plotted and colored according to RNA type (gene [red], repeat [blue], or not available [NA, gray]). (C) Heat map of long hits (n= 841) according to sample and repeat type. Colors represent relative difference of transcript levels in samples [log(CPM + 0.1)]. (D)Long hits repeats (n= 315) according to repeat type. (E) Relative proportion of repeat types in the human genome, n= 5,520,017 (UCSC repeat masker hg38 [Jurka 2000]). (F) Simple repeats (SR) and low complexity (LC) subtypes, n= 806,538 (14.6% of all repeats) with color coded proportions of repeats found in long hits. (G)Long hits repeats, n= 315 relative proportion and coloring according to F. RNAs at replication forks www.rnajournal.org 1107 of genes, 315 as repeats and 121 as unannotated regions (NA). In general, the RNAs that mapped to genes had a low logFC (typically below 2.5) whereas repeats and NAs had considerably higher enrichments (logFC 2.5–10) in the N, M2, and M10 samples over input (Fig. 2C; Supplemental Fig. S2A). These repeats are normally not detectable, or expressed at very low levels, but were enriched in the H3 control compared to input, indicating that the transcripts associate with chromatin and that their enrichment in the N, M2, and M10 samples does not reflect unspecific binding to the streptavidin beads. Strikingly, almost all the enriched long hit repeat-asso- ciated sequences were dinucleotide repeats (CT- or GA-rich) and distributed on all chromosomes, both in intergenic regions and gene bodies (Fig. 2D; Supplemental Fig. S2B,C). It is noteworthy to mention that G-rich telomere sequences (TTAGGG) were also observed to be significantly enriched (NvsIn FDR = 268 × 10 −05 ). GA/CT-rich repeats account for <3% of the total number of repetitive sequences, composing ∼0.33% of the human genome, making it unlikely that the enrichment for this specific type of repeat RNA is by chance (Fig. 2E–G; Supplemental Fig. S2D,E). In order to establish if the transcripts were CT- or GA-rich, and to validate the enrichment independently of genome mapping (avoiding quantitative artifacts derived from multimappers), we performed transcriptome assembly, quantification and DE analysis. Mapping the enriched sequences from the transcriptome assembly analysis back to the genome provided parallel confirmation that the enriched repeats were exclusively GA- and not CT-rich RNAs. Short RNA sequencing confirms GA enrichment and identifies centromeric and telomeric repeats enriched at nascent, 2 and 10 h mature chromatin To get a comprehensive view of RNAs associated with replicated chromatin and to validate the GA-repeat enrichment observed in the long RNA-seq experiments, we additionally performed short RNA sequencing. Anticipating many repeat-containing RNAs, the resulting reads were analyzed by first mapping them to repeat sequences, then unmapped reads were subsequently mapped to hg38. Similar to the long RNA-seq results, the independently analyzed repeat mapping showed no enriched RNAs at the N, M2, or M10 time points when compared to each other, but did present an enrichment of specific repeats in newly replicated and maturing chromatin samples when compared to controls (Fig. 3A; Supplemental Table S2). These repeats were mainly dinucleotide (GA/CT), alpha centromeric repeats and TERRA (telomere repeat) sequences, validating the GA repeat and TERRA enrichment observed in the long hits. The genome mapping analysis also did not show any specific hits at a single time-point, yet an enrichment was observed for repeats at nascent and maturing chromatin when compared to controls (short hits) (Fig. 3B; Supplemental Table S3). In contrast to the long RNA-seq, most of the enriched repeats mapped to a few chromosomes and were concentrated around heterochromatin regions, telomeres, centromeres and pericentromeres (Fig. 3C; Supplemental Fig. S3A). The enriched repeats are annotated mainly as GA/CT rich, centromeric and pericentromeric sequences, such as BSR/beta satellite repeats and CCATT repeats (centromeric repeats) (Fig. 3D; Catasti et al. 1994). CCATT repeats constitute most of HSATII repeats, which are mainly located on Chromosomes 1, 2, 10, and 16 (Tagarro et al. 1994). Repeat masking of transcripts showed that almost all of the unannotated hits (NAs) on Chr 1 and 16, including one hit on Chr 7 and one on Chr 22, were HSATII repeats (Supplemental Fig. S3B;Supplemental Table S4), whereas sequence analysis of the unmasked NA hits (n= 60) showed no strong enrichment of sequence motifs. Taken together, the NCC short RNA sequencing revealed an enrichment of centromeric, pericentromeric (alpha, CCATT, HSATII, BSR) and telomeric repeats (TERRA) and GA-rich sequences shortly after genome replication (nascent chromatin), 2 h and 10 h later (mature chromatin), compared to controls. Although the short hits did not contain as many GA-rich RNA repeats as the long hits, we saw a tendency for most of the long hits to be enriched in the short RNA-seq N, M2, and M10 samples (Supplemental Fig. S3C). The NCC short RNA-seq thus substantiates the findings from long RNA-seq and identified a number of shorter RNA repeats located at telomeres and centromeres, exposing a similar enrichment pattern to GA-rich RNAs. Repeat RNAs are replicated at the time of labeling To further understand why RNA is associated with replicated chromatin, we wanted to investigate if the enriched transcripts in the NCC RNA-seq data were replicated at the time of labeling, an indicator of association in cis.To do so, we compared these data to our previously published NCC DNA-seq (Reveron-Gomez et al. 2018). This work sequenced NCC isolated DNA using the same HeLa S3 cell line and a similar NCC protocol (synchronization, labeling time, pulldown) as used in our study, allowing us to identify which genomic regions were replicated at the time of biotin labeling. Due to large differences between hg19 and hg38 in the centromeric regions where many of the short hits are located, we chose to reanalyze the data by first mapping to Repbase and subsequently to hg38. Of the 737 transcripts observed to be specifically enriched in the N samples, 127 (17%) overlapped with replicated loci identified by NCC DNA-seq (47 out of 309 genes) (Fig. 4A). These transcripts might thus be enriched due to concomitant transcription and replication of the Gylling et al. 1108 RNA (2020) Vol. 26, No. 9 loci. However, as this also roughly corresponds to the proportion expected based on overlap of the NCC DNA- seq data with all transcribed genes in the long RNA-seq data set (∼13%), the interactions are more likely to occur in trans. On the other hand, long hit repeats and NAs showed a strong overlap with replicated regions compared to the overlap in the entire long NCC RNA-seq transcriptome (Fig. 4B). Most of the long hits were GA-rich. To find a common pattern in these sequences, we measured their C/G and GA/CT content, which showed that the vast majority of sequences had >40% C/G, >80% CT/GA content and that 99% were 199 nt or longer (Supplemental Fig. S4A– C). Based on this, we searched the genome for similar sequence compositions and identified ∼10,000 regions across all chromosomes (GA regions, Supplemental Table S5). Of these, 404 overlap with 435 long hit repeats and NAs and ∼37% are replicated in mid/late S-phase (Supplemental Fig. S4D). This high number of replicated regions, however, does not show evidence of GA regions being specifically replicated at the time of labeling since the replicated loci in the NCC DNA-seq data set covers ∼37% of the genome. Mapping the NCC DNA-seq data to Repbase showed alpha, BSR centromeric and TERRA repeats were significantly enriched, indicating that these repeats were replicated at the time of labeling (Fig. 4C; Supplemental Table S6). This is in accordance with previous studies, which shows that pericentromeric and centromeric regions are replicated in mid/late S-phase (O’Keefe et al. 1992; Erliandri et al. 2014). Similarly, repeats and NA regions enriched in short hits overlap >90% with replicated loci identified in NCC DNA-seq, strongly suggesting that NCC-RNA enriched centromeric, pericentromeric and TERRA repeats, are B A CD FIGURE 3. Centromeric, telomeric, and GA-rich repeat RNAs are enriched in N, M2, and M10 samples from short RNA-seq. (A) Repbase mapped RNAs plotted according to log fold-change (logFC) versus average counts per million (CPM) in nascent (N) versus input (nascent) comparison (logCPM NvsIn). Colored points (according to repeat type) highlight significantly enriched repeat RNAs in N, M2, and M10 samples compared to all controls (inputs [N, M2, and M10], H3, and no b-dUTP). ALR = alpha repeat, CCCTAA = telomere repeat. (B) Heatmap of short hits (n= 392) according to sample and repeat type. Colors represent relative difference of transcript levels in samples [log(cpm + 0.1)]. (C)Short hits distribution on chromosomes and coloring according to RNA type. (D) Types of repeats present in short hit repeat RNAs (n= 135). Red nucleotides in repeats correspond to the same repeat sequence (CCATT)n or the reverse complementary (AATGG)n. RNAs at replication forks www.rnajournal.org 1109 replicated at the time of labeling (Fig. 4D). This suggests that repeat RNAs commonly enriched in N, M2, and M10 samples are associated with biotin-dUTP labeled chromatin in cis. It also suggests that the RNA is either being transcribed during locus replication, or reassociates with its locus after fork passage, and that this association persists at least 10 h after replication. The latter could be driven by the reassociation of RNA containing complexes with chromatin early in chromatin restoration. NCC-enriched repeat RNAs do not form R-loops but coding genes prone to R-loop formation are enriched at replication forks R-loops are DNA:RNA hybrids that typically form when a transcribed sequence invades and anneals to the complementary DNA strand in cis (Aguilera and Garcia-Muse 2012). G-rich RNA is prone to make R-loops, and alpha and TERRA sequences have been shown to facilitate specific functions via R-loop formation (Roy et al. 2008; Reddy et al. 2011; Balk et al. 2013; Arora et al. 2014; Groh et al. 2014; Kabeche et al. 2018; Lee et al. 2018). We thus wanted to investigate if the repeat RNAs enriched at nascent, 2 h and 10 h mature DNA could be associating with their parental locus in cis via R-loop formation. To do this, we reanalyzed published DRIP-seq data sets by applying the above-mentioned iterative repeat and reference RNA-seq mapping protocol (Hamperl et al. 2017). The Repbase analysis confirmed TERRA as R-loop forming, in accordance with previous studies (Fig. 5A; Lee et al. 2018). The alpha repeats did not show any enrichment in the R-loop data, suggesting that alpha repeats are not associating with chromatin via R-loops in S- and G1-phase. Since mitotic cells are a small population in asynchronous BA CD FIGURE 4. NCC DNA-seq data show enriched repeat RNAs are replicated at time of labeling. (A) Proportion of replicated transcripts in entire long NCC RNA-seq transcriptome (all long) and in nascent hits, plotted according to transcript type (Genes = UCSC annotated gene in hg38 genome build or NA= not available). Transcripts are replicated if their mapped genomic locus overlaps with a replicated region in NCC DNA- seq. (B)Long hits compared to entire long NCC RNA-seq transcriptome as in A(repeats = UCSC genome browser RepeatMasker annotations for hg38). (C) Repbase analysis of NCC DNA-seq replicate (rep6). Repeat loci plotted according to log fold-change in streptavidin pull-downs (IP) versus input and significantly enriched repeats (P< 0.05) color coded according to repeat type. ALR = alpha repeat, CCCTAA= telomere repeat, BSR= beta satellite repeat. (D) Proportion of replicated transcripts in the entire short NCC RNA-seq transcriptome (all short) and in short hits. Gylling et al. 1110 RNA (2020) Vol. 26, No. 9 cell cultures, our data do not dismiss alpha repeat R-loop formation in mitosis, which has previously been reported (Kabeche et al. 2018). Also, in accordance with previous studies, only “GAA”out of the simple dinucleotide repeats found in both the long and short hits showed statistical significance for R-loop formation, indicating that the GA-rich sequences in general do not form R-loops (Grabczyk et al. 2007; Reddy et al. 2011; Groh et al. 2014). These findings were further strengthened by the genome analysis that showed surprisingly little overlap of long and short hit repeats with R-loop forming regions (Fig. 5B,C). Of the 309 nascent hit genes, 154 were identified in the R-loop data, which included half of the few replicated genes (Fig. 5D). This proportion of R-loop forming genes was expected given that half of all expressed genes from the long RNA-seq overlapped with R-loop forming regions: 3863 of 6867 expressed genes were significantly enriched in the DRIP-seq data. DISCUSSION In this work, we developed a method to investigate the interplay of RNAs with DNA replication and chromatin restoration. We used long- and short RNA sequencing and a dedicated bioinformatics pipeline to reveal that RNAs are associated with nascent and mature chromatin, followed by subsequent characterization of the nature of these RNAs (Fig. 1; Supplemental Fig. S1). The custom bioinformatics pipeline allowed us to investigate repeat reads normally discarded from analysis due to multimapping and their enrichment in different data sets such as the NCC RNA-seq, NCC DNA-seq, and DRIP-seq data. In the long RNA sequencing we identified transcripts arising from genes that are significantly enriched at nascent chromatin, albeit at moderate abundance levels (Fig. 2A; Supplemental Table S1). A subtle enrichment of a transcribed gene at a replication fork could be explained B A CD FIGURE 5. DRIP-seq analysis shows that repeat RNAs enriched in N, M2, and M10 samples do not form R-loops. (A) R-loop analysis of DRIP-seq data mapping to Repbase sequences. Repeats enriched in NCC short RNA-seq Repbase mapping (N, M2, and M10 vs. controls) are highlighted with colors according to repeat type, whereas repeats that are enriched in the DRIP-seq data set are additionally labeled with red text. ALR = alpha repeat, CCCTAA= telomere repeat. (B)Long hits split and colored according to transcript type. Bars indicate number of transcripts in long hits (all long hits) overlapping with NCC DNA-seq (replicated), DRIP-seq (R-loop forming), or both (replicated and R-loop). (C) Same as Bbut for short hits. (D) Forty-seven genes out of 309 (nascent hits) overlap with replicated regions whereas 154 are prone to R-loop formation; 22 genes are present in all three data sets (nascent hits, replicated loci, and R-loops). RNAs at replication forks www.rnajournal.org 1111 by replication-transcription collisions, which can prevent release of the transcribed product. This is also concordant with the fact that only a small proportion of active replication forks would cover a single gene during 15 min labeling in mid S-phase, making a high enrichment unlikely. However, the majority of enriched genes are not replicated at the time of labeling, suggesting that they associate with replication forks in trans (Fig. 4A). It is imaginable that the enriched RNAs associate with replication forks due to specific RNA folds common between the otherwise diverse transcripts, but we are unable to address this hypothesis due to current limitations in computational tools. Thus, our study does not support a model where specific RNAs serve general functions during active DNA replication, but does not rule out site-specific functions such as epigenetic or transcriptional regulation. Interactions in trans could be envisioned through association with proteins, triple-helix formation or R-loop formation. There is, however, no overrepresentation of nascent hit genes in R-loop data compared to all transcribed genes (Fig. 5D). We found three types of repeats associated with mid S-phase replicated DNA, associations that appear to be maintained into G1. GA-rich repeat RNAs were most strongly enriched in the long RNA-seq whereas alpha repeat RNAs were identified in the short RNA-seq and TERRA RNA in both (Figs. 2, 3; Supplemental Figs. S2, S3;Supplemental Tables S1–S4). We found the repeats and NA transcripts from the short RNA-seq to map to centromeric, pericentromeric and telomeric regions with most repeats mapping to Chr 16, 21, and Y. Repetitive and centromeric regions are difficult to sequence due to their size, repetitive nature, and polymorphism, and although hg38 does have an improved annotation of these areas compared to hg19, it is based on a “reference model”(Miga et al. 2014; Rosenbloom et al. 2015). This, and the fact that HeLa cells have a variable number of normal and abnormal chromosomes, may lead to a skewed mapping of repeats to chromosomes. The observed enrichment of repetitive RNAs could be explained by the sheer abundance of repetitive DNA in the genome. However, many other types of repeats are replicated at the time of labeling and transcribed, such as LINE and Alu repeats, which are not enriched in our NCC RNA-seq data (Supplemental Table S6). This suggests that the observed repeat RNAs are indeed retained in chromatin after being transcribed and copurifying with labeled DNA. This could be as a consequence of improper release of these specific types of repeat RNAs or interaction with chromatin either via protein or DNA binding. Alpha repeat transcription and stalling of Pol II has been suggested to contribute to chromatin establishment and function of the centromeres via interaction and recruitment of histones and histone modifiers (Bergmann et al. 2011; Chan et al. 2012; Biscotti et al. 2015; Catania et al. 2015; Blower 2016; Molina et al. 2016). TERRA RNAs have likewise been proposed to facilitate chromatin compaction and play an important role in telomere maintenance during replication, but their transcription is believed to be initiated in less heterochromatic subtelomeric regions (Azzalin et al. 2007; Schoeftner and Blasco 2008; Deng et al. 2009; Porro et al. 2014; Rippe and Luke 2015; Montero et al. 2018). We found alpha and TERRA repeat RNAs to be associated with newly replicated DNA shortly after fork passage and the association to be maintained for at least 10 h after replication (Fig. 4). This suggests that alpha and TERRA RNAs are transcribed and inefficiently released, possibly due to Pol II stalling, on centromeres, pericentromeres and telomeres shortly after chromatin disruption by the replication fork. How transcription occurs or is even initiated in heterochromatin regions and whether stalling of Pol II and transcript retention is important for both repeat RNAs functions, is not known. Alpha RNA has been shown to induce chromatin changes via R-loop formation, which is also a well-known mechanism for TERRA RNA during alternative telomere lengthening (Balk et al. 2013; Castellano-Pozo et al. 2013; Pfeiffer et al. 2013; Arora et al. 2014; Velazquez Camacho et al. 2017; Kabeche et al. 2018; Lee et al. 2018). We found alpha repeat-containing RNAs to have little overlap with R-loop data, as opposed to TERRA RNAs, suggesting that alpha repeats do not form R-loops in general (Fig. 5A–C). It is also possible that some repeat RNAs are bound in complexes with chromatin regulators that reassociate with replicated DNA during chromatin maturation. Regardless, it is intriguing that repeat RNAs are present already at a very early time point after DNA replication, arguing that they can contribute to later steps in chromatin restoration. We report a considerable enrichment for GA-rich RNA sequences in replicated DNA shortly (15 min) after replication fork passage and a maintained enrichment into G1 phase (10 h) (Fig. 4). Only a fraction of the genomic GA regions that match the sequence composition of long hit repeats and NA sequences, and overlap with replicated regions, are present in long hits (404/3802) (Supplemental Fig. S4). It is currently unknown if the remaining GA regions are transcribed or transcribed and then rapidly degraded. However, we cannot rule out that the actual number of genomic regions resembling the GA-rich hits is lower than 10,000 since our search was based on nucleic acid composition and not specific sequence motifs. A particular property of GA-rich RNA sequences is their incapacity to form stable secondary structures. However, sequences with several (four or more) stretches of Gs (may also be interspersed with other nucleotides), have the ability to form tertiary structures known as G-quadru- plexes (G4) (Gellert et al. 1962; Mukundan and Phan 2013). If loops and bulges are taken into consideration, all long hit repeat sequences are highly likely to produce G4 structures. By forming G4 structures the RNA would become more stable, which could explain why the GA-rich Gylling et al. 1112 RNA (2020) Vol. 26, No. 9