Author response: Human ORC/MCM density is low in active genes and correlates with replication time but does not delimit initiation zones (2021)

Article Figures and data Abstract Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Eukaryotic DNA replication initiates during S phase from origins that have been licensed in the preceding G1 phase. Here, we compare ChIP-seq profiles of the licensing factors Orc2, Orc3, Mcm3, and Mcm7 with gene expression, replication timing, and fork directionality profiles obtained by RNA-seq, Repli-seq, and OK-seq. Both, the origin recognition complex (ORC) and the minichromosome maintenance complex (MCM) are significantly and homogeneously depleted from transcribed genes, enriched at gene promoters, and more abundant in early- than in late-replicating domains. Surprisingly, after controlling these variables, no difference in ORC/MCM density is detected between initiation zones, termination zones, unidirectionally replicating regions, and randomly replicating regions. Therefore, ORC/MCM density correlates with replication timing but does not solely regulate the probability of replication initiation. Interestingly, H4K20me3, a histone modification proposed to facilitate late origin licensing, was enriched in late-replicating initiation zones and gene deserts of stochastic replication fork direction. We discuss potential mechanisms specifying when and where replication initiates in human cells. Introduction In human cells, DNA replication initiates from 20,000 to 50,000 replication origins selected from a five- to tenfold excess of potential or ‘licensed’ origins (Moiseeva and Bakkenist, 2018; Papior et al., 2012). Origin licensing, also called pre-replicative complex (pre-RC) formation, occurs in late mitosis and during the G1 phase of the cell cycle. During this step, the origin recognition complex (ORC) binds DNA and, together with Cdt1 and Cdc6, loads minichromosome maintenance complexes (MCM), the core motor of the replicative helicase, as inactive head-to-head double hexamers (MCM-DHs) around double-stranded DNA (Bell and Kaguni, 2013; Evrin et al., 2009; Remus and Diffley, 2009). A single ORC reiteratively loads multiple MCM-DHs. However, once MCM-DHs have been assembled, ORC does not maintain contact with the MCM-DH and neither ORC, nor Cdc6, nor Cdt1 are required for origin activation (Fragkos et al., 2015; Hyrien, 2016; Powell et al., 2015; Remus et al., 2009; Rowles et al., 1999; Sun et al., 2014; Yeeles et al., 2015). During S phase, CDK2 and CDC7 kinase activities in conjunction with other origin-firing factors convert some MCM-DHs into pairs of active CDC45-MCM-GINS helicases that nucleate bidirectional replisome establishment (Douglas et al., 2018; Moiseeva and Bakkenist, 2018). MCM-DHs that do not initiate replication are dislodged from DNA during replication. In Saccharomyces cerevisiae, origins are genetically defined by specific DNA sequences (Marahrens and Stillman, 1992). In multicellular organisms, no consensus sequence for origin activity has been identified and replication initiates from flexible locations. Although mammalian origins fire at different times through S phase, neighboring origins tend to fire at similar times, partitioning the genome into ~5,000 replication timing domains (RTDs) (Rivera-Mulia and Gilbert, 2016a). RTDs replicate in a reproducible order through S phase (Pope et al., 2014; Zhao et al., 2017). One model for this temporal regulation suggests that RTDs are first selected for initiation, followed by stochastic origin firing within domains (Boulos et al., 2015; Pope et al., 2014; Rhind and Gilbert, 2013; Rivera-Mulia and Gilbert, 2016b). A cascade or domino model suggests that replication first initiates at the most efficient (master) origins and then spreads to less efficient origins within an RTD (Boos and Ferreira, 2019; Guilbaud et al., 2011). Various processes and factors contribute to origin specification such as transcription, DNA sequences, histone variants, histone modifications, and nucleosome dynamics (Akerman et al., 2020; Cayrou et al., 2015; Long et al., 2020; Petryk et al., 2016; Prioleau and MacAlpine, 2016; Smith and Aladjem, 2014). For example, we proposed H4K20me3 to support the licensing of a subset of late-replicating origins in heterochromatin (Brustel et al., 2017). Recently, the histone variant H2A.Z has been implicated in ORC recruitment at early origins through deposition of H4K20me2 by histone methyltransferase SUV420H1 (Long et al., 2020). Furthermore, binding sites for the origin-firing factor Treslin-MTBP often feature a nucleosome-free gap adjacent to H3K4me2 (Kumagai and Dunphy, 2020). Different approaches have been developed to characterize mammalian origins. Origins have been mapped at the single-molecule level by optical methods or at the cell-population level by sequencing various purified replication intermediates, such as short nascent strands, replication bubbles, and Okazaki fragments (Hulke et al., 2020). Strand-oriented sequencing of Okazaki fragments (OK-seq) reveals the population-averaged replication fork direction (RFD) allowing to map initiation and termination (Chen et al., 2019; McGuffee et al., 2013; Petryk et al., 2016; Smith and Whitehouse, 2012). Bubble-seq (Mesner et al., 2013), single-molecule analyses (Demczuk et al., 2012; Lebofsky et al., 2006; Norio et al., 2005), and OK-seq (Petryk et al., 2016; Wu et al., 2018) studies of human cells all suggest that replication initiates in broad but circumscribed zones consisting of multiple, individually inefficient sites. OK-seq revealed both, early-firing initiation zones (IZs), which are precisely flanked on one or both sides by actively transcribed genes, and late-firing IZs distantly located from active genes (Petryk et al., 2016). Recently, an excellent agreement was observed between early-firing IZs determined by OK-seq and by EdUseq-HU, which identifies nascent DNA synthesized in early S phase cells in the presence of EdU and hydroxyurea (Tubbs et al., 2018). Furthermore, high-resolution Repli-seq identified both early and late IZs consistent with OK-seq IZs (Zhao et al., 2020). Chromatin immunoprecipitation followed by sequencing (ChIP-seq) was used to map ORC and MCM chromatin binding. In Drosophila, ORC often binds next to open chromatin marks found at transcription start sites (TSSs) (MacAlpine et al., 2010), but MCMs, initially loaded next to ORC, are more abundantly loaded and widely redistributed when cyclin E/CDK2 activity rises in late G1 (Powell et al., 2015). In human cells, ChIP-seq of single ORC subunits identified from 13,000 to 101,000 ORC potential binding sites (Dellino et al., 2013; Long et al., 2020; Miotto et al., 2016). These studies consistently demonstrated a correlation of ORC-DNA binding with TSSs, open chromatin regions, and early replication timing (RT). ChIP-seq of Mcm7 in HeLa cells suggested that MCM-DHs bind regardless of the chromatin environment, but are preferentially activated upstream of active TSSs (Sugimoto et al., 2018). We and others previously used the Epstein–Barr virus (EBV), whose replication in latency is entirely dependent on the human licensing machinery, to compare ORC and MCM binding and replication initiation sites (Chaudhuri et al., 2001; Dhar et al., 2001; Papior et al., 2012; Ritzi et al., 2003; Schepers et al., 2001). A five- to tenfold excess of potential origins were licensed per genome with respect to 1–3 mapped initiation event(s) (Norio, 2001; Norio and Schildkraut, 2004; Papior et al., 2012). These findings support the model that human replication initiates in zones, which comprise multiple, individually inefficient sites. Here, we present the first comparative survey of four different pre-RC components, replication initiation events, transcription activity, and RT in the genome of the human lymphoblastoid Raji cell line by combining ChIP-seq, OK-seq, RNA-seq, and previously published Repli-seq data (Sima et al., 2018). We find that, in pre-replicative (G1) chromatin, ORC and MCM are broadly distributed over the genome with high ORC density better correlating with early RT than MCM density. ORC/MCM are depleted from actively transcribed gene bodies and enriched at active gene promoters. ORC/MCM density is homogeneous over non-transcribed genes and intergenic regions of comparable RT. Furthermore, regions of similar RT show a similar ORC/MCM density, be they IZs, termination zones, undirectionally replicating regions (presumably lacking initiation events), or randomly replicating regions. These findings suggest that ORC/MCM densities do not solely determine IZs and that a specific contribution of the local chromatin environment is required. Indeed, we previously showed that IZs are enriched in open chromatin marks typical of active or poised enhancers (Petryk et al., 2016). We further show that a subset of non-genic late IZs is enriched in H4K20me3, confirming previous finding that H4K20me3 enhances origin activity in certain chromatin environments (Brustel et al., 2017; Shoaib et al., 2018). These findings support the cascade model for replication initiation: the entire genome (except transcribed genes) is licensed for replication initiation. Additional process and factors like adjacent active transcription and epigenetic marks are required to specify master zones of higher replication initiation efficiency. The distributed licensing pattern allows the stochastic activation of secondary origins, possibly triggered by approaching replication forks. Results Moderate averaging is a suitable approach for ORC and MCM-DH distribution analysis We used centrifugal elutriation to obtain a G1-enriched, pre-replicative population of human lymphoblastoid Raji cells (Papior et al., 2012). Propidium iodide staining followed by FACS (Figure 1—figure supplement 1a) and western blot analyses of cyclins A, B, and H3S10 phosphorylation (Figure 1—figure supplement 1b) confirmed the cell cycle stages of elutriated fractions. To ensure unbiased mapping of ORC and MCM, we simultaneously targeted Orc2, Orc3, Mcm3, and Mcm7 using validated ChIP-grade antibodies (validated in Papior et al., 2012; Ritzi et al., 2003; Schepers et al., 2001). ChIP efficiencies and qualities were measured using the EBV latent origin oriP as reference (Figure 1—figure supplement 1c). Raji cells contain 50–60 EBV episomes, allowing an easy detection of ORC/MCM at oriP (Adams et al., 1973). The viral protein EBNA1 recruits ORC to oriP’s dyad symmetry element, followed by MCM-DH loading. We detected both ORC and MCM at the dyad symmetry element in G1, whereas a population containing S-G2-M-phased cells depict a reduction in MCM levels, as expected (Figure 1—figure supplement 1c; Papior et al., 2012; Ritzi et al., 2003). ChIP-seq of two replicates for ORC subunits (Orc2, Orc3) and of three replicates for MCM proteins (Mcm3, Mcm7) resulted in reproducible, but dispersed, ChIP-seq signals as exemplified by the well-characterized replication origin Mcm4/PRKDC (Figure 1a; Ladenburger et al., 2002; Schaarschmidt et al., 2002). We employed the MACS2 peak-calling program (Feng et al., 2012; Zhang et al., 2008), but found that the obtained results were too dependent on the chosen program settings and that ORC and MCM distributions were too dispersed to be efficiently captured by peak calling (data not shown), requiring an alternative approach. Figure 1 with 3 supplements see all Download asset Open asset Moderate averaging represents a valid approach for origin recognition complex/minichromosome maintenance complex (ORC/MCM) chromatin immunoprecipitation followed by sequencing (ChIP-seq) analysis. (a) Sequencing profile visualization in UCSC Genome Browser (hg19) at the Mcm4/PRKDC origin after reads per genomic content normalization: two samples of Orc2 and Orc3, and three samples of Mcm3 and Mcm7, are plotted against the input in three replicates. The profiles are shown in a 10 kb window (chr8: 48,868,314–48,878,313); the mapped position of the origin is indicated as green line. (b) The profile of ORC/MCM ChIP-seq after 1 kb binning at the same locus. The reads of replicates were summed and normalized by the total genome-wide ChIP read frequency followed by input division. Y-axis represents the resulting relative read frequency. (c) Correlation plot between Orc2 and Orc3 relative read frequencies in 1 kb bins. (d) Correlation plot between Mcm3 and Mcm7 relative read frequencies in 1 kb bins. (e) Heatmap of Pearson correlation coefficients r between all ChIP relative read frequencies in 1 kb bins. Column and line order were determined by complete linkage hierarchical clustering using the correlation distance (d = 1 r). Refer to Figure 1—figure supplement 3 for data representation without input division. Consequently, we summed up the reads of the ChIP replicates at different binning sizes and normalized the signals against the mean read frequencies of each ChIP sample and against input, as is standard in most ChIP-seq analyses. We computed the Pearson correlation coefficients between ORC/MCM ChIPs and obtained good correlations at 1 kb bin size and only marginal improvement at larger sizes (Figure 1—figure supplement 2a). When working in 1 kb bins, we still detected the enrichment of ORC/MCM at the MCM/PRKDC origin (Figure 1b), indicating that we do not lose local, biologically relevant signals. In line with a previous report (Teytelman et al., 2009), the input control was significantly underrepresented in DNase hypersensitive (HS) regions, at TSSs, and at early RTDs (Figure 1—figure supplement 2b–d). As sonication-hypersensitive regions correlate with DNase HS regions (Schwartz et al., 2005), we carefully compared our results obtained with and without input normalization. For example, we still detect enrichment of ORC/MCM at the MCM/PRKDC origin when we omit input normalization (Figure 1—figure supplement 3). As will become apparent, similar conclusions were obtained in further analyses performed with or without input normalization. The reliability and reproducibility of our ChIP experiments is reflected by the high Pearson correlation coefficients of the relative read frequencies of Orc2/Orc3 (r = 0.866, Figure 1c) and Mcm3/Mcm7 (r = 0.879, Figure 1d). The correlations between ORC and MCM were only slightly lower (Mcm3/Orc2/3: r = 0.775/0.757; Mcm7/Orc2/3: r = 0.821/0.800, Figure 1e). Hierarchical clustering based on Pearson correlation of ChIP profiles clustered ORC and MCM profiles together. Similar results were obtained using non input-normalized data (Figure 1—figure supplement 3b–d). To compare our ChIP-seq data to previously published Orc2 ChIP-seq from asynchronously cycling K562 cells (GSE70165; Miotto et al., 2016), we calculated the relative read frequencies of our ORC ChIPs around an aggregate of K562 Orc2 peaks (>1 kb) and found substantial enrichment (Figure 1—figure supplement 3e). Miotto et al., 2016 reported that Orc2 co-localizes with DNase HS sites present at active promoters and enhancers. In line with these observations, we found a significant enrichment of ORC at DNase HS regions > 1 kb, compared to regions deprived of DNase HS sites, with or without input normalization (Figure 1—figure supplement 3f, g). These results further validate our data. ORC/MCM are enriched in IZs dependent on transcription We next compared the relative read frequencies of ORC/MCM to active replication initiation units. Using OK-seq in Raji cells (Wu et al., 2018), we calculated the RFD (see Materials and Methods) and delineated preferential replication IZs as ascending segments (ASs) of the RFD profile. RFD profiles present upshifts that define origins to kilobase resolution in yeast (McGuffee et al., 2013), but in mammalian cells these transitions are more gradual, extending over 10–100 kb (Chen et al., 2019; McGuffee et al., 2013; Petryk et al., 2016; Tubbs et al., 2018; Wu et al., 2018). We analyzed ASs > 20 kb, allowing to assess ChIP signals up to 10 kb within ASs (see Materials and Methods). Using the RFD shift across the ASs (ΔRFD) as a measure of replication initiation efficiency, we further required ΔRFD > 0.5 to select the most efficient IZs. In total, we selected 2957 ASs, with an average size of 52.3 kb, which covered 4.9% (155 Mb) of the genome (Figure 2a, green bars, Table 1). In total, 2451 (83%) of all AS located close to genic regions (ASs extended by 20 kb on both sides overlapped with at least one annotated gene). Performing RNA-seq in asynchronously cycling Raji cells, we determined that 673 ASs (22.8% of all ASs) were flanked by actively transcribed genes (transcripts per kilobase per million [TPM] >3) on both sides (type 1 AS), with less than 20 kb between AS borders and the closest transcribed gene. In total, 1026 ASs (34.7%) had only one border associated to a transcribed gene (type 2 AS; TPM >3). Also, 506 ASs (17.1%) were devoid of proximal genes (non-genic AS) (Table 1). The slope did not change considerably between the different AS types, although type 1 ASs were on average slightly more efficient, followed by type 2 ASs, then non-genic ASs (Figure 2—figure supplement 1a). Type 1 and 2 ASs located within early RTDs, while non-genic ASs were predominantly late replicating (Figure 2—figure supplement 1b), as previously observed in GM06990 and HeLa (Petryk et al., 2016). Figure 2 with 2 supplements see all Download asset Open asset Origin recognition complex/minichromosome maintenance complex (ORC/MCM) enrichment within ascending segments (ASs) depends on active transcription. (a) Top panel: example of an replication fork direction (RFD) profile on chr1: 178,400,000–182,800,000, covering 4 Mb. Detected ASs are labeled by green rectangles (irrespective of length and RFD shift). Middle and bottom panels: representative Mcm3 (blue) and Orc2 (red) chromatin immunoprecipitation followed by sequencing (ChIP-seq) profiles after binning for the same region. (b–e) Average input-normalized relative ChIP read frequencies of Orc2, Orc3, Mcm3, and Mcm7 at AS borders of (b) all AS (L > 20 kb and ΔRFD >0.5; n = 2957), (c) type 1 ASs with transcribed genes at both AS borders (n = 673), (d) type 2 ASs oriented with their AS border associated to transcribed genes at the right (n = 1026), and (e) non-genic ASs in gene-deprived regions (n = 506). The mean of ORC and MCM relative read frequencies is shown ±2 × SEM (lighter shadows). The dashed grey horizontal line indicates relative read frequency 1.0 for reference. For type 1 and 2 ASs, yellow bars mark the AS borders associated to transcribed genes. Refer to Figure 2—figure supplement 2 for analysis without input division. Table 1 Characterization of different AS subtypes. NumberGenome coverage (%)Average length (kb)All AS29574.952.3Genic AS24514.152.3Type 1 AS6731.150.7Type 2 AS10265.250.2Non-genic AS5060.850.7Only AS ≥20 kb with ΔRFD > 0.5 were considered. Genic ASs: ASs extended 20 kb on both sides is overlapped by genic region(s) irrespective of transcriptional activity; type 1 and type 2 AS: ASs flanked by expressed genes (TPM ≥3) within 20 kb on both sides (type 1) or one side (type 2); non-genic: no annotated gene ±20 kb of AS borders; AS: ascending segment; RFD: replication fork per kilobase per To the between ORC/MCM densities and replication initiation, we computed the relative read frequencies of ORC/MCM around all AS aggregate ORC and MCM on enriched within ASs compared to regions (Figure Figure 2—figure supplement without input To the of transcriptional activity, we this for the different AS (Figure data in Figure 2—figure supplement 2b–d). activity in AS regions was associated with ORC/MCM ASs Figure and a MCM from transcribed regions (Figure right In in type 2 ASs, ORC/MCM at non-transcribed AS borders (Figure ORC/MCM enrichment was detected within non-genic ASs (Figure AS borders associated with transcriptional activity were enriched in ORC/MCM (Figure both and right borders is in line with previously detected at AS borders (Petryk et al., 2016). non-genic AS borders only showed a local in ORC/MCM (Figure Figure both but the of this A sequence analysis revealed distributions of sequences at AS borders (data not sequences nucleosome and ORC but also Okazaki border detection at (Figure 2—figure supplement 1a) as as ORC and MCM are depleted from transcribed gene bodies and enriched at TSSs with previous OK-seq studies (Chen et al., 2019; Petryk et al., 2016), the average RFD profile of active genes revealed ASs upstream of TSSs and of transcriptional termination and RFD segments across the active gene bodies (Figure supplement 1a). on transcriptional activity as genes an RFD profile (Figure supplement 1a). When our ORC/MCM ChIP-seq data in to transcription, we observed that the ORC relative read distribution was significantly enriched at active TSSs, as demonstrated in (MacAlpine et al., and human cells (Dellino et al., 2013; Miotto et al., 2016). ORC relative read distribution was but significantly higher upstream of TSSs and of than within active genes (Figure These were of input normalization Figure with Figure supplement The of ORC from gene bodies was significant for of actively transcribed genes (Table 1). to ORC, Mcm3 and Mcm7 at TSSs were less but from gene bodies were more (Figure Figure supplement 1b), with and of transcribed gene bodies significantly depleted from Mcm3 and Mcm7, (Table 1). was homogeneous from to that transcription ORC and MCM-DH complexes (Figure In at genes, ORC/MCM were enriched at TSSs and were not depleted from gene bodies (Figure Figure supplement 1c). transcriptional activity as did not have on ORC/MCM at TSSs (Figure Figure supplement 1d). ORC/MCM within gene bodies was slightly more with transcription when normalized for input (Figure but this was less without input normalization (Figure supplement 1e). ORC/MCM upstream of TSSs and of were indicating that the local ORC enrichment at TSSs did not in more MCM upstream than of active genes. Figure 3 with 1 supplement see all Download asset Open asset Origin recognition complex (ORC) is enriched at active transcription start sites (TSSs) while minichromosome maintenance complex (MCM) is depleted from actively transcribed genes. ORC/MCM relative read frequencies around TSSs or transcriptional termination sites for (a) active genes (transcripts per kilobase per million [TPM] >3) and (b) inactive genes (TPM genes larger than kb without adjacent gene within kb were considered. from TSSs or are indicated in of ORC and MCM frequencies are shown ±2 × SEM (lighter shadows). The dashed grey horizontal line indicates relative read frequency 1.0 for reference. (c) ORC/MCM relative read frequencies at TSSs dependent on transcriptional activity × (d) ORC/MCM relative read frequencies upstream of TSSs and within the gene dependent on transcriptional activity × TSSs 3 kb from activity was as no (TPM (TPM (TPM and high (TPM were performed by followed by are indicated to the previous transcriptional Refer to Figure supplement 1 for analyses without input division. that Mcm3 and Mcm7 are significantly depleted from transcribed gene bodies is consistent with their active by transcription in G1, as previously proposed in (Powell et al., and human cells and 2018). process to IZs flanked by active genes. In ORC/MCM density across non-genic AS borders (Figure that ORC/MCM are not to non-genic replication IZs. ORC/MCM genomic distributions are broad and correlate with RT but not IZs RT is a of genome that is with gene and chromatin which the of origins and timing of origin firing et al., 2009). In has been reported that the of MCM-DHs loaded at origins correlates with RT profiles from stochastic origin firing et al., 2015; et al., In human cells, ORC binding data have also been used to RT profiles et al., 2016). To the between firing and ORC/MCM density in human cells, we used Raji Repli-seq data from et al., and RT to ORC/MCM relative read frequencies and RFD slope (Sima et al., 2018). We analyzed four different of RFD pattern in Figure supplement as previously defined in Petryk et al., ascending RFD segments that RFD segments that zones segments of high > over that unidirectionally replicating regions where replication in the same a of initiation and segments of RFD regions over replicating by initiation and observed in late-replicating gene deserts (Figure supplement 1c). We calculated relative Orc2 and Mcm3 (Figure Figure supplement for Orc3 and Mcm7) read frequencies in 10 kb against RT in intergenic regions gene bodies (TPM or active gene bodies (TPM right We all or to ASs, and in were normalized by that each is the probability density of ChIP frequency at a RT Figure 4 with 2 supplements see all Download asset Open asset Origin recognition complex/minichromosome maintenance complex (ORC/MCM) correlate with replication timing and transcriptional activity but are homogeneously distributed the genome and to replication fork direction (RFD) 3 × of of Orc2 (a) and Mcm3 (b) chromatin immunoprecipitation relative read frequency RT over kb to the of RT The analysis was performed in 10 kb bins. are normalized by and the probability density of ChIP relative frequencies at a replication The is indicated on are computed for intergenic regions genes (transcripts per kilobase per million TPM and expressed genes (TPM right start sites and transcriptional termination sites proximal regions were not (see Materials and The show all or to ascending replication initiation, replication unidirectionally replicating no initiation, no and RFD RFD regions, initiation and bottom The of per is indicated in each Figure supplement 1 for Orc3 and Mcm7 analyses. Refer to Figure supplement for with Figure expressed genes showed lower ORC/MCM densities than genes and intergenic regions (Figure Figure supplement was significant in early- and regions, as demonstrated by (Figure supplement 2a, The was more for MCM than ORC, as in Figure In the difference between genes and intergenic regions was at significant In all ORC/MCM densities from early to late RT but this RT was in expressed genes, for MCM, as