Weitao Wang

Kyle N Klein

Karel Proesmans

Hongbo Yang

Claire Marchal

Xiaopeng Zhu

Tyler Borrman

Alex Hastie

Zhiping Weng

John Bechhoefer
email: johnb@sfu.ca

Chun-Long Chen
email: chunlong.chen@curie.fr

David M Gilbert
email: gilbert@bio.fsu.edu

Nicholas Rhind
email: nick.rhind@umassmed.edu

Genome-wide mapping of human DNA replication by optical replication mapping supports a stochastic model of eukaryotic replication

come   L'archive ouverte pluridisciplinaire

Introduction

The human genome is replicated from approximately 50,000 distinct initiation events in each cell cycle [START_REF] Masai | Eukaryotic chromosome DNA replication: where, when, and how?[END_REF]. However, identifying the location and firing timing of human replication origins is an ongoing challenge [START_REF] Hyrien | Peaks cloaked in the mist: the landscape of mammalian replication origins[END_REF]. A major part of this challenge is the inefficient and heterogeneous nature of mammalian origin firing [START_REF] Berezney | Heterogeneity of eukaryotic replicons, replicon clusters, and replication foci[END_REF][START_REF] Taylor | Increase in DNA replication sites in cells held at the beginning of S phase[END_REF]. Even in cases where mammalian origins are thought to occur at relatively well-defined sites, they appear to fire with as little as 10% efficiency [START_REF] Dijkwel | Initiation sites are distributed at intervals in the Chinese hamster dihydrofolate reductase origin of replication but are used with very different efficiencies[END_REF]; in the more common case where initiation appears to occur within broad initiation zones (IZs), spanning tens of kilobases [START_REF] Hamlin | A revisionist replicon model for higher eukaryotic genomes[END_REF], they fire with less than 1% efficiency [START_REF] Demczuk | Regulation of DNA Replication within the Immunoglobulin Heavy-Chain Locus During B Cell Commitment[END_REF]. Therefore, population-based bulk approaches to mapping origins, which average origin-firing behavior across large populations of cells, must contend with low signal-tonoise ratios. As a result, although a number of bulk approaches have been used to map origins in the human genome [START_REF] Mesner | Bubbleseq analysis of the human genome reveals distinct chromatin-mediated mechanisms for regulating early-and late-firing origins[END_REF][START_REF] Besnard | Unraveling cell type-specific and reprogrammable human replication origin signatures associated with G-quadruplex consensus motifs[END_REF][START_REF] Cayrou | Genome-scale analysis of metazoan replication origins reveals their organization in specific but flexible sites defined by conserved features[END_REF][START_REF] Foulk | Characterizing and controlling intrinsic biases of lambda exonuclease in nascent strand sequencing reveals phasing between nucleosomes and G-quadruplex motifs around a subset of human replication origins[END_REF][START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF][START_REF] Macheret | Intragenic origins due to short G1 phases underlie oncogeneinduced DNA replication stress[END_REF][START_REF] Dellino | Genome-wide mapping of human DNAreplication origins: levels of transcription at ORC1 sites regulate origin selection and replication timing[END_REF][START_REF] Petryk | Replication landscape of the human genome[END_REF], there is little concordance among them [START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF][START_REF] Mesner | Bubbleseq analysis of the human genome reveals distinct chromatin-mediated mechanisms for regulating early-and late-firing origins[END_REF].

The heterogeneous nature of mammalian origins may be influenced by the non-site-specific nature of the mammalian origin recognition complex (ORC) DNA binding. Unlike budding yeast ORC, which binds a specific 17-bp sequence [START_REF] Theis | The ARS309 chromosomal replicator of Saccharomyces cerevisiae depends on an exceptional ARS consensus sequence[END_REF][START_REF] Eaton | Conserved nucleosome positioning defines replication origins[END_REF], mammalian ORC shows only a weak affinity for AT-rich sequences in vitro [START_REF] Vashee | Sequence-independent DNA binding and replication initiation by the human origin recognition complex[END_REF][START_REF] De Carli | High-Throughput Optical Mapping of Replicating DNA[END_REF] and no discernible sequence preference in vivo [START_REF] Miotto | Selectivity of ORC binding sites and the relation to replication timing, fragile sites, and deletions in cancers[END_REF]. G-quadruplex (G4) sequences have been proposed to function as origins [START_REF] Cayrou | The chromatin environment shapes DNA replication origin organization and defines origin classes[END_REF][START_REF] Valton | G4 motifs affect origin positioning and efficiency in two vertebrate replicators[END_REF][START_REF] Prorok | Involvement of G-quadruplex regions in mammalian replication origin activity[END_REF], but this hypothesis is controversial [START_REF] Foulk | Characterizing and controlling intrinsic biases of lambda exonuclease in nascent strand sequencing reveals phasing between nucleosomes and G-quadruplex motifs around a subset of human replication origins[END_REF][START_REF] Miotto | Selectivity of ORC binding sites and the relation to replication timing, fragile sites, and deletions in cancers[END_REF]. A simple explanation for the observed patterns of ORC binding and replication initiation is that ORC binds non-specifically to accessible, nucleosome-free DNA [START_REF] Lubelsky | Pre-replication complex proteins assemble at regions of low nucleosome occupancy within the Chinese hamster dihydrofolate reductase initiation zone[END_REF]. In particular, ORC binding and replication initiation correlate well with DNase I accessibility across metazoan genomes [START_REF] Miotto | Selectivity of ORC binding sites and the relation to replication timing, fragile sites, and deletions in cancers[END_REF]Gindin et al., 2014b). Furthermore, MCMs can slide after loading, further dispersing sites of initiation [START_REF] Gros | Post-licensing Specification of Eukaryotic Replication Origins by Facilitated Mcm2-7 Sliding along DNA[END_REF][START_REF] Edwards | MCM2-7 complexes bind chromatin in a distributed pattern surrounding the origin recognition complex in Xenopus egg extracts[END_REF][START_REF] Pourkarimi | Spatiotemporal coupling and decoupling of gene transcription with DNA replication origins during embryogenesis in C. elegans. Elife 5[END_REF]. For these reasons, we will refrain from referring to mammalian sites of replication initiation as origins but rather refer to specific instances of replication initiation as initiation events occurring at initiation sites and to areas enriched in initiation sites as initiation zones.

Although the specific locations of replication initiation are unclear, the general characteristics of mammalian sites of initiation have been described. Initiation is enriched in euchromatic promoter and enhancer regions, consistent with its correlation with accessible chromatin [START_REF] Ganier | Metazoan DNA replication origins[END_REF][START_REF] Petryk | Replication landscape of the human genome[END_REF][START_REF] Pourkarimi | Spatiotemporal coupling and decoupling of gene transcription with DNA replication origins during embryogenesis in C. elegans. Elife 5[END_REF][START_REF] Cayrou | The chromatin environment shapes DNA replication origin organization and defines origin classes[END_REF]. Although initiation tends to be enriched in euchromatic regions, initiation must also occur in heterochromatin. Initiation in heterochromatin appears to be even more heterogeneous, making heterochromatic initiation sites even less well understood [START_REF] Petryk | Replication landscape of the human genome[END_REF][START_REF] Cayrou | The chromatin environment shapes DNA replication origin organization and defines origin classes[END_REF].

Much effort has also been put into understanding the relative timing of DNA replication across the genome [START_REF] Rhind | DNA replication timing[END_REF][START_REF] Fragkos | DNA replication origin activation in space and time[END_REF][START_REF] Marchal | Control of DNA replication timing in the 3D genome[END_REF]. The replication timing of a genomic region is primarily determined by the timing of initiation within that region. Based on population replication timing data, metazoan genomes can be divided into regions of similar local replication timing referred to as constant timing regions [START_REF] Hiratani | Global reorganization of replication domains during embryonic stem cell differentiation[END_REF][START_REF] Farkash-Amar | Global organization of replication time zones of the mouse genome[END_REF][START_REF] Hansen | Sequencing newly replicated DNA reveals widespread plasticity in human replication timing[END_REF]. Constant timing regions consist of smaller units whose replication timing can be regulated during development; these smaller units are called replication domains [START_REF] Hiratani | Global reorganization of replication domains during embryonic stem cell differentiation[END_REF]. Replication domains tend to align with topologically associated domains (TADs) and are thought to represent discrete functional units of the genome [START_REF] Pope | Topologically associating domains are stable units of replication-timing regulation[END_REF]. The replication timing of replication domains correlates with the level of expression of the genes they contain and their chromatin modifications. However, whether the replication timing is a cause, a consequence, or an independent correlate of chromatin structure is a matter of debate [START_REF] Rhind | DNA replication timing[END_REF].

In budding yeast, where replication origin locations are well defined, origin firing is stochastic, with each origin firing with a specific probability, independent of neighboring origins [START_REF] Czajkowsky | DNA Combing Reveals Intrinsic Temporal Disorder in the Replication of Yeast Chromosome VI[END_REF][START_REF] Yang | Modeling genome-wide replication kinetics reveals a mechanism for regulation of replication timing[END_REF][START_REF] De Moura | Mathematical modelling of whole chromosome replication[END_REF]. Such stochastic firing leads to reproducible replication profiles at the population level, because more efficient origins are more likely to fire early and therefore, on average, have early replication times; by contrast, inefficient origins usually fire late or are passively replicated [START_REF] Rhind | Reconciling stochastic origin firing with defined replication timing[END_REF]. The heterogeneous and inefficient nature of metazoan replication initiation is also consistent with stochastic initiation. Simulations with a stochastic firing model, in which initiation is regulated only by a local!firing-probability function, faithfully reproduce experimental genome-wide replication timing profiles, suggesting that no deterministic timing program is required (Gindin et al., 2014b). Furthermore, if the local initiation rate is predicted by DNase I hypersensitivity, the simulation closely matches experimental results, consistent with the observed correlation between promoters and enhancers, which are DNase I hypersensitive, and initiation frequency (Gindin et al., 2014b). On the other hand, the reproducible replication timing of individual replication domains measured in single cells has led to the suggestion that replication within those domains initiates at defined times in most cells in the population [START_REF] Dileep | Single-cell replication profiling to measure stochastic variation in mammalian replication timing[END_REF]. Furthermore, neighboring initiation sites have been proposed to show both cooperative firing and lateral inhibition [START_REF] Cayrou | Genome-scale analysis of metazoan replication origins reveals their organization in specific but flexible sites defined by conserved features[END_REF][START_REF] Guilbaud | Evidence for Sequential and Increasing Activation of Replication Origins along Replication Timing Gradients in the Human Genome[END_REF][START_REF] Löb | 3D replicon distributions arise from stochastic initiation and domino-like DNA replication progression[END_REF], neither of which are consistent with strictly stochastic models. Therefore, whether metazoan initiation timing is stochastic or deterministic, or some combination of the two, is still very much an open question [START_REF] Bechhoefer | Replication timing and its emergence from stochastic processes[END_REF].

A powerful solution to the problems of heterogeneity and low signal-to-noise ratios is singlemolecule analysis, which allows the identification of sites of replication initiation on individual DNA fibers [START_REF] Técher | Replication dynamics: biases and robustness of DNA fiber analysis[END_REF]. Traditional single-molecule approaches-such as fiber autoradiography [START_REF] Huberman | On the mechanism of DNA replication in mammalian chromosomes[END_REF], DNA combing [START_REF] Herrick | Single molecule analysis of DNA replication[END_REF][START_REF] Anglana | Dynamics of DNA replication in mammalian somatic cells: nucleotide pool modulates origin choice and interorigin spacing[END_REF][START_REF] Kaykov | Molecular Combing of Single DNA Molecules on the 10 Megabase Scale[END_REF] and SMARD [START_REF] Norio | Visualization of DNA replication on individual Epstein-Barr virus episomes[END_REF]-have provided critical insight into the location and firing kinetics of mammalian replication origins. However, these techniques are restricted to the analysis of at most a few genomic loci. Nanopore-sequencing-based methods offer the potential to map replication genome-wide at high resolution [START_REF] Müller | Capturing the dynamics of genome replication on individual ultra-long nanopore sequence reads[END_REF][START_REF] Georgieva | Detection of Base Analogs Incorporated During DNA Replication by Nanopore Sequencing[END_REF][START_REF] Hennion | FORK-seq: replication landscape of the Saccharomyces cerevisiae genome by nanopore sequencing[END_REF], but the current throughput of nanopore sequencing does not allow genome-wide analysis of metazoan genomes and current average read lengths are on the order of only 30 kb. Optical mapping technologies provide ample throughput on individual long (150 kb to 2 Mb) DNA fibers [START_REF] Lam | Genome mapping on nanochannel arrays for structural variation analysis and sequence assembly[END_REF] and can be combined with nucleotide analog incorporation to map DNA replication [START_REF] De Carli | High-Throughput Optical Mapping of Replicating DNA[END_REF]. However, the current approach relies on DNA replication in extracts [START_REF] De Carli | High-Throughput Optical Mapping of Replicating DNA[END_REF], precluding in vivo analysis of replication initiation.

We have developed Optical Replication Mapping (ORM), a single-molecule technique to investigate spatial and temporal distribution of origin firing that combines the Bionano Genomics approach to mapping long individual DNA molecules [START_REF] Lam | Genome mapping on nanochannel arrays for structural variation analysis and sequence assembly[END_REF] with in vivo fluorescent nucleotide pulselabeling [START_REF] Wilson | Spatio-temporal re-organization of replication foci accompanies replication domain consolidation during human pluripotent stem cell lineage specification[END_REF][START_REF] Panning | Spatio-temporal organization of DNA replication in murine embryonic stem, primary, and immortalized cells[END_REF] to directly visualize sites of replication initiation within human cells. This approach affords us excellent signal-to-noise characteristics and deep, genomewide coverage, allowing us to identify initiation sites active in as few as 0.1% of human cells. We have used ORM to identify and analyze human sites of DNA replication initiation. We have further applied ORM to asynchronous human cells to obtain cell-type-specific replication profiles.

Results

Optical Replication Mapping of Early-Firing Human Initiation Sites

We mapped sites of early replication in the human genome by ORM (Figure 1A,B and Methods). HeLa S3 cells were synchronized in mitosis by mitotic shake off and arrested at the beginning of S phase with aphidicolin. This protocol allows early initiation events but prevents elongation of replication forks from these initiation events and also prevents subsequent initiation events by activation of the intra-S-phase checkpoint [START_REF] Feijoo | Activation of mammalian Chk1 during DNA replication arrest: a role for Chk1 in the intra-S phase checkpoint monitoring replication origin firing[END_REF]. Cells were then electroporated with fluorescent ATTO647-dUTP, released from aphidicolin, allowed to complete replication, and harvested for analysis on the Bionano Saphyr platform. As a result, cells incorporate fluorescent nucleotides in replication forks that initiated during the aphidicolin arrest and that initiated immediately after release. However, the transfected nucleotides are rapidly depleted, preventing incorporation at later initiation sites [START_REF] Wilson | Spatio-temporal re-organization of replication foci accompanies replication domain consolidation during human pluripotent stem cell lineage specification[END_REF]. In our initial experiment, 87% of cells were synchronized in mitosis and, after release from aphidicolin, 94% incorporated the fluorescent label, with 100% of labeled cells showing an early-replication pattern of replication foci (Table S1, Figures S1A-C), indicating a high degree of synchrony and <1% of contaminating asynchronous S phase cells.

We collected a dataset of 728 Gb across 2.9 million fibers ranging in length from 150 kb to 1.2 Mb, with an average fiber length of 250 kb and an average genome coverage of 206 fold (Table S1 and Figure S1D). On these fibers, we identified fluorescent incorporation signals and mapped them across the genome (Figure 1A, see Methods for details). Technical replicates within our first dataset and three additional biological replicates all show high reproducibility, with the variation dominated by counting noise in infrequently labeled regions of the genome (r = 0.81-0.86, Figures S1E,F). Therefore, we combined the four experiments into one dataset of 11.8 million fibers ranging from 150 kb to 2.2 Mb in length, averaging 285 kb in length and constituting 1,109-fold coverage of the human genome (Table S1). The coverage is uniform, with a coefficient of variation of 12% (Figure S1G).

The replication-incorporation signal we detect is punctate in nature and uniform in intensity (Figures 1A and S2A). From these observations, we conclude that these incorporation signals consist of individual fluorescent nucleotides. Furthermore, the distributions of signal intensity (Figure S2A) and inter-signal distances (Figure S2B) are consistent with sparse labeling by a pool of fluorescent nucleotides that is depleted by replication, as expected for a one-time bolus of nucleotides delivered by electroporation. The distribution of the signals is predicted by inefficient, random incorporation of fluorescent nucleotides at an initial rate of about 1/900 thymidines that is depleted with a half length of about 75 kb (Figure S2B). Given the lengths of the replication tracks and the sparse nature of the signal within them, we estimate that we can localize initiation events with a resolution of about 15 kb (14.2±0.1 kb, see Methods for details).

The fluorescent signals are enriched in replication IZs, as identified in HeLa cells by replication timing peaks (Figure 1B). To test whether the incorporation signals truly reflected initiation events, we modified our protocol to allow for a delay, ranging from 5' to 90', after aphidicolin release and before transfection (Table S1). When we aggregate replication signals around early IZs, we see a time!dependent movement of incorporation tracks away from the IZs at an estimated rate of 1.65±0.31 kb per minute (Figure 1C), consistent with prior estimates of mammalian replication fork rates [START_REF] Jackson | Replicon clusters are stable units of chromosome structure: evidence that nuclear organization contributes to the efficient activation and propagation of S phase in human cells[END_REF][START_REF] Conti | Replication fork velocities at adjacent replication origins are coordinately modified during DNA replication in human cells[END_REF][START_REF] Chagin | 4D Visualization of replication foci in mammalian cells corresponding to individual replicons[END_REF]. This result is consistent with the signal in our 0' data representing initiation events and the signals after various release times representing forks moving away from such initiation events.

To facilitate visualization of labeled replication forks, we segmented the fluorescent-nucleotide incorporation signal into discrete replication tracks, identifying replication forks active during our labeling period (see Methods for segmentation details). In our combined 0' dataset, we observed 977,746 replication tracks averaging 19.5±30.9 kb in length. The observed average track length of 20 kb is shorter than the 75 kb that we estimate is replicated during the labeling period (Figure S2B) due to the sparsity of the labeling. In particular, 58% of our replication tracks are labeled with only one signal. To test whether these solo signals are sparsely labeled replication tracks or spurious background noise, we calculated the enrichment of solo signals at early replication sites. Solo signals are enriched at such sites, demonstrating that they, too, mark regions of replication initiation (Figure S2C), and that noise, if any, is very low with this method. Although many of the tracks are labeled with only one signal, few are not labeled at all. Our modeling of the rate of label incorporation and label depletion suggest that over 90% of initiation events in early S-phase will incorporate at least one label within 15 kb of replication (Figures S2A and B and Supplemental Mathematical Methods Section 4). Therefore, almost every early-initiation event will be visualized at a resolution of 15 kb.

For fibers on which we identify more than one track, the average distance between track centers is 111±78 kb (Figure S2D), consistent with previous measurements [START_REF] Cayrou | Genome-scale analysis of metazoan replication origins reveals their organization in specific but flexible sites defined by conserved features[END_REF][START_REF] Jackson | Replicon clusters are stable units of chromosome structure: evidence that nuclear organization contributes to the efficient activation and propagation of S phase in human cells[END_REF]. We find enrichment of replication tracks around previously identified sites of human replication initiation (Figure 1D,E). In particular, we examined the distribution of replication tracks around the MYC and Top1 loci, both characterized as early!firing origins in the human genome [START_REF] Tao | Major DNA replication initiation sites in the c-myc locus in human cells[END_REF][START_REF] Keller | The origin recognition complex marks a replication origin in the human TOP1 gene promoter[END_REF], and found a pronounced enrichment at these loci, particularly in the intergenic regions surrounding these genes. More generally, replication tracks are enriched in the earliest replicating parts of the genome, as defined by replication timing profiling [START_REF] Hansen | Sequencing newly replicated DNA reveals widespread plasticity in human replication timing[END_REF], as expected for early-initiation events (Figures 1F,G). Interactive display of the ORM data is available on the ORM Browser <http://orm.nucleome.org> (Figure S3, Zhu et al. in prep, personal communication, J. Ma).

We observe 0.97 million initiation events in our combined 0' dataset. From the median incorporation track length of 75 kb (estimated from our modeling of label-incorporation kinetics, which compensates for the sparseness of labeling, Figure S2B), we infer that 73 Gb of DNA is labeled, which is ~2% of the 3.36 Tb in the dataset. Therefore, we are labeling the DNA replicated in the first ~2% of S phase and during that time, about 1000 initiation events occur per genome (977,746 initiation events observed in 1,109 genome equivalents).

Determination of Genome-Wide Replication Kinetics in Asynchronous Cells

The experiments described above provide unprecedented single-molecule mapping of initiation events at the onset of S phase. However, most cell types are not amenable to such precise cell-cycle synchronization. To map replication kinetics in unperturbed cells, we labeled asynchronous HeLa cells and H9 human embryonic stem cells, and mapped replication incorporation tracks by ORM (Figure 2A). We identified 412,113 replication tracks in two biologically-independent HeLa replicates totaling 1.4 Tb of data and 299,595 tracks in one H9 dataset totaling 738 Gb of data (Table S1). Using the same analysis as for our synchronous datasets, we infer a similar nucleotide-labeling frequency of 1/1025 and 1/850 thymidines, respectively. The replication tracks average 23.9±35.5 kb in length in the HeLa data and 27.5±40.4 kb in length in the H9 data, which is comparable to the length of tracks in the synchronized data. Thus, forks released from aphidicolin arrest synthesize at about the same rate as untreated replication forks. The tracks are uniformly distributed across the genome, as predicted for asynchronous replication forks, with an average density of 1.3±0.5%. In particular, in contrast to our synchronous dataset, and as expected, we see no enrichment at replication timing peaks in early-or late-replicating regions (Figure 2B).

In order to use this data to infer replication kinetics, it is necessary to determine the polarity of each incorporation track, so that the direction of each mapped replication fork is known. To determine track polarity, we took advantage of the fact that, as the labeled nucleotides are consumed, the label density diminishes, a strategy that has been used in other fiber!analysis approaches [START_REF] Müller | Capturing the dynamics of genome replication on individual ultra-long nanopore sequence reads[END_REF][START_REF] Huberman | On the mechanism of DNA replication in mammalian chromosomes[END_REF][START_REF] Hennion | FORK-seq: replication landscape of the Saccharomyces cerevisiae genome by nanopore sequencing[END_REF]. We calculated a fork direction index (FDI) as the ratio of integrated fluorescent signal in the left half of the track to that in the right half, for all tracks with 3 or more signals.

As expected, ORM replication tracks tend to be oriented in the direction of replication inferred from replication timing profiles and Okazaki fragment mapping [START_REF] Hansen | Sequencing newly replicated DNA reveals widespread plasticity in human replication timing[END_REF][START_REF] Petryk | Replication landscape of the human genome[END_REF] (Figures 2C,D). Furthermore, genome-wide analysis of replication fork polarity around replication timing peaks shows that the polarity signal in ORM data is comparable to that in published OK-seq data (r = 0.59, Figures 2D,E). Importantly, the polarity signal in ORM data is cell-type specific (Figure 2F). These results

show that ORM data can be used to characterize replication kinetics in unsynchronized cells, demonstrating its applicability to any cells that can be pulse labeled with fluorescent nucleotides.

Genome-Wide Mapping of Early-Firing Human Initiation Zones

To map the genome-wide distribution of early-firing human initiation sites, we analyzed the incorporation signal in our combined 0' dataset and identified 4,930 IZs, defined as peaks of replication signal density (Figure 3A, see Methods for details). The IZs are mostly between 20 and 40 kb, with an average length of 32.3±15.9 kb (Figure S2E), comparable to IZs mapped by OK-seq [START_REF] Petryk | Replication landscape of the human genome[END_REF]. The initiations zones cover about 5.3% of the genome and contain about 19.8% of the ORM signal. Much of the rest of the ORM signal is directly adjacent to defined IZs, but much is distributed elsewhere at lower density, particularly in late-replicating regions on the genome (Figure 3A). Plotting the initiation efficiency of 50 kb genomic bins show that there is a continuous distribution of initiation efficiency between the higherefficiency regions we define as IZs and lower-efficiency regions, with no obvious break between them (Figure S2F).

The median genome-wide inter-initiation-zone distance is 273 kb. However, since our dataset is biased towards early replicating IZs, we suspect that we underestimate the number of IZs in latereplicating parts of the genome. Consistent with that contention, limiting our analysis to the earliest quarter of the genome leads to a median inter-initiation-zone distance of 168 kb. The inter-initiation-zone distances in the other three quarters of the genome are 257, 258 and 337 kb, respectively. We suspect that the lower density of observed IZs in the later-replicating regions of the genome results from having less data rather than fewer initiation events in these regions (Figure 1G).

A major advantage of ORM is that it records the number of fibers that both contain and do not contain replication tracks, allowing for the calculation of the absolute efficiency with which replication initiates at each locus in the genome. Although such an analysis is possible with prior methods, it is considerably more difficult and thus rarely reported [START_REF] Demczuk | Regulation of DNA Replication within the Immunoglobulin Heavy-Chain Locus During B Cell Commitment[END_REF]. Such an analysis shows that early initiation is quite heterogeneous. Examining the IZs with average replication times in the first quarter of S phase (S50 < 25%), we find that they initiate replication, on average, on only 7.5% of molecules (Figure 3B). Even the 5% earliest IZs are active on only 11% of molecules. We therefore conclude that, in aphidicolin arrested cells, replication initiates in only a stochastically-selected subset of IZs, including only 11% of the most efficient early IZs.

We investigated the frequency of early initiation of replication across the genome, examining all 4,930 ORM IZs. We find that the local efficiency of initiation in these zones is generally below 15%, with occasional efficiencies up to 20% and a median efficiency of 3.7% (Figure 3B). All of the IZs with efficiencies of more than 10%, and 99% of the IZs with efficiencies of more than the median of 3.7%, are in early-replicating regions of the genome (Figure 3B), consistent with the results in Figure 1G. We also find a strong concordance between the IZs identified by ORM and those identified by OK-seq [START_REF] Petryk | Replication landscape of the human genome[END_REF]. In particular, 68% of the 6,166 early IZs identified by OK-seq show at least 3% initiation efficiency by ORM and 100% of loci that show greater than 10% efficiency by ORM are identified by OK-seq as IZs (Figure 3C).

Visual inspection of the distribution of replication tracks suggests that initiation in human cells is not concentrated at discrete, frequently-used origins. Instead, initiation appears to be distributed across broad IZs, as has been observed at specific loci and in other genome-wide analyses [START_REF] Petryk | Replication landscape of the human genome[END_REF][START_REF] Mesner | Bubbleseq analysis of the human genome reveals distinct chromatin-mediated mechanisms for regulating early-and late-firing origins[END_REF][START_REF] Besnard | Unraveling cell type-specific and reprogrammable human replication origin signatures associated with G-quadruplex consensus motifs[END_REF][START_REF] Hamlin | A revisionist replicon model for higher eukaryotic genomes[END_REF][START_REF] Dijkwel | Initiation sites are distributed at intervals in the Chinese hamster dihydrofolate reductase origin of replication but are used with very different efficiencies[END_REF][START_REF] Demczuk | Regulation of DNA Replication within the Immunoglobulin Heavy-Chain Locus During B Cell Commitment[END_REF][START_REF] Lubelsky | Pre-replication complex proteins assemble at regions of low nucleosome occupancy within the Chinese hamster dihydrofolate reductase initiation zone[END_REF][START_REF] Tao | Major DNA replication initiation sites in the c-myc locus in human cells[END_REF]. Even the IZ at the Top1 locus appears to have initiation events distributed over tens of kilobases (Figures 1E,S4A). However, given the sparse labeling of our replication tracks, we cannot identify initiation sites with resolution higher than about 15 kb. Nonetheless, 93.5% of IZs contain at least two incorporation tracks with centers that are more than 15 kb apart, suggesting that these IZs contain at least two, and possibly many more, initiation sites (Figures 3D and S4A,B). Further, the 6.5% of IZs that do not have obviously non-overlapping replication initiation sites tend to be smaller (16 kb v. 33 kb), less efficient (1.2% v. 4.5%) and later replicating (S50 0.55 v. 0.32) than those that do (Figure S4B). Furthermore, a reanalysis of 66 potentially discrete initiation site identified by OK-seq [START_REF] Petryk | Replication landscape of the human genome[END_REF] is consistent with the conclusion that there are no isolated, discrete initiation sites in the human genome (Figure S5).

Another way to distinguish between initiations that occur at a single site and those that are distributed across a broad region is to examine the distribution of ORM signals across the IZs, the expectation being that a single initiation site will have narrower distribution than distributed initiation sites. We measured the signal at the edges of IZs as a fraction of the signal at the middle of the IZ and fit to that data two models, one based on a single initiation site and another based on a uniform distribution of initiation sites across the IZ (Figure S4C). For IZs less that 55 kb long, which includes over 90% of IZs, the distribution signal is consistent with a uniform distribution of initiation sites throughout the IZs. Taking all of our data together, we find no evidence of isolated or high-efficiency initiation sites in the human genome.

Initiation Events in Neighboring Initiation Zones are not Clustered

Regardless of how initiation events are locally distributed, it is important to know whether local initiation events are coordinated. With our deep, single-molecule data on long fibers, we can quantitatively test whether initiation events in neighboring IZs are correlated. A simple, unbiased test for correlation between events is to determine whether the spacing between the event deviates from an exponential distribution. Items distributed randomly on a line have an exponential distribution of inter-item distances [START_REF] Birnbaum | Statistical Methods for Poisson Processes and Exponential Populations[END_REF]; therefore, correlated clusters of initiation events would deviate from an exponential distribution of inter-initiation distances. We measured 18,275 inter-initiation distances on 13,769 fibers containing multiple initiation events (Figure 3E). We measure few initiation events less than 60 kb apart, due to the facts that the tracks themselves average 40 kb in length and that early IZs have a median distance of 168 kb. However, above 60 kb, there is an exponential distribution of inter-initiation distances (r = 0.99, Figure 3E, inset). This distribution indicates that, on this length scale, initiation events are distributed randomly.

As an orthogonal test of whether initiation events in neighboring IZs are positively or negatively correlated, we compared the frequency of labeling at each IZ with the frequency with which its neighboring IZs were labeled on the same fiber (Figure 3F). Neighboring IZs can be co-labeled as a result of individual initiation in each zone or of initiation in one zone leading to passive incorporation of label in a neighboring zone by the same replication fork. Correlation can also arise because of heterogeneous uptake of labeled nucleotides. Basically, if one cell has a very high level of nucleotide uptake, all IZs (and any two neighboring IZs, in particular) will be more likely to incorporate label, making them appear to be correlated (see Methods for a more detailed explanation). Therefore, we fit a model to our initiation-zone co-labeling data that takes into account the previously determined incorporation track length (Figure S2B) and nucleotide uptake heterogeneity (Figure S1C), and assumes that individual initiation events are independent. This model fits the data well, suggesting there is no measurable positive or negative correlation between initiation in neighboring zones (Figure 3F).

Comparison of Different Origin Mapping Datasets

We compared our optical replication mapping data to four published genome-wide HeLa replicationinitiation-mapping datasets: OK-seq [START_REF] Petryk | Replication landscape of the human genome[END_REF], SNS-seq [START_REF] Picard | The spatiotemporal program of DNA replication is associated with specific combinations of chromatin marks in human cells[END_REF], Ini-seq [START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF] and Orc1 ChIP-seq [START_REF] Dellino | Genome-wide mapping of human DNAreplication origins: levels of transcription at ORC1 sites regulate origin selection and replication timing[END_REF]. By visual inspection, replication tracks appear to be enriched around the IZs identified by OK-seq (Figure 3A). Across the genome, we see colocalization of OK-seq IZs and ORM signal. We see similar colocalization of ORM signal with Ini-seq signal, a cell-free initiation mapping approach and, to a lesser extent, Orc1 ChIP-seq peaks and SNS-seq peaks, which maps initiation events by sequencing short nascent strands produced by replication initiation (Figure 3A). To determine whether the apparent co!"ocalization of initiation mapping data is robust, and to quantify its extent, we measured the correlation between the five data sets (Figure 3G). We find that ORM replication tracks correlate well Ini-seq (r = 0.59), to a lesser extent with OK-seq (r = 0.49), and even less well with SNS-seq (r = 0.36) and Orc1 ChIP-seq (r = 0.31). These correlations are further confirmed by ROC analysis (Figure S6A).

Characterization of Initiation Zones Mapped by ORM

To characterize the genomic character and chromatin context of our IZs, we investigated the enrichment of various genomic and epigenomic annotations, all derived from HeLa cells (Figure 4). Comparing our IZs with various chromatin states from the ENCODE project [START_REF] Bernstein | An integrated encyclopedia of DNA elements in the human genome[END_REF][START_REF] Ernst | ChromHMM: automating chromatin-state discovery and characterization[END_REF], we find that they are enriched in enhancers and low activity intergenic regions and depleted in transcription units and polycomb-repressed chromatin (Figure S7A). Examining specific histone modifications, our IZs are enriched in H2AZ, which has been implicated in ORC binding, and enhancerspecific modifications, such as H3K4me1 and H3K27ac, and they are depleted in, although adjacent to, transcription-elongation-associated modifications, such as H3K79me2 and H3K36me2 (Figures 4 andS7B,D), in agreement with previous studies [START_REF] Petryk | Replication landscape of the human genome[END_REF][START_REF] Long | H2A.Z facilitates licensing and activation of early replication origins[END_REF][START_REF] Pourkarimi | Spatiotemporal coupling and decoupling of gene transcription with DNA replication origins during embryogenesis in C. elegans. Elife 5[END_REF]. Consistent with the enrichment of initiation events in enhancer regions and their depletion in transcribed regions, IZs are enriched in DNase I hypersensitive sites and depleted in RNA pol II (Figures 4 andS7B,D).

Although metazoan ORC has only a weak affinity for AT-rich sequences in vitro [START_REF] Vashee | Sequence-independent DNA binding and replication initiation by the human origin recognition complex[END_REF][START_REF] De Carli | High-Throughput Optical Mapping of Replicating DNA[END_REF] and no discernible sequence preference in vivo [START_REF] Miotto | Selectivity of ORC binding sites and the relation to replication timing, fragile sites, and deletions in cancers[END_REF], sequences that form G4 quadruplex structures have been proposed to correlate with initiation sites in vertebrate cells [START_REF] Cayrou | The chromatin environment shapes DNA replication origin organization and defines origin classes[END_REF][START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF][START_REF] Valton | G4 motifs affect origin positioning and efficiency in two vertebrate replicators[END_REF][START_REF] Prorok | Involvement of G-quadruplex regions in mammalian replication origin activity[END_REF]. We tested whether IZs are enriched in G4 motifs, using a GC-content-adjusted background model to avoid detecting G4 enrichment as a trivial consequence of the GC-rich nature of enhancers. Of the 737,735 G4 sequences computationally identified in the human genome [START_REF] Lombardi | Thermodynamically stable and genetically unstable G-quadruplexes are depleted in genomes across species[END_REF], 68,585 are found in IZs, which is not more than would be expected by chance, given their GC-rich nature (partial correlation r = -0.0513). Moreover, GC-rich regions (with similar GC% as IZs) outside of IZ are no less likely to contain G4 motifs (105,914 with G4 and 72,550 without) than GC-rich regions within IZs (8,920 with G4 and 6,378 without; chi-square test, p = 0.078). Therefore, although we find that IZs are GC rich (Figure S7C), as previously reported [START_REF] Xu | Genome-wide identification and characterization of replication origins by deep sequencing[END_REF][START_REF] Cayrou | The chromatin environment shapes DNA replication origin organization and defines origin classes[END_REF], they appear to contain G4 motifs as a consequence of their GC richness, but do not appear to be enriched for G4 motif DNA sequences, per se.

Early Initiation Events in Late-Replicating Domains

Although 91% of our early initiation events are in early-replicating regions of the genome, 9% occur in late-replicating domains (Figures 1F,3B). By visual inspection, the replication tracks in late-replicating parts of the genome appear to be in the locally earliest-replicating parts of the late-replicating regions (for instance at 56 and 58 Mb in Figure 1F), suggesting that they reflect the location of the IZs normally used to replicate these regions. To test this suggestion, we calculated the enrichment of replication tracks around IZs identified as replication timing peaks [START_REF] Chen | Replication-associated mutational asymmetry in the human genome[END_REF] grouped into four quartiles, defined by average replication timing (S50). We find significant enrichment of replication tracks IZs in first three quartiles, (but no enrichment in the fourth quartile, which contains only 3% of the genome-wide replication timing peaks), demonstrating that our ORM signal in late-replicating regions is associated with IZs (Figure 5A). Moreover, the signal-to-noise ratios in the first three quartiles is similar (SNR = 3.2-2.8, Figure 5A), demonstrating that, although there is much less signal in the third quartile, it is enriched over background to the same extent as in the earlier quartiles. In contrast, no enrichment is seen in any quartile in asynchronous cells (Figure 2A). These results confirm that the enrichment of early initiation events we see at late!replicating IZs cannot be attributed to contaminating asynchronous background signals and, therefore, presumably results from bona fide early initiation events in IZs that, on average, replicate late in S phase.

Computational Modeling Suggests that Early-and Late-Replicating Initiation Zones are Distinguished Only by Relative Initiation Probability

The distribution of early initiation events we observe-frequent firing in early-replicating IZs and rare firing in late-replicating IZs-is reminiscent of the stochastic regulation of origin-firing timing observed in budding yeast [START_REF] Czajkowsky | DNA Combing Reveals Intrinsic Temporal Disorder in the Replication of Yeast Chromosome VI[END_REF][START_REF] Yang | Modeling genome-wide replication kinetics reveals a mechanism for regulation of replication timing[END_REF][START_REF] De Moura | Mathematical modelling of whole chromosome replication[END_REF]. Models of stochastic replication timing posit that the only difference between early-and late-firing origins is their relative probability of firing. More!efficient origins are more likely to fire first and thus, on average, have earlier firing times, but sometimes they still fire late, whereas less-efficient origins are less likely to fire early and thus, on average, have later firing times, but do sometimes fire early [START_REF] Rhind | Reconciling stochastic origin firing with defined replication timing[END_REF].

Such models predict that measuring the relative probability of initiation in early-S phase is sufficient to determine the kinetics of initiation for all IZs throughout S phase. If this prediction holds true, we should be able to calculate the complete replication timing profiles of cells, even though we have assayed firing probabilities in only the first 2% of S phase. To test this prediction, we used Replicon, a stochastic replication simulator that calculates replication profiles from the probability of initiation at each point in the genome (Gindin et al., 2014a). We find that we can, indeed, accurately recapitulate genome-wide replication timing profiles with our ORM data (r = 0.85, Figure 5B). Further, our ORM data predicts replication timing better than other initiation data sets, such as, OK-Seq (r = 0.81), Ini-seq (r = 0.69) or SNS-seq (r = 0.48, Figure S6B). ORM predicts replication kinetics equally as well as H3K4me1 (r = 0.85), which marks open enhancer and promoter chromatin, and DNase I hypersensitivity (r = 0.85), which had previously been identified as the best predictor of initiation (Gindin et al., 2014a). These result support the hypothesis that the primary difference between early-and late-replicating IZs is their relative initiation probability, influenced by chromatin accessibility, and not a qualitative difference in their initiation times.

Discussion

We have developed Optical Replication Mapping (ORM), which, for the first time, allows genome-wide, high-throughput, single-molecule replication mapping in human cells. The approach is based on in vivo labeling of replication forks and in vitro detection and genomic mapping of incorporated label on long, single DNA fibers using the Bionano Saphyr platform. Single flow cells using the current mapping approach yielded an average 150-fold coverage of the human genome on fibers from 150 kb to 2 Mb, with an average length of 300 kb (Table S1). Altogether, we have collected 2,550-fold coverage with an average fiber length of 284 kb. Our method is similar to the recently described HOMARD approach, which was used to map DNA replication in Xenopus egg extracts [START_REF] De Carli | High-Throughput Optical Mapping of Replicating DNA[END_REF].

Direct, single molecule, measurements of initiation frequency

Detection of DNA synthesis on individual DNA fibers is the most direct way to localize sites of DNA replication. Since the preponderance of the literature describes low-frequency, heterogeneous selection of replication initiation sites in mammalian cells, single-molecule measurements may be the only means to accurately map the relative frequencies of initiation events in the human genome. However, to date, single-molecule replication-mapping technologies, such as DNA combing and SMARD, which have throughput of a few hundred fibers, have provided insufficient depth of coverage for metazoan genomewide analyses. Recently, nanopore sequencing has been used to map replication incorporation tracks genome-wide on single DNA fibers, but the current technology only allows for relatively small datasets of relatively short fibers [START_REF] Georgieva | Detection of Base Analogs Incorporated During DNA Replication by Nanopore Sequencing[END_REF][START_REF] Müller | Capturing the dynamics of genome replication on individual ultra-long nanopore sequence reads[END_REF][START_REF] Hennion | FORK-seq: replication landscape of the Saccharomyces cerevisiae genome by nanopore sequencing[END_REF]. The largest dataset from a recent study describing the DNAscent approach to nanopore-based replication mapping in budding yeast was 3.8 Gb (the equivalent of just over 1x coverage of the human genome) with an average fiber length of 32 kb [START_REF] Müller | Capturing the dynamics of genome replication on individual ultra-long nanopore sequence reads[END_REF]. Therefore, ORM, for which current typical datasets are 150times larger, provides a unique opportunity for single-molecule, genome-wide mapping of metazoan replication at high coverage on long fibers.

Using ORM, we have mapped the replication initiation events occurring in the first 2% of S phase in aphidicolin-synchronized HeLa cells. These initiation events, as others have reported [START_REF] Petryk | Replication landscape of the human genome[END_REF][START_REF] Long | H2A.Z facilitates licensing and activation of early replication origins[END_REF][START_REF] Pourkarimi | Spatiotemporal coupling and decoupling of gene transcription with DNA replication origins during embryogenesis in C. elegans. Elife 5[END_REF][START_REF] Ganier | Metazoan DNA replication origins[END_REF][START_REF] Cayrou | The chromatin environment shapes DNA replication origin organization and defines origin classes[END_REF], are enriched in the enhancer regions of active genes and are enriched in active chromatin modifications, in particular DNase I HS, H2AZ and H3K4me1 (Figure 4). However, unlike previous reports, we can directly estimate the frequency with which initiation occurs in each IZ because we measure all fibers, both labeled and unlabeled. We find that the distribution of early initiation events is heterogeneous. Although 99% of early IZs predicted by OK-seq analysis are detected in our ORM data, only 6% of them, on average, initiate at the beginning of S phase in any given cell. Furthermore, although most of the initiation events that we observe are in early-replicating IZs, a small but significant minority are in late-replicating IZs, and the frequency of initiation is correlated with the replication timing of the IZ (Figure 3B). This observation is not due to contamination of our synchronized G1/S cells with non-synchronized cells because the signal in asynchronous cells shows no enrichment at specific sites, let alone previously mapped late initiation sites (Figure 2A). These results suggest that, in any given cell, a small number of potential initiation sites are stochastically selected from a large number of potential initiation sites, and that the probability of a site being selected is correlated with the replication time of that site, with early-replicating sites having a higher probability of initiating than late-replicating sites (Figure 3B).

We have mapped replication initiation sites across the genome with over 1,000-fold genome coverage, allowing us to identify initiation sites used in as few as 0.1% of cells. The initiation events that we map appear to be spread over broad IZs, consistent with previous observations [START_REF] Petryk | Replication landscape of the human genome[END_REF][START_REF] Mesner | Bubbleseq analysis of the human genome reveals distinct chromatin-mediated mechanisms for regulating early-and late-firing origins[END_REF][START_REF] Besnard | Unraveling cell type-specific and reprogrammable human replication origin signatures associated with G-quadruplex consensus motifs[END_REF][START_REF] Hamlin | A revisionist replicon model for higher eukaryotic genomes[END_REF][START_REF] Dijkwel | Initiation sites are distributed at intervals in the Chinese hamster dihydrofolate reductase origin of replication but are used with very different efficiencies[END_REF][START_REF] Demczuk | Regulation of DNA Replication within the Immunoglobulin Heavy-Chain Locus During B Cell Commitment[END_REF][START_REF] Lubelsky | Pre-replication complex proteins assemble at regions of low nucleosome occupancy within the Chinese hamster dihydrofolate reductase initiation zone[END_REF][START_REF] Tao | Major DNA replication initiation sites in the c-myc locus in human cells[END_REF][START_REF] Anglana | Dynamics of DNA replication in mammalian somatic cells: nucleotide pool modulates origin choice and interorigin spacing[END_REF]. However, the sparse labeling and relatively long length (~40 kb) of our incorporation tracks allow us estimate the location of initiation sites at a resolution of only about 15 kb. Nonetheless, in 93.5% of IZs, we observe incorporation track centers that are more than 15 kb apart, suggesting that they initiated from distinct sites (Figures 3D andS4). Furthermore, the distribution of signal across IZs is not consistent with a single initiation site per zone (Figure S4C). We cannot distinguish between frequent, discrete, origins scattered throughout each IZ and diffuse IZs in which initiation can occur almost anywhere. However, we can rule out the hypothesis that human replication initiates at isolated, efficient, well-defined replication origins. In particular, we see no evidence for unique initiation sites at the lamin B or Top1 loci (Figure S4A), which have been proposed to be isolated origins [START_REF] Keller | The origin recognition complex marks a replication origin in the human TOP1 gene promoter[END_REF][START_REF] Abdurashidova | Start sites of bidirectional DNA synthesis at the human lamin B2 origin[END_REF]. We suspect that the inefficient nature of human IZs and the low signal-to-noise ratio of previous initiation-mapping technology made it difficult to comprehensively map initiation sites across IZs. However, our data suggests that human IZs contain multiple low-efficiency initiation sites, as has been previously suggested [START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF][START_REF] Petryk | Replication landscape of the human genome[END_REF][START_REF] Mesner | The dihydrofolate reductase origin of replication does not contain any nonredundant genetic elements required for origin activity[END_REF][START_REF] Hamlin | A revisionist replicon model for higher eukaryotic genomes[END_REF].

Human Replication Initiation Events are not Clustered

Our single-molecule data also allows us to directly measure correlations between initiation sites on individual chromosomes. A number of previous reports have claimed that initiations events are clustered, in that initiation events are more likely to occur next to one another on chromosomes than would be expected by chance, and that these clustered events are regularly spaced between 75 and 150 kb apart in human cells [START_REF] Jackson | Replicon clusters are stable units of chromosome structure: evidence that nuclear organization contributes to the efficient activation and propagation of S phase in human cells[END_REF][START_REF] Cayrou | Genome-scale analysis of metazoan replication origins reveals their organization in specific but flexible sites defined by conserved features[END_REF][START_REF] Lebofsky | DNA Replication Origin Interference Increases the Spacing between Initiation Events in Human Cells[END_REF][START_REF] Huberman | Direction of DNA replication in mammalian cells[END_REF][START_REF] Blow | Replication origins in Xenopus egg extract are 5-15 kilobases apart and are activated in clusters that fire at different times[END_REF][START_REF] Marheineke | Control of replication origin density and firing time in Xenopus egg extracts: role of a caffeine-sensitive, ATR-dependent checkpoint[END_REF]). However, these observations have been largely anecdotal because they have recorded observations only on fibers that have multiple initiation events, and so have not rigorously measured their frequency relative to single initiation events, which is necessary to determine if multiple initiation events happen more often than would be expected by chance. To test whether initiation events tend to occur in clusters, we examined the frequency and distribution of multiple tracks on single fibers. If events on a line are clustered with a characteristic inter-event distance, the distribution of inter-event distances will be bell-shaped, with a peak at the most frequent inter-event distance. If, however, the events are distributed randomly, with no inter-event correlation, the distribution of inter-event distances will be exponential [START_REF] Birnbaum | Statistical Methods for Poisson Processes and Exponential Populations[END_REF]. We find that, between 60 and 600 kb, the distribution of inter-initiation-event distances is exponential (r = 0.99, Figure 3E), inconsistent with any significant clustering of initiation.

As an orthogonal test of initiation clustering, we directly examined co-replication of neighboring IZs. We examined every pair of IZs on each fiber we collected and asked whether the frequency with which they were replicated was correlated (Figure 3F). Using a model that takes into account both passive replication predicted by incorporation track lengths (Figures S2B) and nucleotide-uptake heterogeneity (Figure S1C) but otherwise assumes no correlation among initiation events, we can account for the observed frequency of replication in neighboring IZs (Figure 3F), consistent with the conclusion that human replication initiation sites are not clustered in early S phase. Thus, although initiation events with similar replication timing tend to be clustered because regions of similar replication timing are grouped together (Figure 1F), initiation events are no more likely to occur next to one another, or in collinear clusters, than expected by chance (Figures 3E,F).

Stochastic Usage of Initiation Zones in Human Cells

The heterogeneous, uncorrelated nature of replication initiation that we observe is reminiscent of stochastic models that have been proposed to explain the regulation of replication timing in budding yeast. Such models propose that the average timing of initiation in a population is regulated by the probability of a site initiating in an individual cell. In particular, neither hierarchical regulation of initiation timing nor correlation between initiation at neighboring sites is predicted by such models. Previous analyses of human replication kinetics have suggested positive correlations, in which one initiation event, or the forks it produces, would increase the probability of other nearby initiation events [START_REF] Guilbaud | Evidence for Sequential and Increasing Activation of Replication Origins along Replication Timing Gradients in the Human Genome[END_REF][START_REF] Löb | 3D replicon distributions arise from stochastic initiation and domino-like DNA replication progression[END_REF]. However, a more recent, higher resolution replication timing study found no evidence of this so-called "domino model" [START_REF] Zhao | High-resolution Repli-Seq defines the temporal choreography of initiation, elongation and termination of replication in mammalian cells[END_REF]. In any case, we find that a stochastic model assuming independent origin initiation can explain our observed data (Figure 3F).

To test whether our results are compatible with a stochastic model of replication kinetics, we used a computational model of DNA replication that assumes stochastic origin firing (Gindin et al., 2014a). Using our initiation-site mapping data to initialize the model produces predicted replication timing profiles that give an excellent fit (r = 0.85) to experimentally produced profiles (Figure 5B). Since the computational model does not contain any correlation between initiation events, this analysis suggests that any positive or negative correlations between initiation events must be of sufficiently low magnitude or restricted spatial distribution as to be not required to explain observed genome-wide patterns of replication timing. Our work, in combination with work from budding and fission yeast [START_REF] Kaykov | The spatial and temporal organization of origin firing during the Sphase of fission yeast[END_REF][START_REF] Yang | Modeling genome-wide replication kinetics reveals a mechanism for regulation of replication timing[END_REF][START_REF] Patel | DNA replication origins fire stochastically in fission yeast[END_REF][START_REF] De Moura | Mathematical modelling of whole chromosome replication[END_REF], supports the hypothesis that stochastic control of initiation timing is a fundamental feature of DNA replication, conserved across eukaryotes.

General Application of ORM to Map Replication Kinetics

In addition to mapping replication initiation sites, we also used ORM to map ongoing replication forks in unperturbed, asynchronous cells. We inferred fork direction from the asymmetric nature of incorporation tracks caused by the reduction of incorporation as labeled nucleotides are depleted, an approach that has been used before for both radio-and fluoro-labeled nucleotides [START_REF] Müller | Capturing the dynamics of genome replication on individual ultra-long nanopore sequence reads[END_REF][START_REF] Huberman | On the mechanism of DNA replication in mammalian chromosomes[END_REF][START_REF] Hennion | FORK-seq: replication landscape of the Saccharomyces cerevisiae genome by nanopore sequencing[END_REF]. Averaged over our nearly-500-fold fiber coverage of the genome, our inferred fork directions allow us to infer genome-wide replication kinetics that agree with independently determined replication-timing and replication-fork-directionality profiles (Figures 2C,D). In future work, it should be possible to use sequential labeling with nucleotides of different colors to unambiguously determine fork direction, thereby greatly increasing the resolution and sensitivity of asynchronous ORM. In addition to identifying replication initiation sites throughout the S phase of any cell type, sequential labeling will allow investigation of other important aspects of DNA replication, recombination and repair. We envision ORM being used to measure replication fork speed, replication fork arrests and reversal, and sister-chromatid exchange, under both normal growth and replication stress conditions. Such application will make ORM a central technique for studying DNA replication, DNA repair and genome instability.

Limitations

The primary limitation of current ORM technology is the in vivo pulse-labeling of replication. Other single-fiber DNA-replication-mapping approaches, such as DNA combing [START_REF] Bensimon | Alignment and sensitive detection of DNA by a moving interface[END_REF] and SMARD [START_REF] Norio | Visualization of DNA replication on individual Epstein-Barr virus episomes[END_REF], use thymidine analogs to label replicated DNA, allowing kilobase-resolution mapping of replication kinetics in a wide range of cell types and organisms. However, neither the two standard thymidine analogs-BrdU and EdU-are compatible with ORM technology; BrdU because anti-BrdU antibodies only recognize single-stranded DNA, but double-strand DNA is required for fiber stretching in the Bionano optical mapping chips, and EdU because the copper catalyst required to conjugate the EdU to a fluorophore produces hydroxyl radicals, which breaks the DNA fiber into sub-hundred-kb pieces. We have therefore labeled our cells by electroporating them with fluorescently-conjugated dUTP, which limits us to using cell lines. Moreover, the conjugated nucleotides are not efficiently incorporated leading to the sparse labeling and limited resolution that we see. DNAreplication labeling with a number of other nucleoside analogs has been demonstrated recently [START_REF] Neef | An azide-modified nucleoside for metabolic labeling of DNA[END_REF][START_REF] Rieder | Alkene-Tetrazine Ligation for Imaging Cellular DNA[END_REF][START_REF] Mitter | Conformation of sister chromatids in the replicated human genome[END_REF]. Although none of them are currently compatible with ORM, we predict that improved nucleoside-labeling technology will allow ORM to by applied to most cells, tissues and organisms at kilobase resolution, making ORM a generally applicable technique for replication mapping.

Method Details

Cell Labeling and Sample Preparation

For asynchronous HeLa datasets, cells were grown to ~70% confluency, washed with warm PBS, trypsinized, and electroporated (Lonza, Nucleofection kit SE, HeLa S3 HV program) in the presence of 40 µM Aminoallyl-dUTP-ATTO-647N (Jena Bioscience NU-803-647N). For G1/S synchronized data sets, cells were synchronized at metaphase by incubation with 0.05 µg/ml nocodazole for 4 hours followed by shake off. The percentage metaphase cells was quantified by metaphase spread. Briefly, 3X volumes of ddH20 was added to an aliquot of media from shake off and was incubated at 37˚C for 15 minutes. Cells were fixed by adding several drops of ice-cold methanol:glacial acetic acid (3:1). Cells were spun down and pellet were resuspended in ice-cold methanol:glacial acetic acid (3:1) and spotted on to a clean slide. Slides were allowed to dry, spotted with Vectasheild (Vector Laboratories H-1000) containing DAPI and covered with coverslips. Metaphase and interphase cells were counted on a Nikon Eclipse Ti microscope. Metaphase cells were released into fresh media containing 10 µg/ml aphidicolin (Millipore Sigma 178273) for 16 hours. For 0-minute release datasets, cells were washed three times with warm PBS (or PBS plus aphidicolin, in the case of B.0) to remove aphidicolin, trypsinized, and electroporated (Lonza, Nucleofection kit SE, HeLa S3 HV program) in the presence of 40 µM Aminoallyl-dUTP-ATTO-647N (Jena Bioscience NU-803-647N). Cells were recovered in warm media plus aphidicolin for 30 minutes in a 37˚C water bath with agitation every 5 minutes to prevent cells from settling in the tube. Cells were then spun down and washed 3 times with warm PBS to remove aphidicolin. For timed release experiments, cells were washed three times with warm PBS to remove aphidicolin, warm media was added, and cells were returned to 37˚C/5% CO2 incubator for the indicated times (5, 10, 20, 30, 45, 60, or 90 minutes). Cells were then washed again with PBS, trypsinized, and electroporated in the presence of fluorescent nucleotides as above. In all cases, cells were allowed to recovered overnight in media at 37˚C/5% CO2. After recovery, all cells were detached with trypsin and snap frozen in liquid nitrogen. H9 cells were grown to ~70% confluency, washed with warm PBS, dissociated to single cells with Gentle Cell Dissociation Reagent (StemCell Technologies 07174) and electroporated (Lonza, Nucleofection kit P3, hESC program) in the presence of 40 μM Aminoallyl-dUTP-ATTO-647N (Jena Bioscience NU-803-647N). Cells were recovered overnight in StemPro media plus 10uM Y-27632 dihydrochloride (ROCK inhibitor, Sigma Y0503) at 37˚C/5% CO2.

Nucleotide labeling efficiency was assayed by plating a small aliquot of cells on coverslips after electroporation. Following overnight recovery, cells were fixed with 4% formaldehyde for 10 minutes at room temperature. Cells were then washed twice with PBS followed by permeabilization and counterstaining with PBS plus 0.2% Triton X-100 plus DAPI for 5 minutes at room temperature. Coverslips were washed twice with PBS, spotted with Vectasheild, and placed on slides. Images were captured on a DeltaVision (GE Life Sciences) microscope. Replication foci patterns were categorized as previously described [START_REF] O'keefe | Dynamic organization of DNA replication in mammalian cell nuclei: spatially and temporally defined replication of chromosome-specific alphasatellite DNA sequences[END_REF].

Optical Replication Mapping

Two Bionano mapping approaches were used: Nick, Label, Repair and Stains (NLRS) and Direct Label and Stain (DLS) (Table S1). For NLRS, cell pellets were thawed on ice, resuspended in cold Cell Buffer (Bionano 80004) and embedded into low-melting point agarose plugs using the CHEF Mammalian Genomic DNA Plug Kit (Bio-Rad #1703591) according to the manufacturer's instructions, with a target concentration of 1,000,000 cells per plug. Cells immobilized in plugs were lysed and treated with Proteinase K and RNase A. The plugs were then melted and solubilized, and the resulting DNA was further cleaned by drop dialysis. The DNA was fluorescently labeled with green fluorophores at Nt.BspQI sites and counterstained with YOYO-1 using Bionano Genomics' NLRS labeling kit (Bionano 80001). Throughout the entire process, the samples were shielded from light whenever possible, in order to maintain integrity of the fluorophores. For DLS, DNA was labeled with green fluorophores using the direct labeling enzyme DLE-1 (Bionano Genomics 80005) according to manufacture's protocol. The rest of processing steps were the same as NLRS labeling.

Labeled DNA samples were loaded onto a Saphyr chip (Bionano Genomics 20319) and run on the corresponding Saphyr instrument (Bionano Genomics IN-011-01), following the manufacturer's instructions. Briefly, the DNA was coaxed into an array of parallel nanochannels via electrophoresis, thereby elongating the DNA to a uniform contour length, allowing for accurate measurement of distances along the molecules. Molecules were imaged to collect the YOYO-1 DNA signal in a false-color blue channel, the Nt.BspQI or DLE-1 labels in the green channel, and the in vivo replication incorporation tracks in the red channel.

Images were converted to digitized molecules and the positions of the green and red labels on each fiber were determined using Bionano Access. The Nt.BspQI and DLE-1 labels were aligned to human reference genome hg19 by creating an in silico digestion of the reference FASTA file. Alignment of the fibers' Nt.BspQI and DLE-1 sites allowed mapping of the red replications signal to the genome. The pipeline produces one .bnx file containing the intensity and location of each green and red signal on each fiber and one .xmap file containing the genomic location of each fiber. We observed abnormal ORM signal enrichment in some specific genomic positions (hereafter called hotspots), which only occurred in the samples mapped by NLRS technique but not those mapped by DLS technique. Further investigation confirmed that most (99.36%) these hotspots are located within < 650 bp distance to a BspQI site used for the mapping in the NLRS technique. Detailed analysis showed that these abnormal ORM signals are likely signals of false positive detection due to the strong BspQI signals, since the green signals (for mapping) associated with the hotspot regions have stronger intensity distribution than the green signal outside the hotspot region while the red signals inside the hotspot regions have intensity weaker than the red signals outside of the hotspots. In order to filter these false positive ORM signals, we calculated the ORM signal count (mean = 1.71 and s.d. = 6.35 for all windows containing at least one red signal) within all 300 bp replication fork. The mean relative fork direction (RFD) within a given window is defined as (R-L)/(R+L), where R and L is the number of rightward (FDI>0) and leftward (FDI<0) ORM tracks, respectively.

Identification of Initiation Zones

Initiation zones were identified based on the ORM signal density along the genome. We first obtained the ORM signal density by calculating the percentage of fibers containing ORM signal in 10 kb sliding windows with 1 kb steps. LOESS smoothing (with α = 0.75 and the polynomial degree = 2) was then performed in successive 160 kb windows, with 8 kb overlaps to avoid discontinuity of fitting. Within the 8 kb overlaps, the final smoothed values were the averages, weighted by their distance to the corresponding window, of the smoothed values from two adjacent windows. IZs were then defined as the peaks on the smoothed ORM signal density profile. The size of each IZs was defined by the smallest window that containing at least 40% of ORM signals within each peak. The value of 40% was chosen because, when aggregated around IZ peaks, the ORM signal density decreases with distance from the peak at a similar rate until about 40% then decreases much more slowly. This process was performed on all 4 biological replicates of HeLa 0' data as well as the combined dataset. Only IZs identified in at least 3 replicates and with a relative peak height in the combined dataset (as measured by the difference in signal density between the peak and the edges of the zone) of greater than 0.3% were retained. The precise location of each retained IZ was then determined by analyzing the transition of ORM signals with a K-mean clustering (n=2 if the center cluster < 30 kb, otherwise n=3) on the ∆ORM signal between adjacent 1kb windows to identify the region with sharpest ORM density around each IZ center.

To estimate the minimum number of initiation sites in each IZ, we calculated number of replication track centers separated by more than the resolution of our ability to localize initiation sites, which is about 15 kb. We calculated the position of the center of each replication track that overlaps a given IZ, discarding any that are more than 7.5 kb outside of it. We then ranked the centers by position and, starting at one end, counted the number of centers that are more than 15 kb apart.

Comparative Analysis of ORM Data with Other Datasets

The raw Repli-Seq data of HeLa S3 cells were downloaded from the Encode project <http://genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&g=wgEncodeUwRepliSeq> [START_REF] Hansen | Sequencing newly replicated DNA reveals widespread plasticity in human replication timing[END_REF] and S50 (#he fraction of S phase at which 50% of the DNA is replicated in a defined genome region) was computed [START_REF] Chen | Impact of replication timing on non-CpG and CpG substitution rates in mammalian genomes[END_REF] and used to identify replication timing peaks [START_REF] Chen | Replication-associated mutational asymmetry in the human genome[END_REF] as previous described. Published H9 ESCs replication timing and peak calling data were used [START_REF] Zhao | High-resolution Repli-Seq defines the temporal choreography of initiation, elongation and termination of replication in mammalian cells[END_REF]. ORM signal enrichment around timing peaks was calculated in R (version 3.5.1 <https://www.rproject.org>). To compare the ORM IZs and other origin mapping data, we used pROC package to calculate the ROC (Receiver Operating Characteristic) curve and AUC (Area Under the Curve) values by using all 100 kb bins along the genome as example space and the bins overlapped with the ORM IZs as true positive space [START_REF] Robin | pROC: an open-source package for R and S+ to analyze and compare ROC curves[END_REF].

Correlations between ORM initiation efficiency and published genome-wide HeLa cell replicationinitiation-mapping datasets-OK-seq [START_REF] Petryk | Replication landscape of the human genome[END_REF], SNS-seq [START_REF] Picard | The spatiotemporal program of DNA replication is associated with specific combinations of chromatin marks in human cells[END_REF], Ini-seq [START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF] and Orc1 ChIP-seq [START_REF] Dellino | Genome-wide mapping of human DNAreplication origins: levels of transcription at ORC1 sites regulate origin selection and replication timing[END_REF]-were calculated across the genome in 50 kb adjacent windows. When necessary, genomic coordinates were remapped to hg19 using LiftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver). ORM initiation efficiency was calculated as the number of fibers with ORM signal in a 50 kb genomic window divided by the number of fibers covering that window. The OK-seq initiation efficiency was calculated either genome-wide using the previously described origin efficiency metric (OEM) approach [START_REF] Mcguffee | Quantitative, genome-wide analysis of eukaryotic replication initiation and termination[END_REF] or for each OK-Seq IZ with the DRFD between its right and left extremities [START_REF] Petryk | Replication landscape of the human genome[END_REF]. Correlations were calculated using Spearman's correlation coefficient in R.

The genome-segmentation results (ChromHMM) based on ENCODE data were retrieved from UCSC genome browser <https://genome.ucsc.edu/cgibin/hgFileUi?db=hg19&g=wgEncodeAwgSegmentation>. To compensate for the large size of IZ, relative to many of the ChromHMM segments, we used the IZ's core regions, defined as the region from center of the IZ to the point of highest ORM signal within that IZ.

Histone modification data for HeLa S3-H2A.Z, H3K4me1, H3K4me2, H3K4me3, H3K9ac, H3K9me3, H3K27ac, H3K27me3, H3K36me3, H3K79me2 and H4K20me1-was downloaded from ENCODE <https://www.encodeproject.org>. If available, the data was obtained from G1 phase cells. The average of all replicates was use in our analysis. The RNA Pol II-pS2 <https://www.encodeproject.org/experiments/ENCSR000ECT> and DNase I hypersensitivity data (UW DNase I HS data track) for HeLa S3 were also downloaded from ENCODE. The signal enrichment of each feature was re-normalized such that the background regions show zero log2 enrichment ratio. Published G4 sequence locations were retrieved (Puig [START_REF] Lombardi | Thermodynamically stable and genetically unstable G-quadruplexes are depleted in genomes across species[END_REF]. We analyzed 3 G4 datasets (G4L3, G4L7 and G4L12) with similar results; we therefore used the G4L12 (G3N1-12G3N1-12G3N1-12G3), which contains the largest dataset and give the best correlation with ORM data, in our final analysis. The genome-wide pairwise correlation analysis and enrichment around ORM IZs were performed by deepTools [START_REF] Ramírez | deepTools2: a next generation web server for deep-sequencing data analysis[END_REF] or with R (version 3.5.1 <https://www.r-project.org>). IZs were classified by timing-early (S50<0.25), mid-early (0.25 ~ 0.5), mid-late (0.5~0.75) and late (S50>0.75)using the S50 of the IZ center. The partial correlation analysis between ORM initiation efficiency and G4 sequences with controlling the effect of GC percentage (obtained from UCSC <http://hgdownload.soe.ucsc.edu/goldenPath/hg19/gc5Base>) was performed by using the R package ppcor <https://cran.r-project.org/web/packages/ppcor/index.html>. Features were visualized in IGV [START_REF] Robinson | Integrative genomics viewer[END_REF] and the ORM Browser <http://orm.nucleome.org> (Figure S3, Zhu et al. in prep, personal communication, J. Ma <https://github.com/nimezhu/ORM>). All published data sources are complied in Table S2.

Correlation of Initiation in Neighboring Initiation Zones

To calculate the degree of correlation between initiation events in neighboring IZs, we defined an initiation event as any ORM signal in any IZ on a fiber. We then calculated the correlation as the frequency with which we observed initiation in neighboring IZs on the same fiber divided by the product of the frequency of overall initiation in the two IZs. The non-neighbor correlation was calculated in the same manner, but using all pairwise combinations of IZs on each fiber farther that 200 kb apart. To calculate the expected correlation due to passive replication of one IZ due to initiation in a neighboring IZ, we used the observed label-depletion rate and inferred the frequency with which initiation at one IZ would lead to incorporation of label at a neighboring IZ at a given distance away (Figure S2B). To calculate the correlation expected from heterogeneity of label uptake by cells, we used the observed distribution of cellular labeling (Figure S1C). The two calculations were combined to produce the model (Eq. 34, Supplemental Mathematical Methods) shown in Figure 3F. See Supplemental Mathematical Methods for details.

Replication Timing Simulation

To simulate the replication timing profiles, we used the Replicon simulation code (Gindin et al., 2014a). Replicon uses three sets of user-defined parameters. The first set of parameters is the initiation probability landscape (IPLS), the relative probability of initiating at any point in the genome. In our simulations, we used the ORM the signal distribution or other genomic features with a bin size of 0.5 kb. The second set of parameters control the virtual flow sorter. In agreement with the results from the original paper (Gindin et al., 2014b), we set this to (0.0, 0.17, 0.35, 0.58, 0.92, 0.99, 1.0), although other choices give similar results. The final parameter is the number of forks. Here, we again followed the original paper and set the number of forks equal to 10.24 + x(7.9E-7), where x is the length of the chromosome in base pairs.

Figure 1. Optical Replication Mapping of Synchronized Cells Identifies Early Initiation Events

HeLa S3 cells were synchronized by mitotic shake off, arrested in early S-phase with aphidicolin, electroporated with fluorescent dUTP and released from the aphidicolin arrest, allowing fluorescent nucleotide incorporation around early replication initiation events. A) A representative image of the Bionano data with DNA in blue, Nt.BspQI restriction sites in green and incorporated nucleotides in red. The field of view is 750 kb. B) Enrichment of ORM signals from the combined 0-minute dataset around early replication timing peaks, those that replicate in the first quarter of S phase. C) Enrichment around early replication timing peaks of ORM signals from datasets in which cells were labeled between 5 and 90 minutes after release from synchronization. D, E) Distribution of replication tracks at two previously characterized human replication IZs: Myc (C) and Top1(D). The distribution of replication tracks in shown beneath the HeLa S3 replication timing profile [START_REF] Hansen | Sequencing newly replicated DNA reveals widespread plasticity in human replication timing[END_REF] and genes from the region. The distribution of total fibers in shown in grey at 1% the scale of the replication tracks. Discontinuities in the fiber density distribution, such as the one around megabase 127.5 in D, are probably caused by segmental duplication increasing the copy number of the region (see Figure S1E). Therefore, all analyses are normalized to fiber coverage depth. B) The density of ORM signals from the HeLa asynchronous dataset around replication timing peaks separated into replication timing quartiles, defined by average replication timing (S50, see Methods). C) Labeled ORM replication forks, with leftward-moving forks colored blue and rightward-moving forks colored red, are depicted below the replication timing profile of the Top1 region and its OK-seq relative fork direction (RFD) data, with blue dots indicating leftward replication, red dots rightward, and the vertical position of the dot, between -1 and 1, indicating the fraction of replication moving in that direction. The polarity of ORM replication forks was determined by calculating their fork direction index (FDI, see Methods). D) The RFD calculated from ORM and OK-seq data across 7 Mb of Chromosome 1. E) The enrichment of HeLa fork polarity, as measured by RFD, around early replication timing peaks. F) The enrichment of HeLa and H9 fork polarity around cell-type specific and cell-type general early replication timing peaks.

Figure 3. Mapping of Human Replication Initiation Zones

A) A comparison of ORM tracks, ORM IZs, OK-seq RFD (as in Figure 2C), OK-seq IZs [START_REF] Petryk | Replication landscape of the human genome[END_REF], SNS-seq [START_REF] Picard | The spatiotemporal program of DNA replication is associated with specific combinations of chromatin marks in human cells[END_REF], Ini-seq [START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF] and Orc1 ChIP-seq [START_REF] Dellino | Genome-wide mapping of human DNAreplication origins: levels of transcription at ORC1 sites regulate origin selection and replication timing[END_REF] at the Top1 locus. B) The correlation between ORM initiation efficiency (fraction of labeled fibers/IZ) in the combined 0minute dataset and replication timing (S50) in the 4,930 ORM IZs. Spearman rank correlation coefficient between the two datasets is 0.75. Inset is the same data on a semi-log plot fit with an exponential curve. C) A Venn diagram of the overlap between OK-seq early IZs and 50 kb windows in the genome that have ORM initiation efficiencies of at least 3% or 10%.

D)

The estimated minimal number of initiation sites per IZ. The estimate was made by calculating minimum number of replication tracks in each IZ whose centers are more than 15 kb apart. Most of the IZs for which all track centers are within 15 kb of each other are small and contain few replication tracks (Figure S4B). E) The distribution of the inter-initiation distances in the combined 0' dataset measured as the distance from the middle of neighboring replications tracks for fibers that have multiple tracks. The average of the distribution is 111 kb and the mode of distribution is 57 kb. Inset: The distribution, plotted on a log y-axis from 60 to 600 kb, fit to an exponential curve (r = 0.99). The exponential distribution of the interreplication-track distances indicates that the distribution of initiation events on this length scale is random. F) The correlation between the probability of replication at neighboring IZs. The data is well fit by a model that assumes no correlation between IZs (see Supplemental Mathematical Methods). We also calculated the correlation coefficient among all IZs on each fiber, not just the neighboring ones, and found a coefficient similar to that of neighboring IZs separated by more than 200 kb, consistent with the interpretation that the low level of residual correlation seen at longer distances is a global effect due to the heterogeneity in cellular label uptake. G) A matrix depicting the correlation between initiation sites predicted by ORM, replication timing [START_REF] Chen | Impact of replication timing on non-CpG and CpG substitution rates in mammalian genomes[END_REF], OK-seq [START_REF] Petryk | Replication landscape of the human genome[END_REF], SNS-seq [START_REF] Picard | The spatiotemporal program of DNA replication is associated with specific combinations of chromatin marks in human cells[END_REF], Ini-seq [START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF] Orc1 ChIP-seq [START_REF] Dellino | Genome-wide mapping of human DNAreplication origins: levels of transcription at ORC1 sites regulate origin selection and replication timing[END_REF], DNase I hypersensitivity [START_REF] Bernstein | An integrated encyclopedia of DNA elements in the human genome[END_REF]G4 motifs (Puig Lombardi et al., 2019).

Figure 4. Characterization of Human Initiation Zones

Enrichment of histone modifications and other genomic features relative to ORM IZs, separated into replication timing quartiles. The upper panels show the genome-normalized relative signal around all IZs. The lower panels show heat maps of the same signal at each IZ. The y-axis scale for the upper panels and the heat-map scale for the lower panels are the same genome-normalized enrichment in relative units.

Figure 5. Probability of Early Initiation Predicts Replication Kinetics Genome-Wide

A) Enrichment of ORM signals from the combined 0-minute dataset around replication timing peaks separated into replication timing quartiles. Signal-to-noise ratio (SNR) was calculated as maximum signal over average signal between 400 and 500 kb away. B) Comparison between experimentally determined replication timing (S50) and replication timing predicted from ORM data using a stochastic model (r = 0.85, Gindin et al., 2014a).

Figure S1. Characterization of ORM Labeling

A) Micrographs of representative HeLa cells electroporated with ATTO-647-dUTP during an aphidicolin arrest, released, allowed to recover overnight, and fixed. The left panel is stained with DAPI, the middle panel visualizes the incorporated fluorescent nucleotide and the right panel is a merger of the two channels. B) Flow cytometry analysis of S-phase progression after ATTO-647-dUTP electroporation. C) Flow cytometry analysis of ATTO-647-dUTP uptake. Cells arrested in aphidicolin at the beginning of S phase were electroporated with ATTO-647-dUTP, or mock electroporated, incubated on ice and analyzed by flow cytometry. The distribution of labeled cells was fit with a log-normal distribution with a mean of 1141 ± 7.7and a coefficient of variation of 0.88 ± 0.01. D) The distribution of fiber lengths of 0-minute dataset B.0, which was collected with the DLS mapping approach and has slightly longer fibers that datasets collected with the NLRS approach (Table S1). The average fiber length is 338 kb and the fiber N50 is 369 kb.

E) The correlation of labeling between the three biological replicates that make up the C.0 0-minute dataset. The number of signals in each 10 kb bin across the genome is plotted and the correlation coefficient is reported. F) The correlation of labeling between the four biological replicates (A.0, B.0, C.0, D.0) of the 0-minute dataset and one asynchronous dataset, as in panel D. The correlations between the biological replicates are higher than those between the technical replicates because the biological replicates are larger, reducing the counting noise in infrequently labeled regions of the genome. G) The depth of genome coverage in the combined 0-minute dataset in 1 kb bins. The aneuploid character of the HeLa genome is evident in the distribution of coverage into four peaks corresponding to the haploid, diploid, triploid and tetraploid regions of the genome. The coverage data was fit with four Gaussian curves with coverage maxima at about 376x (haploid), 705x (diploid), 1047x (triploid) and 1400x (tetraploid, which was fixed at 1400x, because the unconstrained fit had a very large variation). The individual Gaussians and the complete fit are shown. The coefficients of variation of the individual Gaussians and the percent of the genome inferred to have that ploidy are shown in the legend.

Figure S2. Characterization of ORM Label Distribution

A) The distribution of the intensities of the incorporated fluorescence signals in the D.0 0-minute dataset. The fit (obtained from the first four terms of Eq. 5, Supplemental Mathematical Methods) predicts that about 80% of observed signals are single fluorophores and that the other 20% are multiple fluorophores sufficiently close together that they are not resolved by the Saphyr optics. This estimate of 80% single fluorophores is consistent with an average inter-signal distance of 4 kb and 1.3 kb resolution of the Saphyr, both of which parameters can be inferred from the distribution of inter-signal distances (see panel B and Methods). The distribution of intensities in the other datasets are similar, although differences in the Saphyr optical calibration on different runs introduces variation into the absolute value of the measured fluorescent intensities. B) The distribution of inter-signal distances in the combined 0-minute dataset. The fit to the data (Eq. 14, Supplemental Mathematical Methods) between 10 and 500 kb predicts an initial labeling frequency of 1 in every 877±17 thymidines and a depletion half length of 74.5±0.7 kb. Similar fits for the asynchronous HeLa and H9 datasets predict labeling frequencies of 1/1025 and 1/850 and depletion half lengths of 57 and 48 kb, respectively. The similar labeling densities suggest the nucleotide uptake is similar in all three experiments, whereas the shorter depletion half length is consistent with previous reports that the number of forks increases during S phase, which would consume nucleotides more quickly [START_REF] Yang | How Xenopus laevis embryos replicate reliably: Investigating the random-completion problem[END_REF][START_REF] Goldar | Universal temporal profile of replication origin activation in eukaryotes[END_REF]. C) The enrichment of ORM signals in solo-signal replication tracks from the combined 0-minute dataset around early replication-timing peaks, those that replicate in the first quarter of S phase. D) The distribution of inter-replication-track lengths. E) The distribution of IZ lengths. F) The fraction of 50 kb genomic bin with initiation efficiency greater than indicated on the x axis. 50% of bins have an initiation efficiency greater that 0.3% and 10% of bins have an initiation efficiency of greater that 4.3%.

Figure S3. The ORM Genome Browser

The ORM Genome Browser allows interactive visualization the HeLa ORM data. Shown is a screen shot of the Fibers track of synchronous data. It shows fibers as gray bars, ORM signals as yellow hash marks, and inferred replication tracks as black lines. Only fibers labeled with ORM signals are displayed because only ~5% of fibers are labeled; displaying all fibers is impractical. The browser can also display only the replication tracks, to provide a more easily-visualized view of replication initiation. It can also show the fibers and the replication tracks from the asynchronous data. Two low-efficiency IZs are shown, 2353 (3%) and 2354 (2%), because higher-efficiency IZs contain too many labeled fibers to fit on one screen. Note that, although the signals and replication tracks are concentrated around the IZs, many lay outside of them and none are concentrated in discrete areas.

Figure S4. Distribution of Replication Tracks within Initiation Zones

A) The distribution of replication tracks in the merged 0 minute dataset at the Top1 locus. The Top1 IZ has an estimated minimum of five initiation sites because the five replication track centers indicated in red are all 15 kb away from each other. B) The distribution of replication tracks at four examples of IZs for which our estimate of the minimum number of initiation sites is 1. C) The distribution of signal across IZs. The ratio of signal frequency at the IZ center to the IZ boundary is plotted versus IZ length. This value is expected to decrease more quickly in IZs that predominantly have a single initiation site (Eq. 25) than if initiation is distributed across the IZ (Eq. 28, Supplemental Mathematical Methods). The distribution across IZs shorter than 55 kb is consistent with a uniform distribution of initiation sites. At longer lengths, we actually see more signal at the edges of the IZs than it the center. One possible explanation for this phenomenon is that larger IZs may actually be two smaller IZs fused together such that there is more initiation sites towards the edges of the fused IZ and less initiation in the center, which is actually between the two constituent IZs. See Supplemental Mathematical Methods for details.

Figure S5. Reanalysis of Potentially Discrete OK-Seq Initiation Zones

We reexamined the 66 OK-seq IZs that were reported to be less that 5 kb wide [START_REF] Petryk | Replication landscape of the human genome[END_REF]. A,B) 53 are in regions of noisy OK-seq data. Of those, 24 are in late-replicating regions and appear to be in regions with extensive bi-direction replication. Panel A is an example of one such zone. However, since there is little ORM data in these regions, we can say little more about them. 29 are in early replicating regions, but none of them correlate with numerous ORM segments. Panel B is an example of one such zone. We conclude that they are not active IZs in our ORM data and probably not active IZs in the OK-seq data, either. C) 6 are robust transitions that correlate with numerous ORM segments. We conclude that they are IZs in both the OK-seq and ORM data. However, they show broadly dispersed ORM segments, therefore we do not believe they are unusually constrained IZs. Instead, we conclude that they are outliers in the OK-seq data that were identified as unusually narrow due to experimental variation. Panel C is an example of one such zone. D,E) 7 are robust transitions that do not correlate with numerous ORM segments. They could be IZs present in the OK-seq HeLa cell line, but absent in the ORM HeLa cell line. Alternatively, they could be translocation break points in the OK-seq HeLa cell line relative to the hg19 reference sequence. Such breakpoints would explain both the sharpness of the transition and the absence of these putative IZs from the ORM data. Panel D is an example of one such zone. Panel E show two such zones that can be explained by an inversion between them.

Figure S6. Correlations between ORM Data and Other Datasets

A) ROC analysis of the association between the ORM IZs and other origin mapping data shown in Figure 3G. As shown in the ROC curves, the ORM IZs are better correlated with OK-seq IZs (AUC=0.72) and Ini-seq IZs (AUC=0.73). The bias of the SNS ROC curve towards high true positive rates only at high false negative values is consistent with that dataset having more false positive signal, whereas the bias of the replication timing and ORC datasets relatively high true positive rates only at low false negative values is consistent with those datasets having fewer, but more accurate true positives. Areas under the ROC curves (AUC) are shown in the legend.

B)

Comparison between experimentally determined HeLa replication timing (S50) and replication timing predicted from ORM, DNase I hypersensitivity [START_REF] Bernstein | An integrated encyclopedia of DNA elements in the human genome[END_REF], OK-Seq [START_REF] Petryk | Replication landscape of the human genome[END_REF], Iniseq [START_REF] Langley | Genome-wide identification and characterisation of human DNA replication origins by initiation site sequencing (ini-seq)[END_REF] and SNS-seq [START_REF] Picard | The spatiotemporal program of DNA replication is associated with specific combinations of chromatin marks in human cells[END_REF] data using a stochastic model (Gindin et al., 2014a). The Spearman correlation coefficients with replication timing are shown in the legend.

Figure S7. Enrichment of Histone Modifications and Other Genomic Features in Initiation Zones

A) The enrichment of chromatin regions defined by ChromHMM in IZs relative to the genome in general [START_REF] Ernst | ChromHMM: automating chromatin-state discovery and characterization[END_REF]. The faction of IZs or genomic sequence with each of the indicated annotation is shown for all IZs and the IZs that replicate in the first quarter of S phase (S50 < 0.25). B) The enrichment of histone modifications and other genomic features relative to ORM IZs. The upper panels show the genome-normalized relative ChIP-seq signal around all IZs. The lower panels show heat maps of the same signal at each IZ. The left panel includes enhancer-enrich features; the right panel shows transcription-enriched features.

C) The enrichment of GC content relative to ORM IZs. The upper panels show the average % GC content signal around all IZs. The lower panels show heat maps of the % GC content at each IZ. D) Correlation heat maps at various resolution. The left panel shows 100 kb resolution, which does not resolve enhancers, promoters and transcription units. Therefore, features associated with all three correlate with ORM signal. The center panel shows 1 kb resolution, which resolves enhancers, promoters and transcription units. However, the correlation is dominated by replication timing, creating a correlation between ORM IZs and transcription units, which both tend to replicate early. The right panel shows 1 kb resolution for the earliest-replicating quarter of the genome. Here, enhancer-enriched features, such as H3K4me1, H3K9ac and H3K27ac hypersensitivity, are most strongly correlated, while promoter-enriched features, such as RNA Pol II, H3K4me3, and are more weakly correlated and elongation-enriched features, such as H4K20me1, H3K79me2 and H3K36me3, are anti-correlated.

F)

  Distribution of replication tracks across the right arm of Chromosome 20, as in D. G) Distribution of replication timing of replication tracks compared to replication timing of each 50 kb bin in the genome.

Figure 2 .

 2 Figure 2. Optical Replication Mapping of Asynchronous Cells Identifies Cell-Type-Specific Genome-Wide Replication KineticsA) Asynchronous HeLa S3 and H9 human embryonic stem cells were electroporated with fluorescent dUTP, allowing fluorescent nucleotide incorporation into ongoing replication fork. B) The density of ORM signals from the HeLa asynchronous dataset around replication timing peaks separated into replication timing quartiles, defined by average replication timing (S50, see Methods). C) Labeled ORM replication forks, with leftward-moving forks colored blue and rightward-moving forks colored red, are depicted below the replication timing profile of the Top1 region and its OK-seq relative fork direction (RFD) data, with blue dots indicating leftward replication, red dots rightward, and the vertical position of the dot, between -1 and 1, indicating the fraction of replication moving in that direction. The polarity of ORM replication forks was determined by calculating their fork direction index (FDI, see Methods). D) The RFD calculated from ORM and OK-seq data across 7 Mb of Chromosome 1. E) The enrichment of HeLa fork polarity, as measured by RFD, around early replication timing peaks. F) The enrichment of HeLa and H9 fork polarity around cell-type specific and cell-type general early replication timing peaks.

Acknowledgements

We are grateful to Feng Yue for his contribution to the collection of the Bionano data, Jian Ma for help developing ORM browser and Peiyao Zhao maintaining the ORM browser server. The work was funded by NIH grants HG010658 to DMG and GM125872 to NR. CLC was supported by the I. Curie YPI program, the ATIP-Avenir program from CNRS, Plan Cancer from INSERM, the CNRS 80|Prime interdisciplinary program, ANR and INCa. WW was supported by a COFUND IC-3i International PhD fellowship.

Materials Availability

No newly generated materials are associated with this paper.

Data and Code Availability

All raw and processed ORM data, including the identified IZs, is available at figshare <https://doi.org/10.6084/m9.figshare.c.5321048>. All ORM analysis code is available at GitHub <https://github.com/cl-chen-lab/ORM>, as are scripts to process ORM data for display on the ORM browser at <https://github.com/ClaireMarchal/ORM2BED> and <https://github.com/nimezhu/ORM>.

Experimental Model and Subject Details

This study used human HeLa S3 (RRID:CVCL_0058) and H9 cell lines (RRID:CVCL_9773). Both lines are derived from female donors. HeLa S3 cells were grown in DMEM plus 10% Cosmic Calf Serum (GE Life Science SH30087) and Pen/Strep. H9 hESCs were grown in feeder free conditions on Geltrex matrix (Thermo Fisher A14133) coated dishes in StemPro (Thermo Fisher A100701) media according to manufacturer's specifications.

Methods

Lead Contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Nick Rhind (nick.rhind@umassmed.edu).

windows by combing all 0' data, and removed those with >=20 ORM signals (the cutoff was estimated by comparing with the distribution of DLS samples without such hotspots). The removed hotspot signals correspond to 4.7% of total signals in C.0 sample and 3.6% in D.0 sample, and to similar degree (~3-4%) for the other samples from these datasets (i.e. C.5, C.10, D.20, D.30, D.45, D.60 and D.90).

Modeling Signal-Intensity and Inter-Signal-Distance Distributions

To model the signal-intensity distributions (Figure S2A), we assumed that signals represent individual fluorophores that can be resolved with some finite resolution. We modeled both the photon output and the inter-signal distance (which controls the probability of having two fluors closer together than the resolution of the Saphyr) of each fluor as Poisson processes, resulting in the depicted fit (obtained from the first four terms of Eq. 5, Supplemental Mathematical Methods), which predicts that 20% of the fluors are closer than the resolution of the Saphyr. See Supplemental Mathematical Methods for details.

To model the distribution of inter-signal distances (Figure S2B), we assume that the incorporation rate is a Poisson process for which the Poisson statistic is determined by the average amount of label taken up by cells. We further assume that there is an exponentially decreasing amount of label over time as it is consumed by replication. Combining these two assumptions results in the depicted fit (Eq. 14, Supplemental Mathematical Methods), with an inferred initial incorporation rate of about 1/900 thymidines and a depletion half length of about 75 kb for the merged 0-minute dataset. See Supplemental Mathematical Methods for details.

Combining the models in Figures S2A andB, we calculate that 20% of the signals in Figure S2B (corresponding to the 20% of the fluors that we infer from Figure S2A are not individually resolved by the Saphyr) are located within 1.3±0.3 kb of a neighboring fluor, consistent with the reported resolution of the Saphyr and the observation that there are almost no inter-signal distances below 1.3 kb in our datasets.

Identification of Replication Tracks

After mapping the fibers on the human genome, the distances between the genomic mapping signals were recalibrated to match their reference genome positions, compensating for heterogeneity in fiber stretching. We then mapped signals in the red channel between two adjacent genomic mapping signals as nucleotideincorporation signals. The nucleotide-incorporation signal is discontinuous and appears as clusters of neighboring signals (Figure 1A). To segment the ORM signals into replication-incorporation tracks, we used a two-step process based on Gaussian mixture models. We first fit a Gaussian mixture model using the R package mixtools (https://www.rdocumentation.org/packages/mixtools) with 3-gaussian distributions to the inter-signal distance distribution, with the first peak corresponding to the distances between signals within the same segment, the second peak corresponding to the distances between signals in adjacent segments and the third peak corresponding to signals in non-adjacent segments. The tail of the first Gaussian distribution (16,384 kb) was comparable amongst various data sets and was therefore selected as the cutoff value (cutoff1) to merge adjacent signals into primary segments. We then measured the distance between adjacent primary segments and fit a 2-Gaussian distribution to find the second cutoff value (cutoff2, 32,768 kb), which we used to merge all the adjacent primary segments into final replication tracks.

To estimate the resolution with which we can identify the initiation site of any given replication track, we defined the maximum likelihood distribution of initiation sites for a given replication track and, using numeric simulations, estimated the standard deviation of that distribution as 14.2±0.1 kb. See Supplemental Mathematical Methods for details.

To determine the direction of replication tracks in our asynchronous data, we defined a fork direction index (FDI): FDI = log2(IntensitySUM_Left / IntensitySUM_Right), where IntensitySUM_Left and IntensitySUM_Right is the total fluorescent signal in the left and right half, respectively, of each ORM track with at least 3 signals. A positive (negative) FDI corresponds to a rightward (leftward) moving

Quantification and Statistical Analysis

Statistical analyses were performed in R. Statistical tests and parameter are reported in the text and figure legends.

Additional Resources

Interactive visualization of the ORM data is available on the ORM Browser <http://orm.nucleome.org>.

Author Contributions

Conceptualization